*Article* **Statistical and Comparative Analysis of Multi-Channel Infrared Anomalies before Earthquakes in China and the Surrounding Area**

**Yingbo Yue 1,2, Fuchun Chen 1,\* and Guilin Chen <sup>1</sup>**


**\*** Correspondence: fuchun.chen@mail.sitp.ac.cn

**Abstract:** Abundant infrared remote sensing images and advanced information processing technologies are used to predict earthquakes. However, most studies only use single long-wave infrared data or its products, and the accuracy of prediction is not high enough. To solve this problem, this paper proposes a statistical method based on connected domain recognition to analyze multi-channel anomalies. We extract pre-seismic anomalies from multi-channel infrared remote sensing images using the relative power spectrum, then calculate positive predictive values, true positive rates and probability gains in different channels. The results show that the probability gain of the singlechannel prediction method is extremely low. The positive predictive value of four-channel anomalies is 41.94%, which is higher than that of single-channel anomalies with the same distance threshold of 200 km. The probability gain of the multi-channel method is 2.38, while that of the single-channel method using the data of any channel is no more than 1.26. This study shows the advantages of the multi-channel method to predict earthquakes and indicates that it is feasible to use multi-channel infrared remote sensing images to improve the accuracy of earthquake prediction.

**Keywords:** earthquake prediction; infrared remote sensing; multi-channel; pre-seismic anomaly; relative power spectrum; connected region

#### **1. Introduction**

Earthquake prediction is a complex and challenging theme. Since scientists discovered pre-seismic infrared abnormal phenomena in the 1980s [1], scholars around the world have undertaken various relevant research. Qiang et al. put forward a relatively reasonable theoretical mechanism, which indicates that the main causes of the abnormal temperature increase before the earthquake are gases released from the earth's crust and the change of the electric field, based on various observations and experiments [2]. The change of water content in the earth's surface soils is also able to cause infrared anomalies before earthquakes [3]. Some experiments proved that the infrared radiation of the rock changes when it is pressed by stress [4]. Their studies provide theoretical support for the development of earthquake prediction using infrared radiation data.

Infrared remote sensing images have a definite advantage of wide-field and continuous observation over ground-based observation and are widely applied to earthquake prediction [5–8]. Most researchers only analyzed single-type data, mainly using long-wave infrared radiation or its products, such as land surface temperature (LST) and outgoing long-wave radiation (OLR) [9–11]. In some case studies, the pre-seismic infrared anomaly is discovered from remote sensing images as a common phenomenon [12–14]. There are only a few studies about the middle-wave infrared anomalies before earthquakes, whose trend is similar to the anomalies in the long-wave channel [5,15]. Other earthquake case studies found that there are abnormal changes in water vapor content in the atmosphere before and after the event [16,17]. Most infrared data used to predict earthquakes are from

**Citation:** Yue, Y.; Chen, F.; Chen, G. Statistical and Comparative Analysis of Multi-Channel Infrared Anomalies before Earthquakes in China and the Surrounding Area. *Appl. Sci.* **2022**, *12*, 7958. https://doi.org/10.3390/ app12167958

Academic Editors: Eleftheria E. Papadimitriou and Alexey Zavyalov

Received: 8 July 2022 Accepted: 5 August 2022 Published: 9 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/).

NOAA satellites, Terra/Aqua satellites and Fengyun satellites [18,19]. Wei et al. found the infrared anomalies before the Ms8.0 earthquake in Sichuan, China, using single-channel images from the FY-2C satellite [20]. Ouzounov et al. found the anomalies before several strong earthquakes in Xinjiang, China, using OLR data from NOAA satellites and the FY-2D satellite [21]. Zhong et al. found the infrared anomaly associated with the 2017 M6.5 Jiuzhaigou earthquake from the data of two Fengyun satellites (FY-2E and FY-2G) [22].

Researchers have proposed different methods that were successfully used to extract pre-seismic anomalies, such as robust satellite techniques [23], interquartile, wavelet transform, Kalman filter methods [24], power spectrum [25] and some artificial intelligent methods [26]. Although abundant remote sensing data and advanced information processing technologies are used to predict earthquakes, there has not been any stable and valid algorithm to eliminate the influence of non-seismic factors, such as seasonal changes, weather conditions and human activities [18]. This is because of the complexity and variability of the earth system and the space environment. As a result, most methods merely work well in a few cases, lacking statistical evidence. Some statistical results in a small region show that there is an infrared anomaly before most earthquakes, but they did not analyze the proportion of the anomalies followed by an earthquake. In a study of 20 earthquakes in the Tibet region, the anomalies of the brightness temperature appeared before 17 earthquakes, and that of long-wave radiation appeared before 16 earthquakes [27]. The accuracy of earthquake prediction is unable to be estimated without PPV. Some statistical studies show a low positive predictive value (PPV). The statistical results based on the robust satellite techniques indicate the true positive rate (TPR) is high but the PPV of pre-seismic infrared anomalies, which is 25.9% in Sichuan province, is too low to put into practice [28]. The accuracy of earthquake prediction in the statistics by Jiao and Shan is 6.01% [29]. The PPV calculated by Filizzola et al. is 7.61% [23]. This means that most prediction results are wrong. Some studies do not even support the feasibility of earthquake prediction based on infrared remote sensing data [30]. Although some researchers obtained the PPV of 76.1% and the TPR of 67.1%, the spatiotemporal occupation of anomalies is high at 43.4%, which may cause a large prediction range and excessive public panic [31].

To explore a more valid algorithm for earthquake prediction using the data of infrared remote sensing, both data and research methods are improved. For one thing, multichannel infrared images are used in this paper, and the four channels could provide more information about the state of the earth's surface and atmosphere. For another thing, the statistical method based on connected region recognition is proposed to analyze the correlation between infrared anomalies and earthquakes, which could recognize spatiotemporal characteristics of anomalies in the long-term and wide-region studies. In this paper, the relative power spectrum method is respectively used to extract the anomalies from data of every channel. Both statistical analysis and case study are used to compare the prediction performance of the data from any single channel. Finally, four-channel anomalies are analyzed statistically. The results show that multi-channel anomalies could provide larger PPV and probability gain. It proves the potential of multi-channel infrared data in earthquake prediction and shows that it is possible to identify anomalies associated with earthquakes using multi-dimensional or multi-source data. The statistical method based on connected region recognition could be used to analyze pre-seismic anomalies from most kinds of remote sensing data. It lays the foundation for more data to be used in earthquake prediction.

#### **2. Materials and Methods**

The full research method is shown in Figure 1. Section 2.1 introduces the data and research range (time and region). The wavelet decomposition and relative power spectrum are used to extract anomalies from the time series on every pixel, which is detailed in Section 2.2. Then, the statistical method based on connected region recognition is proposed to analyze the correlation between anomalies and earthquakes, which is detailed in Section 2.3.

**Figure 1.** The flow of the research method.

#### *2.1. Data and Study Area*

Feng Yun-2G is the fifth operational satellite of the first generation of geostationary meteorological satellites launched by China on 31 December 2014, which can observe the parameters of the earth's surface and atmosphere widely and continuously. The main remote sensor is the visible and infrared spin scanning radiometer with 5 channels, including one visible channel (VIS), one medium-wave infrared channel (IR4), one water vapor channel (IR3) and two long-wave infrared channels (IR1 and IR2). The detailed parameters of each channel are shown in Table 1.

**Table 1.** The specific information of each channel.


In this paper, the original data are the images from four infrared channels of the Feng Yun-2G satellite, which are provided by the National Satellite Meteorological Center. (http://satellite.nsmc.org.cn/PortalSite/Data/Satellite.aspx, accessed on 25 June 2022). One full-disk data file is generated every half hour or one hour. The file structure is shown in Figure 2, including file attributes and scientific datasets. Every file includes the geographic location of the image center, cloud classification, five-channel images and calibration tables and is saved as the hierarchical data format (HDF). The red digits are the number of the datasets, while the white digits are the size of one dataset. Moreover, the center also provides a lookup table so that the users could calculate the longitude and latitude of every pixel in the image.

Due to the high transmittance of the earth's atmosphere in the atmospheric windows, the radiation value from two long wave channels and one middle wave channel mainly depends on the temperature and the emissivity of the earth's surface. The emissivity of the same target is usually constant. The peak wavelength of most objects on the earth's surface is in the long-wave infrared channel, so the brightness temperature data from the long-wave infrared channel could show the temperature change trend of most objects on the earth's surface [32]. The data from the medium-wave infrared channel show the temperature change of the high-temperature targets and are affected by the reflection of solar radiation in the daytime, so data at night are selected in this paper. The water vapor channel is in the strong absorption band of water vapor, which is one of the main infrared absorption gases in the earth's atmosphere [33].


**Figure 2.** The structure of the data file.

China is located between the Pacific seismic zone and Euro–Asia seismic zone, so earthquakes happen frequently there, especially in Qinghai–Tibet Plateau, Yunnan–Guizhou Plateau and Taiwan. In China, the terrain is high in the west and low in the east, as shown in Figure 3. The lines in different colors are the different types of plate boundaries. The topography data are provided by the University of Califonia San Diego. (http: //topex.ucsd.edu/marine\_topo/mar\_topo.html, accessed on 23 July 2022). China has a wide territory and rich soil resources, mainly including 15 types of soil. The variation of water content in the soil during the earthquake preparation period may also lead to anomalies in infrared remote sensing images [3]. The study area of this paper is between 0◦ N to 60◦ N and 70◦ E to 140◦ E, including China and the surrounding area. All images in this region from June 2015 to December 2020 are applied to extract abnormal signals. There are 358 earthquakes with a magnitude over five in the area from 31 January 2016 to 1 January 2021. The strongest earthquake in the study range was the earthquake with a magnitude of 7.3 that happened in Kyushu, Japan on 16 April 2016. According to the theoretical model proposed by Dobrovolsky, the relationship between the radius of earthquake preparation region R and the earthquake magnitude M is shown in Equation (1) [34].

$$R = 10^{0.43\text{M}}\text{ },\tag{1}$$

The earthquake magnitude in this study is between 5 and 7.3, so the radius of the influenced region varies from 141 to 1377 km.

**Figure 3.** The terrain in the study area.

#### *2.2. Anomaly Extraction*

The relative power spectrum is a common method to extract seismic information from infrared remote sensing images and is used in some case studies with significant results. The process is shown in Figure 4 and was programmed using Python in the Spyder platform. The algorithm is introduced in detail as follows [35,36]:


**Figure 4.** Flow chart of anomaly extraction algorithm.

In this paper, the algorithm is used to process the data from every channel of FY-2G, respectively, to compare the difference between channels and improve prediction performance using multi-channel data.

#### *2.3. Statistical Method*

Accuracy of earthquake prediction is the basis of practical application, because wrong predictions may cause serious economic losses and public panic. The PPV of infrared abnormal signals can evaluate the accuracy of its application in earthquake prediction, which is defined as the proportion of abnormal signals related to an earthquake in total anomalies, as shown in Equation (2) [28].

$$\text{PPV} = \frac{\text{N}\_0}{\text{N}\_a} \tag{2}$$

where *N*<sup>0</sup> is the number of abnormal signals related to an earthquake, and *Na* is the number of total abnormal signals in the research area. To calculate the PPV, we propose a method based on the connected domain to identify the spatiotemporal occupancy range of abnormal signals. We deem the point where the relative power spectrum is more than five as abnormal. The connected region of abnormal points is regarded as the spatial–temporal range of a single abnormal phenomenon, and abnormal regions with relatively small areas or short-term are ignored to eliminate the infrared radiation increase caused by human activities. In this paper, an anomaly that covers less than 500 pixels or lasts less than 3 days is ignored. In previous studies, anomalies appeared on average one to three weeks before the earthquake [38]. As a result, the anomaly is regarded to be related to the earthquake, if the distance between the abnormal region and the epicenter is less than the distance threshold (Td) and the anomaly appears within one month before the earthquake.

The TPR is the proportion of earthquakes with pre-seismic infrared abnormalities in total earthquakes, which is defined as Equation (3) [28].

$$\text{TPR} = \frac{N\_1}{N\_c} \tag{3}$$

where *N*<sup>1</sup> is the number of earthquakes with any pre-seismic anomaly, and *Ne* is the number of total earthquakes in the research area. It is used to assess the universality of infrared anomalies before earthquakes.

It provides higher PPV and TPR to increase the distance threshold, but the spatial accuracy of prediction would decrease. Therefore, earthquake prediction needs to predict more earthquakes successfully with lower spatiotemporal occupancy. The probability gain is the ratio of TPR to spatiotemporal occupancy, shown as Equation (4) [39,40]. It can be regarded as a criterion for selecting the optimal distance threshold.

$$\text{Gain} = \frac{\text{TPR}}{\pi} \text{.} \tag{4}$$

where TPR is the true positive rate, τ is the fraction of space–time occupied by the predicted range, and is associated with the distance threshold. Matlab platform was used for statistical analysis and results display.

#### **3. Results and Discussion**

#### *3.1. Channels Comparison*

Some previous research found the epicenter of an impending earthquake may be far from the anomaly [41–45], so both the abnormal region and the region around the anomaly should be considered as the predicted region. The possibility that there is any upcoming earthquake in the predicted region is larger if the area of the predicted region is larger. PPVs with different distance thresholds in four channels are calculated and shown in Figure 5. There are a total of 619 anomalies from 2016 to 2020, including 155 in the IR1 channel, 161

in the IR2 channel, 132 in the IR3 channel, and 171 in the IR4 channel. PPVs vary obviously with the change of the distance threshold. The difference between the PPVs in any two channels is less than 0.08 for the same distance threshold. It means that there is not much difference between the PPVs of anomalies from different channels.

**Figure 5.** Positive predictive values in different channels at different distance thresholds.

Figure 6 shows the TPRs with different distance thresholds in four channels. The TPRs in the two long-wave infrared channels are similar. The TPR in the IR4 channel is the highest among all channels whatever the distance threshold is, while that in the IR1 channel is the lowest. The difference between the TPRs in the IR2 channel and the IR3 channel becomes large with the increase in the distance threshold. The long-wave infrared channels (IR1 and IR2) behave worse on PPV and TPR, although they were earlier used for earthquake prediction.

**Figure 6.** True positive rates in different channels at different distance thresholds.

The infrared anomalies of 358 earthquakes in the study area were analyzed statistically. As shown in Figure 7, the horizontal axis is a four-digit number whose four digits represent the state of four channels (IR1, IR2, IR3, IR4) successively, where 1 means there was an anomaly before the earthquake, and 0 means there was no anomaly before the earthquake. For example, 1111 means there were abnormalities in all four channels. The state of the anomalies in the two long-wave channels is similar. There were only 26 earthquakes with four-channel anomalies within 400 km around the epicenter, and 18 of them were within 200 km.

As shown in Figure 8, the Gain in each channel varies slightly with the distance threshold. The probability gain cannot be improved by changing the distance threshold. The gains in the IR3 channel are highest among those in all four channels for the same distance threshold and only a little higher than that in other channels.

**Figure 7.** Infrared anomalies before earthquakes.

**Figure 8.** Probability gains in different channels at different distance thresholds.

#### *3.2. Case Study*

Considering regional differences, three seismic cases in different provinces were taken as examples. Detailed seismic information is shown in Table 2. The magnitude of the Sichuan earthquake was the largest. All three earthquakes occurred at shallow depths.



#### 3.2.1. Sichuan Earthquake

There was a strong earthquake with a magnitude of 7.0 in Sichuan at 21:19:46 on 8 August 2017 (LT). The epicenter was at 103.82◦ E, 33.20◦ N. The depth was 20 km. At 07:27:52 on the next day (LT), an earthquake with a magnitude of 6.6 occurred in Xinjiang. The epicenter was at 82.89◦ E, 44.27◦ N. The depth was 11 km.

The information on pre-seismic anomalies in different channels is shown in Table 3. In Table 3, the start time and the end time are the numbers of days that the start time and end time of the anomaly relative to the earthquake time. Negative numbers are before the earthquake and positive numbers are after the earthquake. The anomalies of the two long-wave channels appeared 25 days before the earthquake, as shown in Figure 9. On that day, there were discrete high-value points north of the epicenter in the images of the four channels, and the anomaly area of the long-wave infrared channels was larger, so it was identified earlier. As shown in Figure 10, there were high-amplitude and large-area anomalies in four channels on the 19 days before the earthquake. On that day, IR1, IR2 and IR4 obtained their maximum abnormal value. The anomalies were located between the epicenters of the two earthquakes and more close to that of the Sichuan earthquake. Figure 11 shows the largest coverage area of anomalies in the earthquake preparation region. During the abnormal period, the abnormal area of the four channels went through two times of increase and decrease. The abnormal area in the IR3 channel was less than 500 pixels 19 days before the earthquake, while that in the IR4 channel was also too small during the first time of increase and decrease. They are both ignored in Table 3. Anomalies in all channels lasted until 18 days before the earthquake.

**Table 3.** The spatial–temporal characteristics of anomalies before the Sichuan earthquake.


**Figure 9.** Power spectrum images on 14 July 2017.

**Figure 10.** Power spectrum images on 20 July 2017.

**Figure 11.** The cover area of the anomalies before and after the Sichuan earthquake.

#### 3.2.2. Taiwan Earthquake

An earthquake with a magnitude of 5.2 happened in Taiwan on 16 December 2018, 05:21:05 (LT). The location of the epicenter was 23.71◦ N, 121.80◦ E, and the focal depth was 26 km. On the same day, an earthquake with a magnitude of 5.7 broke out at 12:46:07 (LT) in Sichuan. The location of the epicenter was 28.24◦ N, 104.95◦ E, and the focal depth was 12 km.

As shown in Table 4, the IR3 channel appeared abnormal at first, and anomalies in other channels appeared simultaneously and disappeared simultaneously, which were a little later than that in the IR3 channel. In images from the IR3 channel, a high-value region arose on 6 December 2018, shown in Figure 12. In the IR3 channel image, there was a weak anomaly between the two epicenters. Relatively, there were more obvious anomalies in the high latitude region, but they were far from the two epicenters, which is not significant for the prediction of this earthquake. It became an isolated region on 8 December 2018, while the anomalies in the other three channels just started, as shown in Figure 13. The anomalous positions of the two long-wave channels and the medium-wave channel were close, but the anomaly area of the medium-wave channel was smaller. The anomaly area of the IR3 channel was the largest, and it has obvious location deviation from the other three channels. On 9 December 2018, the anomalies in IR1 and IR3 obtained their maximums, as shown in Figure 14. In the abnormal region of the IR3 channel image, the southern part is closer to the epicenter of Taiwan, in which anomaly intensity was significantly higher

than that in the northern part. Similar to the first case, the anomaly occurred between the epicenters of the two earthquakes. The abnormal region was mainly on the land and distributed along the coastline. The epicenter of the Taiwan earthquake was at the edge of the anomaly. The anomaly in the IR3 channel was closest to the epicenter. There was only 1.7 km between the anomaly and the epicenter. The anomalies in four channels all disappeared completely on 12 December 2018. According to Figure 15, the anomalies went through appearance, increase, decrease and disappearance before the earthquake. The anomaly in the IR4 channel intersected with that in other channels, but the distance from the epicenter was far beyond the radius of the earthquake preparation region calculated theoretically.

**Figure 12.** Power spectrum images on 6 December 2018.

**Figure 13.** Power spectrum images on 8 December 2018.

**Figure 14.** Power spectrum images on 9 December 2018.

**Figure 15.** The cover area of the anomalies before and after Taiwan earthquake.



#### 3.2.3. Tibet Earthquake

Tibet is located on the Qinghai–Tibet Plateau with high altitudes and complex terrain. There was an earthquake with a magnitude of 5.6 at 17:22:14 on 19 July 2019. The epicenter was located at 27.67◦ N, 92.89◦ E. The depth was 10 km. Anomalies before the earthquake are shown in Table 5. From 30 days to 20 days before the earthquake, only the IR4 channel showed abnormalities. Anomalies in the IR3 and IR4 channels appeared on the 11th day before the earthquake. IR1 and IR2 showed abnormalities later. The occurrence time and disappearance time of the anomalies in the four channels were relatively close. The relative powers spectrum on 9 July 2019 is shown in Figure 16. Anomalies in all four channels were close in location. The distance between the anomalies and the epicenter is 11.5 km. The

cover area of the anomalies before and after the earthquake is shown in Figure 17. The variation trend of the cover area of anomalies in the four channels is similar during the four-channel abnormal period. For this seismic case, the anomalies of the four channels show obvious similarities in spatiotemporal characteristics and the evolution process.


**Table 5.** The spatial–temporal characteristics of anomalies before the Tibet earthquake.

**Figure 16.** Power spectrum images on 7 July 2019.

**Figure 17.** The cover area of the anomalies before and after Tibet earthquake.

#### *3.3. Multi-Channel Anomalies*

To improve the accuracy of earthquake prediction, we scan the relative power spectrum data in the four channels and identify the region of four-channel anomalies, which means the region where all of the four channels are abnormal at the same time. The statistical result is shown in Table 6. A larger distance threshold can improve PPV and TPR, but can not improve the probability gain. The probability gain is the highest when the distance threshold is 200 km. Figure 18 shows the comparison between the PPV, TPR and gain of multi-channel methods with that of the single-channel method with the distance threshold of 200 km. The PPV and probability gain of the multi-channel method are higher than that of any single-channel method. In particular, the probability gain is roughly doubled, although the TPR is very small. Only 19 earthquakes are associated with the multi-channel anomalies because the multi-channel method only considers simultaneous anomalies of four channels, and earthquakes with anomalies of no more than three channels will be missed. The locations of the epicenters of these 19 earthquakes are shown in Figure 19. Four-channel infrared anomalies were found before the earthquakes in different regions, which indicates the phenomenon is not only suitable for a specific region. The distribution of anomaly time is shown in Figure 20. Eight earthquakes occurred one to two weeks after the anomalies appeared. Only one earthquake occurred three to four weeks after the anomaly appeared. Figure 21 is the Molchan diagram of the different methods. The spatiotemporal occupancy of the four-channel anomaly is much smaller than that of the single-channel anomaly. Although the miss rate is so high, enough small space–time occupations could provide a high value of probability gain.

**Table 6.** Statistical result of multi-channel anomalies.


**Figure 18.** Comparison of single-channel method and multi-channel method.

**Figure 19.** Epicentres of earthquakes following four-channel anomaly.

**Figure 20.** The number of days between anomalies and earthquakes.

**Figure 21.** Molchan diagram of different methods.

To compare the multi-channel method with other methods, the authors, data, study region, study period, earthquake magnitude and their results in previous studies are shown in Table 7 [23,28,29,31,46]. The anomaly time means the number of days from anomaly to earthquake. The multi-channel method could obtain the highest probability gain, although the TPR is the lowest. This means that this method reduces the uncertainty of earthquake prediction more obviously than other methods. The PPV of the multi-channel method is also the highest in the studies on earthquakes of magnitude five and above. It means that the multi-channel method has higher accuracy in earthquake prediction than previous methods for earthquakes with magnitude five or above.



1. OLR means outgoing long-wave radiation; 2. ST means land surface temperature; 3. TIR means thermal infrared radiation; 4. BT means brightness temperature; 5. None means that the parameter was not mentioned in the study; 6. ≈ means the value was estimated according to the point in the images.

#### **4. Conclusions**

Based on the relative power spectrum method, we propose a statistical method based on connected domain recognition to calculate the PPVs, TPRs and probability gains in different channels. The results show that the PPV and TPR could be improved by increasing the distance threshold. The probability gain is low and its change with distance threshold is not obvious. In addition, we also statistically analyzed the multi-channel infrared anomalies before 358 earthquakes. There is at least one channel anomaly within one kilometer of the epicenter within one month before 36.87% of the earthquakes, but there are only 26 earthquakes with four-channel anomalies within 400 km of the epicenter and 18 earthquakes with four-channel anomalies within 200 km of the epicenter.

In the study of three earthquake cases, four-channel anomalies appeared and disappeared before the earthquake. The epicenter is at or some distance from the edge of the anomaly. Due to the low PPV and probability gain of the earthquake prediction method using single-channel data, multi-channel infrared remote sensing images are used for earthquake prediction. The PPV of four-channel anomalies is 41.94%. This is higher than that of single-channel anomalies at the same distance threshold of 200 km. Meanwhile, the method causes a lower TPR. Significantly, the spatial–temporal occupancy of four-channel anomalies is very low, and the probability gain is doubled.

This study shows the difference between pre-earthquake anomalies in multi-channel infrared remote sensing images and indicates that multi-channel infrared remote sensing images may have more advantages in the PPV and the probability gain of earthquake prediction than single-channel data. In earthquake prediction, the PPV, which indicates the reliability of the algorithm, is more important than the TPR. This is because the accurate prediction of a single earthquake can also save a lot of life and property. However, it is still difficult to use the four-channel infrared data to obtain high enough accuracy for the practical application of earthquake prediction. The purpose of this paper is to show the advantages of multi-channel data over single-channel data. The results could be compared with other types of pre-seismic anomalies to study the mechanism of anomalies during the earthquake preparation period and explore a better method for earthquake prediction.

Earthquake prediction needs a lot of remote sensing data and ground-based observation data. In future studies, we can improve the performance of earthquake prediction with the following aspects:


The anomaly extraction algorithm in this paper can be used for time series analysis of other remote sensing data, and the statistical method can be used for other types of wide-field and long-time spatial–temporal data. Some parameters may need to be adjusted for using other data.

**Author Contributions:** Conceptualization, G.C.; methodology, F.C.; investigation, Y.Y.; writing original draft preparation, Y.Y.; writing—review and editing, F.C.; supervision, G.C.; funding acquisition, F.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** Thanks to the National Satellite Meteorological Center of China for providing the data of the FY-2G satellite.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Maksim Gapeev** *<sup>∗</sup>***,† and Yuri Marapulets †**

Institute of Cosmophysical Research and Radio Wave Propagation FEB RAS, Kamchatka Region, Elizovskiy District, Mirnaya Str. 7, 684034 Paratunka, Russia

**\*** Correspondence: owlptr@yandex.ru

† These authors contributed equally to this work.

**Abstract:** In seismically active regions of the Earth, to which the Kamchatka peninsula refers, preseismic anomalies are recorded in different geophysical fields. One of such fields is the acoustic emission of rocks, the anomalies of which are recorded 1–3 days before earthquakes at the distance of the first hundreds of kilometers from their epicenters. Results of joint acoustic-deformation measurements showed that growth of geoacoustic radiation intensity occurs during the increase in the level of deformations in rock masses by more than one order compared to the background values. Simulation studies of the areas with increased deformation are realized to understand the causes of anomalous acoustic-deformation disturbance occurrences before strong earthquakes. The model is based on the assumption that the Earth's crust in the first approximation can be considered as a homogeneous isotropic elastic half-space, and an earthquake source can be considered as a displacements along a rectangular fault plane. Based on these assumptions, deformation regions of Earth's crust were modeled during the preparations of two earthquakes with local magnitudes *ML* ≈ 5 occurred on the Kamchatka Peninsula in 2007 and 2009. The simulation results were compared for the first time with the data of a laser strainmeter-interferometer installed at the Karymshina observation site (52.83◦ N, 158.13◦ E). It was shown that, during the preparation of the both earthquakes, the Karymshina observation site was within the region of shear deformations <sup>≈</sup>10<sup>−</sup>7, which exceeded the tidal ones by an order. On the whole, simulation results corresponded to the results of the natural observations. Construction of an adequate model for the generation of acoustic-deformation disturbances before strong earthquakes is topical for the development of an early notification system on the threat of catastrophic natural events.

**Keywords:** earthquake preparation; areas of increased shear deformations; mathematical simulation; rock deformation; acoustic emission of near-surface rocks

#### **1. Introduction**

It is generally accepted that mechanical processes play a leading role in the preparation of seismic events [1–3]. They cause increased stresses, leading to deformations of the Earth's crust around an earthquake source. These changes in the stress–strain state of rocks lead to the anomaly occurrences, classified as earthquake precursors, in various geophysical fields [4–9].

The increase in rock acoustic emission in kilohertz frequency range is one of the identified pre-seismic anomalous disturbances in geophysical fields. Such anomalies were observed in various seismically active regions of the world: in Armenia [10], in Italy [11], and in Russia on the Kamchatka peninsula, which is a part of the Circum-Pacific Orogenic Belt, also known as the "Ring of Fire" [12–14]. As a result of long-term studies of acoustic emission in Kamchatka, a high-frequency acoustic emission effect was revealed [15]. This effect consists in an increase of geoacoustic radiation intensity with an increase of rock deformation rate. This effect is determined by rock deformations at observation sites and

**Citation:** Gapeev, M.; Marapulets, Y. Modeling Locations with Enhanced Earth's Crust Deformation during Earthquake Preparation near the Kamchatka Peninsula. *Appl. Sci.* **2023**, *13*, 290. https://doi.org/10.3390/ app13010290

Academic Editor: Marco Vona

Received: 8 December 2022 Revised: 22 December 2022 Accepted: 22 December 2022 Published: 26 December 2022

**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/).

manifests the most clearly in the kilohertz frequency range 1–3 days before earthquakes at a distance of the first hundreds of kilometers from their epicenters [16].

The peculiarity of the geoacoustic observations in Kamchatka is the application of broadband piezoceramic hydrophones installed in natural and artificial reservoirs to record the signals. The use of receivers of this type allowed us to expand the frequency range of registration to 0.1 Hz–11 kHz in comparison with standard geophones [16]. To confirm the deformation nature of acoustic emission anomalies, a laser strainmeter-interferometer of unequal shoulder type with a measuring base length of 18 m and sensitivity of 10–11 m was installed in 2005 at the Karymshina site (52.83◦ N, 158.13◦ E), located 41 km southwest of Petropavlovsk–Kamchatskiy [17]. Taking into account the peculiarities of its installation on the surface without optical waveguide, the calculated measurement accuracy of relative deformations was no less than 10−<sup>8</sup> [18]. The data of joint acoustic-deformation observations at this site was studied in detail [13]. It is shown that with the deformation intensification at the observation site, the near-surface sedimentary rock deformation rate also increases simultaneously. These rocks have a polydisperse fluid-saturated porous structure of low strength. Such activation of deformations is accompanied by relative micro-displacements of rock fragments, their interaction and, as a consequence, generation of acoustic emission of increased intensity. Such effects are observed the most clearly at the final stage of earthquake preparation a few days before their onset. Two cases of high-frequency geoacoustic responses to the activation of deformation process before two earthquakes with local magnitudes *ML* ≈ 5 were detected and studied in detail [13]. These earthquakes occurred in 2007 and 2009 at the epicentral distances of about 150 km to the Karymshina site. It was reasonable to make model studies of the deformation fields that occurred during the seismic event preparation and to correlate them with the real deformation levels recorded at the observation site. Comparison of the results of the simulation and the natural experiment is topical to understand the causes of anomalous acoustic-deformation disturbances at the final stage of earthquake preparation.

#### **2. Research Significance**

Investigation of pre-seismic anomalies in geophysical fields is of high practical significance for the development of the methods of early notification on seismic hazard. Geophysical field intensity variations in seismically active regions are known to be associated with stress–strain state near earthquake sources during their preparation. In this case, pre-seismic anomalies are often recorded at the distances of hundreds of kilometers from preparing earthquake epicenters. That corresponds to the cases of acoustic emission of rock anomalies recorded in kilohertz frequency range. However, signals at such frequencies cannot propagate from preparing earthquake epicenters due to strong attenuation. Thus, they occur as the result of medium response at the recording site on the change of its stress-strain state. It appears that the changes of stress–strain state near earthquake sources propagate at the distances up to hundreds of kilometers. Undoubtedly, that requires theoretical, model and experimental confirmation. Such investigations of the fields of the Earth's crust stress and deformation are topical in order to obtain new knowledge on the physical processes occurring during earthquake source preparation.

The basic concepts of plate tectonics, theory of elasticity and rock plastic deformations were described by D. Turcotte and G. Shubert [19]. V. Nikolaevsky [20] formulated the basic concepts of deformation and destruction of fractured rocks under static and dynamic load. C. Sholz [1] described in detail the modern understanding of earthquake mechanics. The theory of earthquake-induced seismic wave propagation and modern method of observation data interpretation and processing were considered in detail by K. Aki and P. Richards [21]. Analytical solutions of the mathematical model, describing the Earth's crust deformations in medium elastic approximation under earthquake source effect, were considered in detail by Yu. Okada [22,23]. I. Dobrovolskiy [24] introduced the notion of precursor manifestation zone and showed that its size is determined by the energy of a preparing earthquake. I. Dobrovolskiy [25] also suggested limiting the dimensions

of this zone by the boundaries, behind which anomalies in deformations do not exceed the deformation process background values of the order of 10−8. A. Alekseev and his co-authors [26] connected the geophysical field anomalies, occurring in seismically active regions, with the appearance of zones of geo-environment nonlinear loosening (dilatancy). Dilatancy zones, which are formed in the vicinity of earthquake sources during the stresses close to destructive values for rocks, were modeled in the approximation of elastic halfspace [26]. Three-dimensional block visco-elastic model of the south-eastern Tibetan Plateau was constructed [27]. Earth's crust model taking into account friable-plastic transitions in the Earth's crust was constructed. This model is based on Maxwell and Kelvin–Voigt models [28]. Visco-elastic fault-based model of crustal deformation was proposed for the 2023 Update to the "U.S. National Seismic Hazard Model" [29]. Earth's crust deformations model under earthquake source impact in the form of a rectangular plane was developed. In this model the Earth is a homogeneous elastic sphere [30]. Some complication of this model was made. It was generalized to the case of a layered elastic sphere [31]. The pre-seismic deformations before the Japanese earthquake occurred on 11 March 2011 in the Tohoku region were estimated [32]. An analysis of seismic anomalies associated with Ludian earthquake on 3 August 2014 was made [33]. Post-seismic deformations after Kokoxili earthquake occured on 14 November 2001 in the Northern Tibetan Plateau were estimated using satellite data [34]. A wide review of articles on the estimates of the Earth's crust deformations based on satellite data was made in 2018 [35].

Simulation of stress fields formed around Kamchatka earthquake sources was carried out earlier. The authors of the studies used different approaches to describe an earthquake source as a point source in the form of a single force [36] or a double force [37,38], as a distributed source in the form of a rectangular plane [39]. It was shown in all the cases that regions with deformations, exceeding the background values, occur around earthquake sources at the distances of hundreds of kilometers. Comparative modeling of deformation fields using these source models was carried out. It was shown that it is better to use the models of a distributed source in the form of a rectangular plane, as they describe more accurately the force action in the earthquake source [40].

The proposed paper continues those investigations. Based on the assumption that the Earth's crust can be considered as a uniform isotropic elastic half-space in first approximation and an earthquake source is a shift along a fracture rectangular plane, the authors modeled Earth's crust deformations, occurring around sources during the preparation of two earthquakes in Kamchatka. For the first time ever, the deformation values, obtained during the modeling, were compared with the data from a laser strainmeter-interferometer installed at the distance of several hundreds of kilometers from the earthquake sources.

#### **3. Research Methods**

The source of a tectonic earthquake is formed as a result of release of the stresses accumulated by elastic medium during tectonic deformation [21]. As a result of this release, a break of medium continuity appears. The accumulated elastic energy of deformation turns inelastic. According to this theory, an earthquake source can be described through a displacement along a fault plane [21]. It is notable that a displacement along a fault excites the same seismic waves as some system of forces distributed on the fault with zero total moment. In general, distribution of forces may have different form. However, in the case of an isotropic medium, it can always be chosen as a surface distribution of double pairs of forces.

In accordance with that, limiting ourselves to the consideration of a fault flat plane, the earthquake source model can be represented schematically as follows (Figure 1).

**Figure 1.** Schematic description of the earthquake source model. In the figure: *α* is the dip angle, *β* is the strike angle, *δ* is the angle of displacement direction, *C* is the hypocenter depth, *L* is the plane length, *W* is its width, *N* is the North direction (the axis is aligned with *OY*), and Σ is the fault plane with an equivalent distributed system of double forces with a moment.

Some parameters of this model (*α*, *β*, *δ*, *C*) can be accessed directly, for example, from the Harvard Catalog of Earthquake Mechanics CMT Catalog [41].

The linear dimensions of the fault plane, *L* (km) and *W* (km), as well as the displacement *U* (cm) magnitude can be estimated using the following correlation equations [42]:

$$\begin{aligned} \lg(L) &= 0.75 \cdot M\_W - 3.60, \\ \lg(W) &= 0.75 \cdot M\_W - 1.45, \\ \lg(L) &= 0.75 \cdot M\_W - 0.37, \end{aligned} \tag{1}$$

where *MW* = 2/3(*M*<sup>0</sup> − 16.1) is the moment magnitude, *M*<sup>0</sup> is the scalar seismic moment.

The research objective is the simulation of stress and strain fields caused by energy accumulation during earthquake preparation. It is obvious that this energy is significantly greater than the released energy of elastic deformations at the times of earthquakes.

In a generalized form, the correlation Equations (1) are presented as follows:

$$\log(N\_E) = a \cdot M\_W + b\_\prime \tag{2}$$

where *a* and *b* are some coefficients, *NE* is the characteristic of an earthquake source, calculated taking into account the released energy of elastic deformations.

The efficiency coefficient of elastic deformation energy release is:

$$
\eta = \frac{E}{W'} \tag{3}
$$

where *W* is the total energy of elastic deformations in the area including the earthquake source before the fault activation.

In Equation (2), the moment magnitude is expressed in terms of the earthquake energy *E*, using the Gutenberg–Richter equation: *E* = 101.5*MW* <sup>+</sup>5. Eliminating the logarithm, the following relation is obtained:

$$N\_E = \left(\frac{E}{10^5}\right)^{(2/3)a} \cdot 10^b. \tag{4}$$

In Equation (4), the earthquake energy *E* is replaced by the total energy of elastic deformations *W*, using Equation (3). The following relation is obtained:

$$N\_W = \left(\frac{1}{\eta}\right)^{(2/3)a} \cdot N\_{E\_{\prime}} \tag{5}$$

where *NW* is the earthquake source characteristic calculated taking into account the total energy of elastic deformations. The coefficient (1/*η*)(2/3)*<sup>a</sup>* carries the meaning of an increasing coefficient for correlation Equation (1). This coefficient makes it possible to calculate the stress–strain state of rocks taking into account the total energy of elastic deformations during earthquake preparation.

There are various approaches to evaluate both effective released stress [43,44] and to evaluate the efficiency of elastic deformation energy release. For example, the most accurate approach to estimate the coefficient of efficiency of elastic deformation energy release *η*, which requires reconstruction of tectonic stress in a seismically active region, is described by Yu. Rebetsky [45]. I. Dobrovolsky [24] proposed a less accurate but a simpler variant to calculate the coefficient:

$$
\eta = 10^{0.26M\_W} - 3.93.\tag{6}
$$

Equation (6) was used in further calculations.

#### **4. Simulation of Stress and Strain Fields**

The following model for the formation of regions with increased deformation of the Earth's crust during earthquakes preparation is proposed. The Earth's crust is a homogeneous isotropic elastic half-space. The model of an earthquake source is a dislocation in the form of a rectangular plane with a constant displacement vector (Figure 1). The stress–strain state of the Earth's crust is determined by the accumulated elastic energy in the process of earthquake preparation. Zones of acoustic-deformation anomalies are the areas of daytime surface defined by the equation *z* = 0 with the level of relative deformations exceeding the tidal ones (>10−8). Shear sources of acoustic emission prevail, since rock strength with respect to tangential stresses is less than to compression. Therefore, only shear deformations are taken into account in the simulation.

Using Mindlin's solutions [46,47], Yu. Okada [22,23] obtained compact analytical solutions for the displacement vector and its spatial derivatives in the case of three types of displacement: in the direction of strike, in the direction of dip and expansion.

The Navier equations underlying the model are linear. Therefore, the solution in case of an arbitrarily oriented displacement (not for expansion) can be obtained in the form of a linear combination of solutions for the displacement in the strike and dip directions:

$$\mathcal{U}\_{\rm strike} = \mathcal{U} \cdot \cos(\delta),\tag{7}$$

$$\mathcal{U}\_{dip} = \mathcal{U} \cdot \sin(\delta),\tag{8}$$

where *Ustrike* is the component of the displacement vector along the strike, *Udip* is the component of the displacement vector along the dip, *δ* is the angle of displacement direction.

The Mercator projection was used to convert the geographical coordinates to Cartesian ones. An additional coordinate system was built for each earthquake to simplify the calculations. A system had a center at the earthquake epicenter and was oriented relative to the OZ axis of the original system by the angle of *β* − 90◦. Thus, the axes *OX* and *OY* were parallel to the projections of the sides *L* and *W* of the displacement plane on the plane *z* = 0 (Figure 2).

**Figure 2.** Schematic representation of an additional coordinate system centered at the earthquake epicenter (asterisk). The dotted line is the projection of the fault plane onto the Earth's surface (*z* = 0).

Two earthquakes, before which simultaneous anomalies of acoustic emission and rock deformations were observed at the Karymshina site, were simulated [13]. Earthquake No. 1 occurred on 2 May 2007, at 12:00:48.4 UT, the coordinates are 52.29◦ N, 160.55◦ E, the depth is 28 km, local magnitude *ML* = 5.2. Earthquake No. 2 occurred on 8 October 2009, at 05:25:13.4 UT, the coordinates are 52.84◦ N, 160.15◦ E, the depth is 20 km, local magnitude *ML* = 5.1. These earthquakes are not listed in the CMT catalog due to their low energy. The data were taken from the earthquakes catalog for Kamchatka and the Commander Islands [48]. Unfortunately, it is impossible to obtain the information on the orientation of the displacement plane from this catalog. This information is necessary for the computational experiment. To estimate it, the earthquakes, which occurred in the area of the earthquakes under the study and which were presented in the CMT catalog, were analyzed. The analysis was carried out under the assumption that for the Kamchatka subduction zone, the general directions of force impact from the sources of small-focus earthquakes, located in some small area, are quite constant. For this purpose, the entire observation period from 1 January 1976 to 30 June 2022, presented in the CMT catalog, was considered. Seismic events in the area, close to the simulated earthquakes (latitude interval: [52◦, 53◦], longitude interval: [160◦, 161◦]), were analyzed.

In total, nine small-focus earthquakes were represented in the CMT catalog in this area. More detailed information about these earthquakes, including information about the focal mechanism, is presented in Table 1. Epicenter locations are shown in Figure 3. Numbers of earthquakes in Figure 3 correspond to the ones in Table 1.


**Table 1.** Data on nine small-focus earthquakes that occurred during the period from 1 January 1976 to 30 June 2022 in the area under the study (latitude interval: [52◦, 53◦], longitude interval: [160◦, 161◦]).

<sup>1</sup> Degree measures of the angle are given.

**Figure 3.** Map of Kamchatka peninsula with location of earthquake epicenters presented in the CMT catalog and occurred during the period from 1 January 1976 to 30 June 2022 (black circles) in the area under the study and location of simulated earthquake epicenters (red squares). Location of earthquake epicenters are shown in scaled part of map. The black triangle on the map indicates the location of the Karymshina observation site.

Only one of them (Earthquake No. 4 in Table 1) significantly differed in the orientation of the fault plane. All other earthquakes were very similar in these parameters. It was removed from the sample and the following statistical estimates of the orientation angles were obtained:

$$
\mathbb{R} = 211.25^{\circ}, \; \mathbb{S}(\mathfrak{a}) = 8.48^{\circ}, \tag{9}
$$

$$
\bar{\beta} = 26.25^{\circ}, \; \mathcal{S}(\beta) = 6.55^{\circ}, \tag{10}
$$

$$
\vec{\delta} = 83.38^{\circ}, \; S(\delta) = 12.08^{\circ}, \tag{11}
$$

where *α*¯, *β*¯, ¯ *δ* are average values of angles, *S*(*α*), *S*(*β*), *S*(*δ*) are standard deviations.

The moment magnitude values are required for the application of correlation Equation (1). The relationship between the local magnitude *ML* for Kamchatka earthquakes and the moment magnitude *MW* is [49]:

$$M\_L = M\_W - 0.4.\tag{12}$$

The following parameters of the elastic half-space were taken: the shear modulus, *<sup>μ</sup>* <sup>=</sup> 3.675 · <sup>10</sup><sup>10</sup> N/m2 , the second Lame parameter, *<sup>λ</sup>* <sup>=</sup> 3.675 · <sup>10</sup><sup>10</sup> N/m2 [40]. The simulation was carried out on a grid with the dimensions of 8◦ in latitude and 8◦ in longitude with the step of 0.01◦. Earthquake coordinates were the center of the grid.

#### **5. The Results of Computational Experiment**

#### *5.1. Earthquake No. 1*

Figure 4 shows the example of a simultaneous anomaly of acoustic emission and rock deformations recorded on 1 May 2007, 25 h before the earthquake that occurred on 2 May 2007, at 12:00 UT [13]. It is clear from Figure 4 that during the period from 1 to 9 o'clock, rather sharp compressions of rocks occurred, followed by the releases lasting from 1 to 5 min, which were accompanied by increases in the deformation rate and simultaneous increase in the emission level in kilohertz frequency range. The level of relative deformations during the compression reached the order of 10<sup>−</sup>7, and the deformation rate increased to 10−<sup>8</sup> s<sup>−</sup>1.

**Figure 4.** An example of a simultaneous acoustic-deformation anomaly before earthquake No. 1. (**a**) Variations of rock relative deformation *ε*, (**b**) variations of deformation rate *ε*˙, (**c**) variations of acoustic pressure *Ps*, accumulated over 4 s in the frequency range of 2.0–6.5 kHz.

The simulation results for the zones of relative shear deformations that occurred during earthquake No. 1 preparation are presented in Figure 5.

**Figure 5.** Zones of relative shear deformations on the Earth's surface *z* = 0 simulated for earthquake No. 1. The triangle on the map indicates the location of the Karymshina observation site.

It is clear from Figure 5 that the Karymshina observation site is on the boundary of the region of relative shear deformations of the order from 10−<sup>8</sup> to 10−7. That generally corresponds to the results of the natural experiment with a laser strainmeter-interferometer.

#### *5.2. Earthquake No. 2*

Figure 6 shows an example of a synchronous recording of acoustic emission and rock deformation from 6 October to 8 October 2009 before the earthquake that occurred on 8 October 2009, at 05:25 UT [13].

Figure 6a,b shows that, 35 h before the earthquake, there was a simultaneous anomaly of acoustic emission and rock deformation lasting for about 12 h. Figure 6c,d shows more detailed fragments of the record during the anomaly. For comparison, Figure 6e illustrates the subsequent calm period. The level of relative deformations during the anomaly was 10−<sup>7</sup> and sometimes reached the order of 10<sup>−</sup>6.

Figure 7 represents the simulation results for the zones of relative shear deformations occurring during this earthquake preparation.

It is clear from Figure 7 that the Karymshina site is in the region of relative shear deformations of the order from 10−<sup>8</sup> to 10−7, as in the case of earthquake No. 1, which corresponds to the results of the natural experiment using a laser strainmeter-interferometer.

In both cases presented, the calculated levels of relative shear deformation turned out to be slightly lower than the data of the natural experiments. In the case of earthquake No. 1, the Karymshina site is on the boundary of the deformation region from 10−<sup>8</sup> to 10<sup>−</sup>7, while according to deformation measurements, the relative deformation was about 10−7. In the case of earthquake No. 2, the Karymshina site is in the area of deformations from 10−<sup>8</sup> to 10−7, while according to deformation measurements, the relative deformation was about 10−<sup>7</sup> and sometimes reached the order of 10<sup>−</sup>6.

**Figure 6.** Record fragment of acoustic emission in different ranges and rock deformations from 00:00 on 6 October 2009, to 10:00 on 8 October 2009. (**a**) Variations of acoustic pressure *Ps*, accumulated over 4 s in the frequency range 2.0–6.5 kHz, (**b**) change in the strainmeter base Δ*L*. The red arrow shows the earthquake moment. At the bottom (**c**–**e**), enlarged fragments of the rock relative deformation *ε*, the deformation rate *ε*˙, and sound pressure *Ps*, accumulated over 4 s in the frequency range of 2.0–6.5 kHz, are presented.

**Figure 7.** Zones of relative shear deformations on the Earth's surface *z* = 0 simulated for earthquake No. 2. The triangle on the map indicates the location of the Karymshina observation site.

#### **6. Discussion**

The suggested pre-seismic deformation model, based on the classical theory of elasticity and simplification of the Earth's crust model to isotropic elastic half-space, has its advantages and disadvantages. Such a model has analytical solutions that simplifies calculations, obviates the need for the estimate of numerical solution stability. For example, a more complicated model of a medium in the form of isotropic uniform elastic sphere requires numerical solution of differential equation system [30]. In that case, the differences in estimates turn to be significant mainly for deep earthquakes at the distance of about 5◦, from their epicenters [30]. Both earthquakes considered in this paper are shallow. Their epicenters are at the distance of about 2◦, from the Karymshina observation site. That makes it possible use the Earth's model in the form of elastic half-space, not taking into account its surface curvature.

One more disadvantage of the proposed model is the application of deformation statistical equations, which do not take into account the deformation rate variation. In this respect, rock mass deformation rate affects significantly the activation of acoustic emission before earthquakes [15]. The Earth's crust deformation rate makes it possible to take into account different visco-elastic models of a medium, in particular, the Maxwell visco-elastic model [50]. However, it is very difficult to apply such models to estimate the Earth's crust deformations during earthquake preparation. For example, it is impossible to determine the exact time of earthquake preparation and the function of force impact change in a source. Thus, application of deformation static equations is justified.

When modeling, the authors did not take into account plastic deformations and heterogeneous structure of the Earth's crust. In fact, the Karymshina site is located in the zone of different-rank tectonic faults that may result in the recording of the deformations with the levels exceeding the calculated ones. This fact, also known as "problem of fardistance effect of earthquake sources", was considered by the researchers before. For example, it was proposed, when modeling, to take into account the Earth's crust regions with anomalous regime of the stress state (fault zones, layers with high fluid pressure) [51]. It was shown that, when introducing such regions of postcritical deformation with inelastic properties into a model, the decrease in disturbed deformation level at a large distance is 10<sup>4</sup> times less than in case of medium elastic model. Moreover, during the simulation, only shear deformations, which prevail during acoustic radiation generation, were considered, whereas the strainmeter records rock deformation within its base, not taking into account its type.

However, the suggested model makes it possible to estimate the stress–strain state of the Earth's crust during earthquake preparation at different distances from their sources. The computational experiment, using the proposed model, showed that the deformation during earthquake preparation at the Karymshina observation site exceeded the tidal ones by an order. Overall, the simulation results corresponded to the deformations measured by the laser strainmeter-interferometer. That is the ground to state that the observed joint acoustic and deformation anomalies are associated with the process of earthquake preparation in Kamchatka.

#### **7. Conclusions**

In seismically active regions of the Earth, to which the Kamchatka peninsula refers, pre-seismic anomalies are recorded in different geophysical fields. One of such fields is the acoustic emission of rocks, the anomalies of which are recorded 1–3 days before earthquakes at the distance of the first hundreds of kilometers from their epicenters. The results of joint acoustic-deformation measurements showed that growth of geoacoustic radiation intensity occurs during the increase in the level of deformations in rock masses by more than one order compared to the background values. According to the assumption that the Earth's crust in the first approximation can be considered as a homogeneous isotropic elastic half-space, and an earthquake source can be considered as a displacements along a rectangular fault plane, the authors proposed to simulate the deformation of the

Earth's crust around the source of impending earthquake. Another assumption of the authors was that increased deformations occur at a distance of hundreds of kilometers from the epicenters, and that causes acoustic emission and rock deformations anomalies. The total energy of elastic deformations accumulated during the earthquake preparation was estimated for the simulation. It determines the stress–strain state of the Earth's crust around the epicenter and is significantly greater than the released energy of seismic waves. Based on these assumptions, deformation regions were modeled for two earthquakes with local magnitudes *ML* ≈ 5 occurred on the Kamchatka Peninsula in 2007 and 2009. Simultaneous anomalies of acoustic emission and deformation of near-surface rocks were recorded before these earthquakes.

The simulation results were compared for the first time with the data of the strainmeterinterferometer and geoacoustic system installed at the Karymshina observation site in Kamchatka. It was shown that during the preparation of both earthquakes, the Karymshina observation site was within the region of shear deformations ≈10−<sup>7</sup> , which exceeded the tidal ones by an order. On the whole, simulation results corresponded to the results of the natural observations. Comparison of the results of the simulation and of the natural analysis is topical for understanding the causes of anomalous acoustic-deformation disturbance occurrences at the final stage of earthquake preparation. Construction of an adequate model for generation of such disturbances is relevant for development an early notification system on the threat of earthquakes.

**Author Contributions:** Conceptualization, Y.M. and M.G.; methodology, Y.M.; software, M.G.; validation, Y.M. and M.G.; formal analysis, M.G.; investigation, Y.M. and M.G.; writing—original draft preparation, M.G.; writing—review and editing, Y.M.; visualization, M.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research was carried out as a part of implementation of the Russian Federation state assignment AAAA-A21-121011290003-0.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data on nine small-focus earthquakes, which were used for estimated average dislocation parameters in the region (latitude interval: [52◦, 53◦], longitude interval: [160◦, 161◦]) are available at The Global Centroid-Moment-Tensor Catalog at https://www.globalcmt. org/cgi-bin/globalcmt-cgi-bin/CMT5/form?itype=ymd&yr=1976&mo=1&day=1&otype=ymd&oyr =2022&omo=6&oday=22&jyr=1976&jday=1&ojyr=1976&ojday=1&nday=1&lmw=0&umw=10&lms= 0&ums=10&lmb=0&umb=10&llat=52&ulat=53&llon=160&ulon=161&lhd=0&uhd=90&lts=-9999&u ts=9999&lpe1=0&upe1=90&lpe2=0&upe2=90&list=0 (accessed on 27 June 2022). The data on two earthquakes, for which the simulation was performed, are available at Unified Information System of Seismological Data of the Kamchatka Branch of Geophysical Service Russian Academy of Science at http://sdis.emsd.ru/info/earthquakes/catalogue.php (accessed on 5 June 2022).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

### *Article* **Spatial Correlations of Global Seismic Noise Properties**

**Alexey Lyubushin**

Institute of Physics of the Earth, Russian Academy of Sciences, Moscow 123242, Russia; lyubushin@yandex.ru

**Abstract:** A study of global seismic noise during 1997–2022 was carried out. A property of waveforms known as the Donoho–Johnston (DJ) index was used, which separates the values of the wavelet coefficients into "small" and "large". For each reference point in an auxiliary network of 50 points, a time series was calculated with a time step of one day for the median of the values at the five nearest stations. In a moving time window of 365 days, correlations between the index values at the reference points were calculated. A decrease in the average values of the DJ-index and an increase in correlations were interpreted as a sign of an increase in global seismic danger. After 2011, there was a sharp increase in the maximum distances between reference points with large correlations. The high amplitude of the response of the DJ-index to the length of the day for 2020–2022 could predict a strong earthquake in the second half of 2023. The purpose of this study was to improve the mathematical apparatus for assessing the current seismic hazard according to the properties of seismic noise.

**Keywords:** seismic noise; wavelet-based entropy; wavelet-based Donoho–Johnston index; correlations; day length

#### **1. Introduction**

Information about seismic noise makes it possible to study the processes preceding strong earthquakes [1–3]. Atmosphere and ocean (cyclone movement and the impact of ocean waves on the shelf) are the main sources for the energy of seismic noise [4–10]. At the same time, processes inside the earth's crust are reflected in changes in the properties of seismic noise and studying these properties helps in investigating the structural features of the earth's crust [11–13].

The Donoho–Johnston index (DJ-index) of a random signal can be defined as the ratio of the number of the coefficients of orthogonal wavelet decomposition, which in absolute value exceed the threshold introduced in [14], to the total number of wavelet coefficients. The threshold that separates "large" wavelet coefficients from others was defined first in [14] and is used to shrink noise from signals and images using wavelets. It has turned out that the DJ-index has the ability to most clearly highlight the effects of the spatial correlation of seismic noise, and it outperforms other noise statistics (such as entropy and the width of the multi-fractal spectrum of the singularity) in terms of highlighting predictive effects [15–18]. This article presents a detailed study of the DJ-index of global seismic noise, both the change in time and space of the correlation properties of noise and its response to the irregularity of the Earth's rotation, as well as the possible use of this response to assess the current seismic hazard.

#### **2. Initial Seismic-Noise Data**

The data used were the vertical components of continuous seismic-noise records with a 1 s sampling step, which were downloaded via the address http://www.iris.edu/ forms/webrequest/ (accessed on 1 January 2023) from 229 broadband seismic stations of three networks, http://www.iris.edu/mda/\_GSN (accessed on 1 January 2023), http: //www.iris.edu/mda/G (accessed on 1 January 2023), and http://www.iris.edu/mda/GE

**Citation:** Lyubushin, A. Spatial Correlations of Global Seismic Noise Properties. *Appl. Sci.* **2023**, *12*, 6958. https://doi.org/10.3390/ app13126958

Academic Editor: Roberto Zivieri

Received: 4 April 2023 Revised: 6 June 2023 Accepted: 7 June 2023 Published: 8 June 2023

**Copyright:** © 2023 by the author. 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/).

(accessed on 1 January 2023), which belong to the Incorporated Research Institutions for Seismology.

The duration of the time interval of observations was 26 years, from 1 January 1997 to 31 December 2022. The data were resampled to 1 min time series by computing the mean values within adjacent time intervals of the length 60 s.

An auxiliary network of 50 reference points was introduced. The positions of these points were determined by cluster analysis of the positions of 229 seismic stations using the hierarchical "far neighbor" method, providing compact clusters [19]. The positions of the 229 seismic stations and 50 control points are shown in Figure 1. The numbering of the reference points was carried out in order as the latitude of the point decreases.

**Figure 1.** Positions of 229 broadband seismic stations (blue circles) and a network of 50 reference points (numbered red circles).

#### **3. Wavelet-Based Measure of Time-Series Non-Stationarity**

Let us consider a random signal *x*(*t*) where *t* = 1, ... , *N* is an integer discrete-time index. For a finite sample of the signal, the entropy [20] could be defined as:

$$En = -\sum\_{k=1}^{N} p\_k \cdot \log(p\_k), \ p\_k = c\_k^2 / \sum\_{j=1}^{N} c\_j^2 \tag{1}$$

Here, *ck* are the coefficients of orthogonal wavelet decomposition. Within a finite set of Daubechies wavelet bases [20] with the number of vanishing moments from 1 to 10, the optimal basis is chosen from the minimum of (1).

The DJ-index was introduced in the problem of noise reduction and compression of information by setting all "small" wavelet coefficients to zero and performing an inverse wavelet transform [14,20]. This operation of noise reduction is performed for the optimal wavelet basis which was found from the minimum of entropy. The problem is how to define the threshold which separates "large" wavelet coefficients from "small" ones. In [14], this problem was solved using a formula of the asymptotic probability of the maximum deviations of Gaussian white noise. The DJ threshold is given by the expression *TD J* = *σ* √ 2 · ln *N*, where *σ* is the standard deviation of wavelet coefficients of the first detail level of wavelet decomposition. The first level is the most high frequency and is considered as the most noisy. Thus, in order to estimate the noise standard deviation, it is

possible to calculate the standard deviation of the first-level wavelet coefficients *c* (1) *<sup>k</sup>* . The robust estimate was proposed:

$$\sigma = med \left\{ \left| c\_k^{(1)} \right| \; ; \; k = 1, \dots, N/2 \right\} / 0.6745 \right\} \tag{2}$$

where *N*/2 is the number of wavelet coefficients on the first level, and *med* is the median. Formula (2) gives the relationship between the standard deviation and the median for the Gaussian random variable.

Finally, it is possible to define the dimensionless signal characteristic *γ*, 0 < *γ* < 1, as the ratio of the number of the most informative wavelet coefficients, for which the inequality |*ck*| > *TD J* is satisfied, to the total number *N* of all wavelet coefficients. For stationary Gaussian white noise, which is a kind of ideal stationary random signal, the index *γ* = 0. That is why the DJ-index could be called a measure of non-stationarity.

For each reference point, the index value *γ* is calculated daily as a median of the values in the five nearest operational stations. Thus, continuous time series are obtained with a time sampling of one day at 50 reference points. Figure 2 shows the graphs of the index *γ* for 9 reference points (point numbers are indicated next to each graph).

**Figure 2.** Graphs of daily index *γ* values for 9 reference points. The green lines represent the moving averages in a 57-day window.

Figure 3 shows a graph of average daily values of index *γ* calculated for all 50 reference points. The red lines represent the linear trends of the average values calculated for the time fragments 1997–2003, 2004–2015 and 2016–2022. The boundaries of time fragments were chosen so that the continuity of the broken line at the union of linear trends was approximately observed. For the time intervals 1997–2003 and 2016–2022, linear trends coincided with the average values and it can be seen that during the intermediate time interval 2004–2015 there had been a significant decrease in the global average *γ*. According to the results published in articles [15,17,18], a decrease in the average value of the index *γ*, which is a sign of the simplification of the statistical structure of noise (approximation of

its properties to the properties of white noise), was associated with an increase in seismic hazard. Thus, the graph of average values in Figure 3 can be interpreted as a transition to a higher level of global seismic hazard during 2004–2015.

**Figure 3.** Graph of the average daily values of the DJ-index for all 50 reference points. Red lines present linear trends for 3 time intervals: 1997–2003, 2004–2015 and 2016–2022.

#### **4. Probability Densities of Extreme Values of Seismic-Noise Properties**

Further analysis of the properties of seismic noise was based on the identification of regions that are characterized by the most frequent occurrences of extreme values for certain seismic-noise statistics. To do this, it was necessary to estimate the probability densities of the distribution over the space of extreme values by their values in the network of reference points. The following process was used. Denote by *ζk*, *k* = 1, ... , *m* = 50 coordinates of the reference points. Let *Zk* be the values of any property of the noise at the reference point with the number *k*. We will be interested in maps of the distribution of probability densities of extreme values (minimums or maxima) of values *Z* over the Earth's surface. To construct such maps, consider a regular grid of 50 nodes in latitude and 100 nodes in longitude covering the entire earth's surface. Next, maps will be built for various seismic-noise statistics *Z*, the values of which will be estimated in time windows of various lengths. For each such window, we define the time index *t*, which corresponds to its right end. For each time index *t*, find the reference point with the number *k*∗ *<sup>t</sup>* , at which the extremum of one or another property *Z* is realized, and let *ζ*∗ *<sup>t</sup>* be the coordinate of the reference point with the number *k*∗ *t* .

Let *ρ*(*ζk*,*rij*) be the distance between the points *ζ<sup>k</sup>* and the coordinates of the nodes of the regular grid *rij*. Then, the value of the probability density of extreme values for the quantity *Z* at the node *rij* of the regular grid can be calculated according to Gaussian distribution:

$$p(r\_{ij} \left| t \right) = \exp\left( -\rho^2(\zeta\_t^\*, r\_{ij}) / (2h^2) \right) / P(t) \tag{3}$$

Here, *h* is the distance, which corresponds to an arc of 15 angular degrees, which is approximately equal to 1700 km, and *P*(*t*) is a divisor that provides the normalization condition so that the numerical integral of function (3) over the Earth's surface is equal to 1. Formula (3) gives an "elementary map" of the probability density of the property *Z* extrema distribution for each current window with a time index *t*. This makes it possible to calculate the average probability density for a given time interval from *t*<sup>0</sup> to *t*1:

$$\overline{p}(r\_{i\bar{j}}|t\_0, t\_1) = \sum\_{t=t\_0}^{t\_1} p(r\_{i\bar{j}}|t) / n(t\_0, t\_1) \tag{4}$$

where *n*(*t*0, *t*1) is the number of time windows labeled from *t*<sup>0</sup> to *t*1.

Another way to estimate the variability of the probability density of extreme properties of seismic noise is to construct histograms of the distribution of the numbers of control points, in which the statistics extrema are realized. The construction of such histograms in a moving time window, the length of which should include a sufficiently large number of initial time windows within which the noise properties are estimated, makes it possible to determine the reference points at which the minima or maxima of the seismic-noise properties *Z* are realized most often.

This method has the disadvantage that it does not give a direct distribution of the places of the most probable realization of the noise-property extrema in the form of a geographical map, although each reference point has geographical coordinates. On the other hand, this method makes it possible to compactly visualize the temporal dynamics and, in particular, to determine the time intervals when the places of concentration of extreme values for the noise statistics change sharply. To solve a similar problem of visualizing temporal dynamics using averaged probability densities (4), one would have to build a laborious and long sequence of maps for time fragments consisting of the same number of time windows.

Therefore, in the future, the probability densities of extreme values of noise properties will be presented simultaneously in two forms: as averaged densities according to Formula (4) for two time intervals, before and after 2011, and as a temporal histogram of the distribution of numbers of reference points in which extremes are realized. In this case, the number of base intervals (bins) of histograms is taken equal to the number of reference points, that is, 50, which makes it possible to visualize the dynamics of the emergence and disappearance of bursts of the probability of extreme values at each reference point.

Figure 4a,b show maps of the distribution of the minimum values of the DJ-index, calculated by Formulas (3) and (4) for windows with time stamps before and after 2011. The choice of 2011 as a boundary time mark is connected, as will be seen from the following presentation, with the fact that the correlation properties of seismic noise changed dramatically after two mega-earthquakes: on 27 February 2010, *M* = 8.8 in Chile, and on 11 March 2011, *M* = 9.1 in Japan. Previously, this change has been already detected in [3,18]. Increased attention to the areas where the minimum values *γ* are most often realized, as was already mentioned above when discussing Figure 3, is due to the fact that an increase in seismic hazard is associated with a simplification of the statistical noise structure and a decrease in *γ*. For regions that are quite densely covered with seismic networks, for example, Japan, this property of statistics *γ* makes it possible to estimate the location of a possible strong earthquake [17].

However, the global network of seismic stations is sparse and therefore the selection of places where the probability of occurrence of small values is increased does not necessarily coincide with the known places of the occurrence of strong earthquakes. In Figure 4a,b, the maxima of the probability-density distribution of the minimum values *γ* are located in the vicinity of reference points with numbers 23, 25 and 37, 38 in the western part of the Pacific Ocean, where the epicenters of strong earthquakes and volcanic manifestations are really concentrated, including the largest volcanic eruption, Hunga Tonga-Hunga Ha'apai on 15 January 2022 [21]. Note that in Figure 4c, which shows a diagram of the change in the histogram of the numbers of reference points, in which the minimum value *γ* was realized, one can see the switching of the histogram maximum from the reference point number 23 to the reference point number 38 just after 2011. As for the concentration of probabilities of minimum values *γ* in the vicinity of reference points with number 46 in the southern Indian Ocean (in the vicinity of Kerguelen Island) and in the vicinity of point 40 in the Atlantic Ocean (Saint Helena) after 2011, these features can be associated with a mantle-plume activation regime [22].

**Figure 4.** (**a**,**b**) The probability densities of the distribution of the minimum values of the seismicnoise index *γ* in a network of 50 control points before and after 2011; (**c**) histogram of numbers of control points, in which the minimum values of the index *γ* were realized in a moving time window 365 days long.

#### **5. Estimation of Spatial Long-Range Correlations of Seismic-Noise Properties**

The simplification of the statistical structure of seismic noise manifests itself in a decrease in the average value of the index *γ* (Figure 3), and it is also reflected in an increase in entropy (1) and a decrease in the width of the multi-fractal singularity spectrum support width ("loss of multifractality"). These changes in noise properties are indicators of the transition of a seismically active region (including the whole world) to a critical state [3,17].

Another indication of an increase in seismic hazard is an increase in the average correlation as well as the radius of strong correlations between noise properties in different parts of the system. In order to find out how the spatial scale of correlations changes with time, for each reference point we calculated the correlation coefficient between the index *γ* values at this point and at all other reference points. We performed these calculations in a sliding time window 365 days long with a shift of 3 days. Next, for each reference point, we calculated the average value of the correlation coefficient with the values of *γ* at all other reference points. Since the result of calculating such an average value depends on the time window, we obtained graphs of changes in average correlations for all reference points. Figure 5 shows examples of graphs of changes in correlations for nine reference points.

From the graphs in Figure 5, a general trend towards an increase in the average correlations for each reference point is noticeable. As a result of averaging all the graphs of the type shown in Figure 5 for all 50 reference points, we obtained a general graph of average correlations, as shown in Figure 6.

It can be noticed from the graph in Figure 6 that the time interval 2004–2016, in which the index *γ* decreases (Figure 3), corresponds to an increase in average correlations between the properties of noise at various reference points and that there is a transition from small correlations before 2004 to fairly large correlations after 2016.

**Figure 5.** Graphs of the average correlation coefficients of the index *γ* for 9 reference points with values at other reference points in a moving time window of 365 days with a mutual shift of 3 days.

**Figure 6.** The average correlation coefficient between the values of the DJ-index *γ* at each control point with the values at other points. Red lines show linear trends for 3 intervals of time marks of the right ends of windows with a length of 365 days: 1998–2004, 2004–2016 and 2016–2023. On the last time fragment, the linear trend coincides with the average value.

In order to assess the scale of strong correlations between seismic-noise properties, for each reference point we defined another point that had a maximum correlation with this point, provided that this maximum exceeded the threshold of 0.8, and we calculated the distance between these points. Figure 7 shows the graphs of the changes in the average and maximum values of the distances between points for which the correlation maxima exceeded 0.8.

**Figure 7.** (**a**) The average distance between the reference points for which the maxima of correlations between the values of the DJ-index of seismic noise exceeded the threshold of 0.8; (**b**) the maximum distance between the reference points with the maximum correlation coefficient above the threshold of 0.8.

It can be seen from the graph in Figure 7b that after 2011 there was a sharp increase in the maximum distance at which strong maximum correlations occur. Thus, the growth of average correlations, which is visible in Figure 6, was supplemented by an explosive growth of the maximum scale of strong correlations.

It is of interest to identify those reference points in the vicinity of which strong correlations most often occur. As before, for this purpose we used estimates (3–4) of the probability density of the distribution over space of the maximum correlation values for the time intervals before and after 2011, as well as the distribution histograms of the numbers of control points in which the maximum correlation values occurred.

To select the length of the time window for calculating histograms, we recall that the correlation coefficients between the values at each reference point with the values of the DJindex at other reference points are calculated in windows of 365 days with a shift of 3 days. Let us choose the length of the window for calculating histograms so that the dimensional length of the window is equal to five years. It is easy to calculate that this length will be 488 values of "short" time windows with a length of 365 days with an offset of 3 days. In this case, the histogram is calculated in a window of length 365 + 3 · (488 − 1) = 1826 days, which is approximately equal to five years, given that every fourth year is a leap year. As for the number of base intervals (bins) of histograms, in this case there will not be 50, as in Figure 4c, but 44, since for reference points with numbers 1–6 there are no realizations of correlation maxima.

Estimates of the position and temporal dynamics of places with an increased probability of maximum correlations are presented in Figure 8. Comparison of Figure 8a,b shows that the places of the concentration of maximum correlations have changed significantly since 2011. In particular, there has been a rapid shift of such a center from the Southwest Pacific to Central Eurasia, which can be seen in the histogram diagram in Figure 8c.

**Figure 8.** (**a**,**b**) The probability densities of the correlation maxima for index *γ* values at each reference point with index values at other reference points before and after 2011; (**c**) a histogram of the numbers of reference points, in which the correlation maxima for values of *γ* at each reference point with the values at other points in a moving time window of 5 years were realized.

#### **6. Seismic-Noise Response to the Irregularity of the Earth's Rotation**

The connection between the irregular rotation of the Earth and seismicity was investigated in [23]. The influence of strong earthquakes on the Earth's rotation was studied in [24,25]. Figure 9 shows a graph of the time series of the length of the day for the time interval 1997–2022.

**Figure 9.** Time series plot of the length of a day for the time interval 1997–2022. Day-length data (LOD) are available from the website at: https://hpiers.obspm.fr/iers/eop/eopc04/eopc04.1962-now (accessed on 1 January 2023).

Further, the LOD time series is used as a kind of "probing signal" for the properties of the seismic process [3,15–18]. To study the relationship between the properties of the seismic process and the irregularity of the Earth's rotation, the squared coherences maxima between the LOD and the properties of seismic noise at each reference point were estimated. The coherences were estimated in moving time windows of 365 days with a shift of 3 days

using a vector (two-dimensional) fifth-order autoregressive model [15–18]. In what follows, these maximum coherences will be referred to as the seismic-noise responses to the LOD. Figure 8 shows examples of the LOD response graphs at nine reference points.

We define the integrated response of global seismic noise to the LOD as the average of the responses of the type shown in Figure 10 from all 50 reference points. Let us compare the behavior of the average response to the LOD with the amount of released seismic energy. To do this, we calculate the logarithm of the released energy as a result of all earthquakes in a moving time window 365 days long with a shift of 3 days. In Figure 11a,b, there are graphs of the logarithm of the released energy and the average response to the LOD. It is visually noticeable that bursts of the average maximum coherence in Figure 11b precede the energy spikes in Figure 11a. In order to quantitatively describe this delay, we calculated the correlation function between the curves in Figure 11a,b. The graph of this correlation function for time shifts of ±1200 days is shown in Figure 11c.

**Figure 10.** Plots of maximum squared coherence between index *γ* and LOD in a moving time window of 365 days with mutual shift of 3 days for 9 reference points.

The estimate of the correlation function in Figure 11c between the logarithm of the released seismic energy in Figure 11a and the average response of seismic noise to the LOD had a significant asymmetry and was shifted to the region of negative time shifts, which corresponded to the coherence maxima being ahead of the maxima of seismic energy emissions. The correlation maximum corresponded to a negative time shift of 534 days. This estimate of cross-correlations produced a foundation to propose that the burst shown in Figure 11b by the magenta line may precede a major earthquake with an average delay of 1.5 years.

The curve in Figure 11b can be split into three sections. The first two intervals with time marks of the right ends of windows before and after 2012 differ significantly from each other by their mean values presented by red lines.

**Figure 11.** (**a**) Logarithm of the released seismic energy (in joules) in a sequence of time intervals 365-days long, taken with an offset of 3 days; (**b**) average values of maximum coherences between LOD and daily DJ-index values at 50 reference points; (**c**) the correlation function between the logarithm of the released seismic energy and the average value of the maximum coherence between the LOD and the seismic-noise DJ-index at 50 reference points.

The third time interval in Figure 11b is presented by the magenta line and refers to the processing of the most recent data which correspond to right-hand end windows after 2021, i.e., referring to the time span of 2020–2022. A short third segment was identified based on the results of data processing for 2021–2022 due to an unusually large spike in the response of noise properties to the LOD. The maximum value of the response in Figure 11b was reached at the time corresponding to the date 9 May 2022. The "naive" forecast, in which 534 days are added to this date, gives the most probable date for the future strong event as being 24 October 2023. A more realistic prediction would be to blur this date with some uncertainty interval. However, it is not yet possible to give such an interval due to the short history of using the seismic-noise response to the LOD as a predictor.

Similar to how it was done to identify the places of concentration of the maxima of spatial correlations of seismic noise in Figure 8, we can search for places where the maximum response of seismic-noise properties to the uneven rotation of the Earth is most likely. The results of such a search are shown in Figure 12. Comparison of Figure 12a,b shows that the maximum response to the LOD is concentrated in the subarctic regions.

In addition to the issue of the synchronization of seismic-noise properties, which are reflected in the graphs in Figures 5–7, we also studied the synchronization of seismicnoise responses to the irregularity of the Earth's rotation. Previously, this issue has been considered for the synchronization of seismic-noise responses to LODs in Japan and California [16].

For each reference point, we estimated the correlation of the response at this point with the responses at all other points. It should be taken into account that the correlated responses themselves were calculated in sliding time windows of 365 days with a shift of 3 days. These were "short" time windows. In this case, we needed to take a certain number of consecutive "short" windows and form a "long window" from them, such that it contained a sufficiently large number of estimates of the noise responses to the LOD in the "short" time windows. Next, 250 adjacent "short" windows were selected to form a "long" window. Correlation coefficients between the responses to the LOD at various

reference points were calculated in a sliding "long" window with a shift of 1, that is, in a time window with a length of 365 + 3·(250 − 1) = 1112 days, which is approximately 3 years. A shift of 1 count meant a real shift of 3 days.

**Figure 12.** (**a**,**b**) Probability densities of the seismic-noise index *γ* response maxima to LOD before and after 2011; (**c**) histogram of numbers of reference points, in which the maximum values of the index *γ* response to LOD were realized in a sliding time window of 5 years.

The goal of such calculations for each reference point was the average correlations of responses over all other reference points. As a result of such averaging, 50 dependences of the same type as shown in Figure 5 were obtained. In order to obtain an integrated measure of the correlation of global noise responses to the LOD, the average dependences for all reference points needed to be averaged secondarily. As a result, we obtained the graph shown in Figure 13.

**Figure 13.** The average value of the correlation coefficients between the seismic-noise index *γ* responses per LOD for each of the 50 control points with responses at other points.

Figure 13 shows the change in time of the measure of synchronization of the response of seismic noise to the irregularity of the Earth's rotation. This graph also shows the growth of correlations, similar to the growth of the initial correlations in Figure 6.

#### **7. Conclusions**

The results of a detailed analysis of the properties of seismic noise based on the use of the Donoho–Johnston index, defined as the proportion of the maximum values of the moduli of the wavelet coefficients separating the "noise" coefficients, was presented. The main results of the conducted research were as follows:


It should be noted that the noise level in real observations of the seismic background is inevitably extremely high. Transition to the low-frequency part of the spectrum allows the elimination of most of the man-made noise. Moreover, this transition makes it possible to study that part of the spectrum that is in the intermediate, little-studied region of periods between traditional seismology and gravimetry. The author's previous studies give hope that it is in this region of periods that the effects preceding strong earthquakes can be detected. This work is presented as a further development of these methods. The complexity of the data and processes underlying the formation of low-frequency seismic noise in the Earth's crust does not allow direct use of conventional statistical procedures for estimating confidence intervals, standard deviations and other measures of uncertainty. It is possible that further study of the statistical structure of low-frequency seismic noise will enable us to build such models that will make it possible to estimate the measures of uncertainty of the obtained statistical inferences, for example, using the bootstrap technique.

**Funding:** This research received no external funding.

**Data Availability Statement:** Open-access data were used from the sources: http://www.iris.edu/ mda/\_GSN, http://www.iris.edu/mda/G, http://www.iris.edu/mda/GE, https://hpiers.obspm. fr/iers/eop/eopc04/eopc04.1962-now and https://earthquake.usgs.gov/earthquakes/search/. All data accessed on 1 January 2023.

**Acknowledgments:** This work was carried out within the framework of the state assignmen Institute of Physics of the Earth of the Russian Academy of Sciences (topic FMWU-2022-0018).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

**Rodolfo Console 1,2,†, Paola Vannoli 1,\*,† and Roberto Carluccio 1,†**


**Abstract:** We studied the long-term features of earthquakes caused by a fault system in the northern Adriatic sea that experienced a series of quakes beginning with two main shocks of magnitude 5.5 and 5.2 on 9 November 2022 at 06:07 and 06:08 UTC, respectively. This offshore fault system, identified through seismic reflection profiles, has a low slip rate of 0.2–0.5 mm/yr. As the historical record spanning a millennium does not extend beyond the inter-event time for the largest expected earthquakes (*M* 6.5), we used an earthquake simulator to generate a 100,000-year catalogue with 121 events of *Mw* ≥ 5.5. The simulation results showed a recurrence time (*Tr*) increasing from 800 yrs to 1700 yrs as the magnitude threshold increased from 5.5 to 6.5. However, the standard deviation *σ* of inter-event times remained at a stable value of 700 yrs regardless of the magnitude threshold. This means that the coefficient of variation (*Cv* = *σ*/*Tr*) decreased from 0.9 to 0.4 as the threshold magnitude increased from 5.5 to 6.5, making earthquakes more predictable over time for larger magnitudes. Our study supports the use of a renewal model for seismic hazard assessment in regions of moderate seismicity, especially when historical catalogues are not available.

**Keywords:** numerical modelling; earthquake simulator; statistical methods; earthquake clustering; northern Adriatic sea

#### **1. Introduction**

In regions of moderate seismicity, the recurrence times of damaging earthquakes are commonly longer than the time windows covered by reliable historical observations, even in Italy, which owns one of the longest records of past strong seismic events. These limitations of historical seismic data preclude the application of statistical analysis for seismic hazard assessment, even if they can be partly overcome by geological and paleoseismological investigations on specific faults. The information on past earthquakes is particularly poor for offshore seismic activity (where paleoseismological and geomorphological studies are not feasible), as is the case for the thrust fault systems located seaside of the northern Marche coast (central Italy), recently affected by a seismic sequence that has significantly concerned the population in a large area of central Italy (Istituto Nazionale di Geofisica e Vulcanologia website at http://terremoti.ingv.it; main event of the sequence at http://terremoti.ingv.it/en/event/33301831, accessed on 14 February 2023; ISIDe\_Working\_Group [1], Tertulliani et al. [2].

The use of earthquake simulators has become increasingly popular in generating synthetic earthquake catalogues with a vast number, even reaching millions, of events. This allows for statistical analyses to be conducted on simulated catalogues that are considerably more reliable than those conducted on real catalogues. However, criticism has been expressed regarding the practicality of simulated catalogues. Certain seismologists have pointed out that the algorithms used in earthquake simulators are based on oversimplified physical models and contain arbitrary assumptions, which pose significant challenges for

**Citation:** Console, R.; Vannoli, P.; Carluccio, R. The 2022 Seismic Sequence in the Northern Adriatic Sea and Its Long-Term Simulation. *Appl. Sci.* **2023**, *13*, 3746. https:// doi.org/10.3390/app13063746

Academic Editor: José A. Peláez

Received: 17 February 2023 Revised: 7 March 2023 Accepted: 14 March 2023 Published: 15 March 2023

**Copyright:** © 2023 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/).

creating a dependable representation of actual seismicity. Nevertheless, a wide consensus does exist on the potential of earthquake simulators to provide support for historical observations in the context of seismic hazard assessment, as long as a realistic physical model is available and the results are based on reliable source parameters. Some examples of research in this area include Schultz et al. [3], Christophersen et al. [4], Field [5]. In the frame of the ongoing debate on the usefulness of earthquake simulators, some valuable insights have also been gathered through the application of earthquake simulators to the seismicity of Italy, Greece, California, and Japan (Console et al. [6–9], Parsons et al. [10]). These insights have been gained by comparing the simulated catalogues with real catalogues from the respective regions, and the results have shown that the same algorithms can be applicable in different geographic regions and on different magnitude scales.

In a recent study by Console et al. [6], a new earthquake simulation code was developed, which included an enhanced Coulomb stress transfer among rupturing fault elements and increased production of multiple main shocks, frequently observed in earthquake sources of different mechanisms in Italy. To better understand the clustering behaviour of these events, a specific algorithm was created to detect and list inter-correlated earthquake clusters based on certain empirical rules.

The simulated catalogue produced by this code exhibits spatiotemporal features that imitate those frequently observed in reality. Using a stacking procedure, Console et al. [6] were able to represent the number of moderate-magnitude events that preceded and followed stronger earthquakes in a 100,000-year simulated catalogue. Their results for the central–northern Apennines region showed interesting long-term acceleration and quiescence patterns of seismic activity before and after mainshocks. They also analysed short-term patterns in periods of some weeks preceding and following every strong earthquake, confirming the simulator code's ability to reproduce typical foreshock–aftershock sequences.

In addition, the simulated catalogue for the central–northern Apennines region showed a decrease of b-values lasting some weeks before the occurrence of strong earthquakes, followed by a sudden increase at the time of these earthquakes. This pattern was previously observed by Gulia and Wiemer [11] in actual earthquake sequences. The fault model used in the simulations by Console et al. [6] consists of 43 fault systems divided into 198 quadrilaterals, which provide a suitable approximation of the seismic structures reported in the Database of Individual Seismogenic Sources (DISS) composite sources (DISS\_Working\_Group [12]; http://diss.rm.ingv.it/diss/, accessed on 14 February 2023).

Among these fault systems, the one labelled as ITCS106 (a NW–SE elongated fault system 70 km long in the Adriatic sea at about 25 km from the Italian coast) is believed to be responsible for the recent seismic sequence, which started with two relatively similar large mainshocks on 9 November 2022 at 06:07:25 UTC (*Mw*5.5) and 06:08:28 (*Ml*5.2), respectively. According to the definition adopted by Console et al. [6] this sequence, with two main earthquakes of similar magnitude, whose origin times differ by only one minute, can be retained a multiple shock.

No earthquakes from A.D. 1000 to 2020, with *Mw* ≥ 5.5, are reported in the Parametric Catalogue of Italian Earthquakes (CPTI15; https://emidius.mi.ingv.it/CPTI15-DBMI15/, accessed on 14 February 2023) within a 5 km buffer from the ITCS106 fault system. The lack of historical and instrumental events pertinent to the ITCS106 earthquake source is consistent with the low slip rate of this source, as will be shown in more detail in the following sections. In this paper, we present a case study of application of an earthquake simulator as a contribution to seismic hazard assessment in an area of moderate seismicity, for which historical catalogues are not sufficient for this purpose.

#### **2. Seismotectonic Settings**

The northern Marche coastal belt and the Adriatic offshore area are characterised by a series of NW–SE trending, NE verging folds forming the easternmost edge of the Apennines thrust front (Figure 1; e.g., DISS\_Working\_Group [12], Vannoli et al. [13], Fantoni and Franciosi [14], Kastelic et al. [15]). This Apennines thrust front has progressively migrated toward the east-northeast all through the Tertiary-Quaternary (e.g., Elter et al. [16]) and several geological evidence suggest that the folds are still growing and hence that the blind thrust fronts are active (e.g., Vannoli et al. [13]). Especially in the offshore area the interpretation of seismic reflection profiles is strategic to understanding the geometry of the active faults (e.g., Maesano et al. [17]).

**Figure 1.** Seismotectonic framework of the northern Marche coastal and offshore area. The seismic sequence from 9 November to 14 February 2023 is shown in light blue (http://terremoti.ingv.it/, accessed on 14 February 2023); historical and instrumental earthquakes from CPTI15 are shown with coloured squares (Rovida et al. [18]), earthquakes with *Mw* ≥ 5.5 are labelled (see Table 1); surface projection of Composite Seismogenic Sources from DISS 3.3.0 are shown with orange ribbons (https://diss.ingv.it, accessed on 14 February 2023; DISS\_Working\_Group [12]). Focal mechanisms of the 9 November 2022 earthquake and the 30 October 1930 event are from TDMT (http://webservices.ingv.it/webservices/ingv\_ws\_map/data/tdmt/19591/111604601\_15 4\_tdmt\_reviewer\_solution.pdf, accessed on 14 February 2023), and Vannoli et al. [19]), respectively.

The occurrence of historical and instrumental earthquakes (e.g., 1672, 1690, 1786, 1875, 1916, 1924, 1930, all with *Mw* ≥ 5.5; Table 1; Figure 1) suggests that the thrust faults are also seismogenic. The earthquakes in the area have often generated a tsunami; subsequently, their causative faults have been able to produce significant vertical displacement of the ground surface or sea floor, and are very close to the coast or offshore (Table 1; Vannoli et al. [19]).

A relevant aspect of seismic activity in the study area, as in the rest of Italy (e.g., Vannoli et al. [20]), is the occurrence of seismic sequences with more than one main shock of similar magnitude (Table 1). These sequences can be distinguished from typical aftershock sequences because the difference between the largest and second-largest magnitudes is relatively small and does not follow the so-called "Bath Law" (see Console et al. [7]).

A seismic sequence with two events of similar magnitude also began on 9 November 2022. As a matter of fact, at 06:07:25 UTC a *Mw*5.5 earthquake struck the northern Marche offshore area, and only one minute later, a *Ml*5.2 earthquake struck about 7 km to the south (http://terremoti.ingv.it, accessed on 14 February 2023; Figure 1; Table 1). According to the hypocentral location of these two mainshocks, the ITCS106 fault system is retained the causative source of the sequence. More than 900 aftershocks occurred offshore Pesaro and Senigallia from 9 November 2022 to 14 February 2023 (Figure 1).

The occurrence of this sequence demonstrates that the buried offshore thrust front, included in the Database of Individual Seismogenic Sources since 2013 despite the absence of historical and recent seismicity (DISS\_Working\_Group [12]), is active and seismogenic. Specifically, the ITCS106 Pesaro mare-Cornelia source straddles the Adriatic sea just north of the city of Ancona, and it is part of the Umbro-Marche Apennines outer offshore thrust. ITCS106 is a blind fault system believed capable of infrequent moderate-size earthquakes. Its slip rate, a main ingredient of the earthquake simulators and seismic hazard model, is in the range of 0.20–0.52 mm/yr (DISS\_Working\_Group [12] and reference therein). Therefore, ITCS106 can produce only subtle displacement of the Adriatic sea floor, even when it is accumulated over several seismic cycles.

**Table 1.** Historical and instrumental earthquakes with *Mw* ≥ 5.5 of the study area (Figure 1; CPTI15; Rovida et al. [18]; CFTI5Med; Guidoboni et al. [21]), and tsunami with their intensity (EMTC; Maramai et al. [22]). T: true; SA: Intensity Sieberg–Ambraseys scale; PI: Intensity Papadopoulos and Imamura scale; \* unspecified day; \*\* the earthquake triggered a large landslide that fell into the sea on the eastern side of Monte Conero (CFTI5Med [21]); \*\*\* some historical sources describe the occurrence of a second large earthquake one hour after the first event (CFTI5Med [21]); \*\*\*\* EMS-98 scale; Tertulliani et al. [2].


#### **3. Earthquake Simulation Method**

Our simulator code's algorithm was originally introduced by Console et al. [8], and subsequently modified in different versions with the inclusion of new features, as described by Console et al. [7,9]. The present study employs the following fundamental principles of the simulator algorithm:


The reason for choosing a trapezoidal shape for the seismic sources is to achieve a more precise representation of curved seismic structures. For instance, this helps to reduce gaps that may exist between adjacent rectangular segments when they do not have the same strike. It should be noted that using fault segments of rectangular or trapezoidal shapes to model seismic sources constitutes a simplification of the algorithm adopted in our physics-based simulation code. It is not a constraint that implies a rupture stops at the edges of every segment.

The simulator algorithm has two free parameters that control how a rupture nucleates, propagates, and stops. The first of these parameters is "Strength Reduction" (S-R), which reduces the strength at the edges of growing ruptures through a sort of weakening mechanism. Increasing this parameter favors the growth of ruptures and, consequently, it has the effect of decreasing the b-value of the frequency-magnitude distribution. The second parameter is "Aspect Ratio" (A-R), which is introduced to limit the effect of the S-R parameter when a rupture grows beyond certain specified limits. This parameter assumes a significant role in the frequency-magnitude distribution in the range of large magnitudes, as shown in Figure 5 of Console et al. [8].

It is important to note that neither of the free parameters mentioned above influences the long-term annual seismic moment rate of the simulated seismicity, which depends only on the slip rate assigned to the simulator within the kinematic fault model. In this study, we analysed the simulations computed by Console et al. [6] with the free parameters S-R = 0.2 and A-R = 5, obtaining a synthetic catalogue characterised by a b-value of the Gutenberg–Richter distribution equal to 0.95 and a realistic production of multiple events.

#### **4. Simulation of the Seismicity in the Adriatic Thrust Fault Systems**

In the application of the simulator, the quadrilateral segments representing the ITCS106 fault system were discretized in cells of 1.0 km × 1.0 km. As to the slip rate, we adopted the value of 0.5 mm/year corresponding to the top value of the range reported by the DISS 3.3.0 database for this source. We found in previous applications of the simulator to central and northern Apennines Console et al. [6,7,9] that the largest slip rate is the one that better matches the annual rate of seismic moment obtained from instrumental observations.

We chose for the synthetic catalogues a minimum magnitude of 4.2, which is produced approximately by the rupture of two cells. Magnitude 4.3 is not present in the output catalogue because the rupture of three cells produces an event of magnitude 4.4. The duration of the synthetic catalogue was 100,000 years. A warm up period, the events of which are not included in the output catalogue, is necessary to bring the system in a stable situation. This period should be longer than a few inter-event times of the strongest magnitude. In this regard, we chose a warm up period of 10,000 years.

In this study, still taking into account the interactions from other sources, we limited our analysis to the above-mentioned ITCS106 source in the Adriatic sea, where the November–February 2023 earthquake sequence took place. At first, for this source we considered the magnitude–frequency distribution. Figure 2 shows a clear deviation of the magnitude–frequency distribution obtained by our simulation algorithm from a plain linear Gutenberg–Richter distribution. In fact, the incremental magnitude–frequency distribution shows a fairly linear trend with a b-value slightly greater than 1.0 up to magnitude 5.5, but it exhibits a sharp "bump" around magnitude 6.5, This value could be defined as "characteristic" for the considered fault system, and is related to its dimensions. It can also be noted that the cumulative magnitude distribution of the synthetic catalogue is characterised by a very small b-value (i.e., a b-value slightly larger than 0) in the magnitude range 5.5<M< 6.5.

**Figure 2.** Incremental magnitude–frequency distribution of the earthquakes in the synthetic catalogue obtained from the simulation algorithm (triangles) and the corresponding cumulative distribution (diamonds).

Our results are consistent with those of other studies where the magnitude–frequency distribution is analysed on individual faults, ignoring the contribution of smaller surrounding sources. For instance, Parsons et al. [10], by means of various statistical methods, demonstrated strongly characteristic magnitude–frequency distributions on the San Andreas fault, with a nearly flat cumulative distribution in the 6.5–7.5 magnitude range, and higher rates of large earthquakes than what would be expected from a Gutenberg–Richter distribution. The difference of approximately 1 magnitude unit between the results obtained for the ITCS106 fault system (about 70 km long) and the Saint Andreas fault (several hundreds of km long) is easily justified by their different dimensions.

We focused our attention on the inter-event times between events for the simulated catalogue of 100,000 yrs for earthquakes of magnitude equal to or larger than magnitudes 5.5, 6.0 and 6.5, respectively. The smallest size of 5.5 was the magnitude threshold adopted for the pivot earthquake of a multiple event by Console et al. [6]. According to popular scaling relationships (Leonard [23]) magnitude 6.5 is expected for the largest size of an earthquake that may rupture the entire area of the ITCS106 fault system. Table 2 reports the results of the statistical analysis carried out on the simulated catalogue, and Figure 3 shows the number of inter-event times in bins of 250 yrs.

**Table 2.** Results of the statistical analysis of the simulated 100,000 yrs earthquake catalogue for the DISS composite source labelled ITCS106 for three different magnitudes thresholds. *Tr* = recurrence time, *σ* = standard deviation, *Cv* = coefficient of variation.

**Figure 3.** Density distributions of the number of inter-event times in bins of 250 yrs. The first bin corresponds to the range of 0–250 yrs and so on.

The plot of Figure 3 shows a bi-modal distribution of inter-event times with a local minimum at 1000–1250 yrs. Smaller inter-event times are pertinent to clustering (with particular regard to small magnitude earthquakes), while the maximum at 1500–1750 yrs denotes quasi-periodicity (typical of larger magnitude earthquakes).

As expected, the number of events decreases and the average inter-event time *Tr* increases when the magnitude threshold is increased (Table 2). However, the standard deviation *σ* maintains a value close to 700 yrs. In this way, the coefficient of variation (*Cv*), which is defined by the ratio between *σ* and *Tr*, decreases from 0.90 for *M* ≥ 5.5 to 0.40 for *M* ≥ 6.5. It is well known that *Cv* is equal to 1.0 for a sequence of events following a memoryless Poisson distribution and is equal to 0 for a perfectly periodical sequence. Thus, the earthquake sequence of our simulated catalogue behaves more and more as a quasi-periodical process when, increasing the threshold magnitude of the analysis, the ruptured areas of the considered earthquakes approach the size of whole seismic structure.

Shortest inter-event times (smaller than 1000 yrs) are typical of small magnitude earthquakes (Table 2, and Figure 3), while the opposite situation happens for the inter-event times larger than 2000 years (Figure 3). An exceptionally long inter-event time of 4625 yrs is estimated between two earthquakes of magnitude larger than 6.5 (Figure 3).

The simulations analysed in this study were obtained adopting the largest slip rates reported for the earthquake source in the DISS database. If we adopted smaller slip rates, the time scale of our inter-event times would become inversely longer.

The synthetic 100,000-year catalogue highlights a spatiotemporal behaviour that can hardly be compared with analog features observed in reality, due to their brevity, but can be considered realistic in light of an earthquake cycle hypothesis.

Long-term patterns preceding and following a strong earthquake are shown in Figure 4. This figure is obtained by stacking the number of *M* ≥ 4.2 earthquakes before and after an *M* ≥ 5.2 earthquake within a distance of 20 km between their respective epicenters in bins of 5 years. Here, we can observe an increase in seismic activity starting a couple of centuries before a strong earthquake, and a sharp increase in occurrence rate in the first 5 years after that event. The latter feature is clearly related to aftershock sequences following strong events. After this aftershock phase, the plot shows a long-lasting quiescence with some modest fluctuations. This long-term pattern highlights the existence of a preparatory phase before every strong earthquake. Presumably, during this preparatory phase, stress accumulation causes an increase in the rate of moderate-magnitude events. This hypothesis is supported by a detailed study of the time history of the average stress and its standard deviation during several seismic cycles in the Corinth (Greece) and Nankai (Japan) fault systems (Console and Carluccio [24], Console et al. [25]).

**Figure 4.** Stacked number of *M* ≥ 4.2 earthquakes that preceded and followed an *M* ≥ 5.2 earthquake within an epicentral distance of 20 km in the 100,000 years simulated catalogue.

#### **5. Conclusions**

A seismic sequence initiated by a pair of mainshocks of magnitude 5.5 and 5.2 took place in the Northern Adriatic sea on 9 November 2022 at 06:07 UTC. The two mainshocks were separated by only one minute in time and about 10 km in space (Figure 1). The sequence is believed to be generated by the ITCS106 seismogenic source of the DISS database, belonging to the series of NW–SE trending, NE verging folds forming the easternmost edge of the Apennines thrust front (Figure 1).

Due to its low slip rate, estimated between 0.2 and 0.5 mm/yr in the DISS database, the ITCS106 source is characterised by moderate seismic activity, so as no earthquakes with *Mw* ≥ 5.5 were reported in the Parametric Catalogue of Italian Earthquakes from 1000 to 2020, before the two mainshocks of 9 November 2022. The occurrence of these earthquakes supports the hypothesis that even sources of low slip rate identified and characterised in the DISS database, for which little or no historical information is reported in the historical catalogues, may have seismic potential to be considered relevant in the context of seismic hazard assessment.

To quantify this hypothesis, we performed a statistical analysis of the results obtained by a previously developed simulator algorithm (Console et al. [6]) for the ITCS106 fault system. The results of our simulation consist of a catalogue of 100,000 yrs duration, containing 121, 84 and 57 earthquakes of magnitude equal to or exceeding a magnitude of 5.5, 6.0 and 6.5, respectively. The statistical analysis of the synthetic catalogue, shortly reported in Table 2, puts in evidence that the earthquakes belonging to the lowest magnitude class have an average inter-event time of 814 ± 733 yrs, while those belonging to the highest magnitude class have an average inter-event time of 1713 ± 690 yrs. Consequently, the simulated catalogue exhibits a nearly memoryless Poissonian behaviour for moderatemagnitude earthquakes (whose coefficient of variation is 0.9), and a time-dependent quasiperiodic behaviour for large-magnitude earthquakes (whose coefficient of variation is 0.4). The statistical features obtained by the application of our earthquake simulator appear consistent with the hypothesis of the earthquake cycle. In this respect, our work may be relevant for possible applications related to the recent method of earthquake nowcasting (Rundle et al. [26,27], Varotsos et al. [28], Christopoulos et al. [29]), in which interoccurrence natural time between strong earthquakes is used for the estimation of the current stage of the earthquake cycle (see e.g., Varotsos et al. [30,31]).

Another significant aspect of our results is the acceleration of seismic activity before strong earthquakes lasting a couple of centuries in the ITCS106 seismic source. We conclude by suggesting that statistical parameters derived by physics-based earthquake simulators, even recognizing the limitations connected with the extreme simplicity of simulators, could be of support for the comprehension and modelling of seismic processes in situations where there is a lack of suitable observations for earthquake hazard assessment.

**Author Contributions:** Conceptualization, R.C. (Rodolfo Console), P.V. and R.C. (Roberto Carluccio); Validation, R.C. (Rodolfo Console) and P.V.; Data curation, R.C. (Roberto Carluccio); Writing—review & editing, R.C. (Rodolfo Console), Paola Vannoli and R.C. (Roberto Carluccio). All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The Italian Parametric Earthquake Catalogue (CPTI15) can be found at https://emidius.mi.ingv.it/CPTI15-DBMI15/ (accessed on 14 February 2023), the Database of Individual Seismogenic Sources (DISS) can be found at http://diss.ingv.it (accessed on 14 February 2023).

**Acknowledgments:** The authors are grateful to Jörg Renner, Eleftheria Papadimitriou and anonymous reviewers for their careful reading of the manuscript and useful suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

MDPI St. Alban-Anlage 66 4052 Basel Switzerland www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

Disclaimer/Publisher's Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Academic Open Access Publishing

mdpi.com ISBN 978-3-0365-9481-1