Seismicity Patterns Prior to the Thessaly (M_w6.3) Strong Earthquake on 3 March 2021 in Terms of Multiresolution Wavelets and Natural Time Analysis

Abstract

On 3 March 2021, a strong, shallow earthquake of moment magnitude, M_w6.3, occurred in northern Thessaly (Central Greece). To investigate possible complex correlations in the evolution of seismicity in the broader area of Central Greece before the M_w6.3 event, we apply the methods of multiresolution wavelet analysis (MRWA) and natural time (NT) analysis. The description of seismicity evolution by critical parameters defined by NT analysis, integrated with the results of MRWA as the initiation point for the NT analysis, forms a new framework that may possibly lead to new universal principles that describe the generation processes of strong earthquakes. In the present work, we investigate this new framework in the seismicity prior to the M_w6.3 Thessaly earthquake. Initially, we apply MRWA to the interevent time series of the successive regional earthquakes in order to investigate the approach of the regional seismicity at critical stages and to define the starting point of the natural time domain. Then, we apply the NT analysis, showing that the regional seismicity approached criticality a few days before the occurrence of the M_w6.3 earthquake, when the κ₁ natural time parameter reached the critical value of κ₁ = 0.070.

1. Introduction

On 3 March 2021, an M_W6.3 earthquake occurred in northern Thessaly (Central Greece), close to the cities of Tyrnavos, Elassona and Larisa (Figure 1). The earthquake occurred in a region that is one of the most seismically active in Greece, mainly characterized by normal faulting along NW–SE striking faults, which belong to the Thessaly fault zone [1,2,3,4,5,6,7,8,9,10]. Based on the provided focal plane solutions [11], the mainshock was generated by the activation of an NW–SE striking normal fault (Figure 1) [12]. The mainshock was widely felt in the Thessaly basin and in the surrounding areas, from Athens in the south to the northern borders of Greece.

The Thessaly basin has a well-known history of large earthquakes, with mainshocks presenting typical magnitudes between 6.0 and 7.0 [3,4,5,6,7,8,9,10,11,12,13,14,15]. During the last century, eight major earthquakes, with magnitudes equal to or larger than 6.0, have occurred in this area. The M = 7.0, 1954 Sofades earthquake was the most destructive event, resulting in heavy damages in the broader southern Thessaly region (Papazachos and Papazachou, 2002; [15]). The M = 6.8, 1957 Velestino earthquake and the M = 6.5, 1980 Almyros earthquake, are two more significant events, with similarly destructive consequences.

Strong arguments show that the earthquake generation process can be considered as a critical point phenomenon that culminates with a large event as some critical point is approached [16,17,18,19,20,21,22,23,24,25]. New findings regarding the complex dynamics that characterize various geodynamic phenomena illustrate stimulating features in the framework of new concepts, as that of non-extensive statistical physics [22,23,24,25], multiresolution wavelets analysis [26,27,28] and of the novel time domain, termed as natural time [29,30,31,32,33,34,35,36,37,38,39].

The concept of natural time (NT) has been introduced recently to analyze possible pre-seismic signals [29,30,35]. The analysis of various complex systems in the NT domain enables the optimal extraction of signal information by reducing the uncertainties related to the conventional time, as well as the identification of long-range correlations in the evolution of the system, even in the presence of “heavy tails” [40]. The usefulness of NT analysis has been discussed in a number of applications to known critical phenomena, such as fracturing, earthquakes, the 2-D Ising model and 3-D turbulent flow [24,35], and references therein, and it has been tested experimentally in fracturing experiments in the laboratory by analyzing acoustic emissions time series [23,41].

Furthermore, wavelet-based methods have been introduced to characterize fractal signals [42,43,44,45] and to overcome effects associated with non-stationarities [46,47], a very frequent effect in the time dynamics of an earthquake sequence.

The goal of the present work is to test and evaluate the seismicity patterns in terms of MRWA and NT analyses, as applied in the evolution of seismicity prior to the recent M_w6.3 Thessaly strong event. The recent upgrading of the regional seismological networks [48] provides an accurate catalogue of microseismicity in the area and enables the application of such methodologies. The earthquake catalogs used herein are taken from the Hellenic Unified Seismological Network (HUSN) (http://eida.gein.noa.gr/, last accessed on 27 May 2021), where instruments belonging to the HL (National Observatory of Athens, Institute of Geodynamics, 1997) [49] and HT (Aristotle University of Thessaloniki Seismological Network, 1981) [50] networks provide a complete spatial coverage in the broader area of Greece, with a magnitude of completeness (Mc) down to 2.0. Figure 2 presents the seismic activity observed in the region of Thessaly for a period starting in January 2016, approximately 1900 days before the 3 March 2021 mainshock and within an area of radius, R = 80 km, around its epicenter. The main objective of this study is to investigate the applicability of NT analysis, as presented in the regional seismic activity prior to the M_w6.3 Thessaly earthquake, integrated with the results of MRWA applied to the interevent time series of the successive events, in order to define, with an objective technique, the starting point for the analysis in the NT domain. The description of seismicity evolution with the NT parameters, integrated with the results of MRWA, represents a novel framework that may lead to a better understanding of the evolution of earthquake generation processes.

2. Multiresolution Wavelets Analysis in the Seismicity of Thessaly Region

The temporal evolution of seismicity and the time-scaling properties are of crucial importance [51,52,53] for understanding the correlation properties of seismicity [54]. The analysis of time intervals between successive seismic events can be grouped in exponential or power laws revealing similar behaviours over different scales [55].

A Wavelet Transform (WT) involves the decomposition of a signal function into simpler, fixed building blocks at different scales and positions. In a similar way to Fourier transform (FT), WT operates on a signal and transforms it from the time domain to a different domain, offering a new representation. In Fourier analysis, the sine and cosine functions are used as a basis and localized in the frequency domain with a difficulty to process a function having components that are localized in the time domain. As a result, a small frequency change in FT produces changes everywhere in this domain. On the other hand, wavelet functions are localized both in frequency or scale, and in time, via dilations and translations of the mother wavelet, respectively. This leads to compact representation of large classes of functions in the wavelet domain. This is one of the major advantages of WT: an event can be simultaneously described in the frequency domain as well as in the time domain, unlike the usual Fourier transform where an event is accurately described either in the frequency or in the time domain. As a consequence, MRWA of data with different behaviour on different scales can greatly benefit from the use of WT.

The discrete wavelet transform (DWT) transforms a data vector of length M into a different vector of the same length. A wavelet basis is characterized by a particular set of parameters, called wavelet filter coefficients. In practice, DWT is commonly implemented using dyadic multirate filter banks, which divide the signal frequency band into sub bands [28]. For a point process such as that of the interevent times sequence, the wavelet coefficients can be derived from

$$W_{m,n}^{wav}=2^{-m/2}{\displaystyle \sum }_{i=1}^L{t}_i\psi (2^{-m}i-n)$$

(1)

where the scale variable m and the translation variable n are integers, L represents the total number of interevent times t_i analysed and ψ is the wavelet function. The DWT is evaluated at the points (m, n) in the scale-interval-number plane. Smaller scales correspond to more rapid variations and, therefore, to higher frequencies.

In the current study, we perform MRWA by examining the standard deviation of wavelet coefficients as a function of scale, as described from

$${\sigma }_{wav}\left(m\right)=\sqrt{\frac{1}{N-1}{\displaystyle \sum }_{n=1}^N{\left(W_{m,n}^{wav}-\langle W_{m,n}^{wav}\rangle \right)}^2}$$

(2)

where N is the number of wavelet coefficients at a given scale m and the brackets indicate the average among the coefficients at a scale m. Figure 3 presents the interevent times between two successive events versus the occurrence time of the second event until the major seismic event, in a region within a circle of radius, R = 80 km, around the epicenter and a magnitude threshold, M_th = 2.0. The time period that was covered for MRWA of interevent times spanned from January 2016 until 3 March 2021, when the main event of M_w6.3 occurred. In Figure 4 the time– earthquake magnitude plot for a radius R = 80 km around the epicenter and a magnitude threshold, M_th = 2.0 is presented.

The wavelet choice is dictated by the requirement to identify a rather sharp change in a possible cyclic sequence. Following the same pre-processing approach as in [28], we tested several candidate wavelets at small scales up to m = 4 and received quite similar results. Thus, in the current work, we present results from the analysis using the db4 wavelet.

We investigated the time evolution of the σ_wav(m), using fixed event number windows of 16 events shifting through the entire series. The shift between successive windows was set in two events. Consistently with the length of the time window, we analysed the time variation of the σ_wav(m) for lower scales (m = 1 to 4) since the number of available events is limited. Each calculated value is associated with the time of the last event in the window. Figure 5 shows a representative set of results for the time evolution of the σ_wav(m) using the db4 wavelet with four scales for MRWA, for the seismicity observed in three circles around the epicenter of the mainshock and within a radius of 30 km, 50 km and 80 km, respectively.

An initial comment from Figure 5 is the significant temporal variability in the strength of the multiscale properties of the interevent times. As observed in previous studies [26,27], before the major event of the seismic sequence a traceable decrease in the temporal evolution of the σ_wav, m(t) appeared, especially at lower scales. Plots at Figure 5 dictate the search for a time marker beginning several months before the major event for all the scales analyzed. The sharp decrease at lower scales (m = 1 and m = 2), which is observed before the major event, can be qualified as such a time marker since the decrease is evident for several days and is clearly identifiable.

Translating the result from lower scales in an alternative way, we propose the use of the observed time markers, which appear a few months before the major event, as the initiation point for the natural time analysis that follows. The main purpose is to combine the two independent methods (MRWA and NT analysis) that have been successfully used for the identification of critical stages in earthquake preparation processes, in a joint approach that will maximize the advantages of each one. More specifically, the initial application of MRWA in a broader time period reveals time segments where the NT analysis is then used to investigate for indicators suggesting the entrance to the critical stage.

3. Natural Time Analysis of Seismicity before the Thessaly Mw6.3, March 2021 Earthquake

The analysis of a complex system in the NT domain has been introduced in [29,35]. In the case of seismicity, the natural time χ, defined as ${\chi }_k=k/N$, serves as an index for the occurrence of the k^th event out of N total events. The seismic moment released during the k^th event is then considered, forming the pair (χ_k, M_k) for further analysis (see [28]). The evolution of (χ_k, M_k) is further described by the continuous function F(ω), defined as: $F\left(\omega \right)={\displaystyle \sum }_{k=1}^N{M}_{k}exp\left(i\omega \frac{k}{N}\right)$ (3) where $\omega =2\pi \varphi$ and $\varphi$ stands for the natural frequency.

F(ω) is normalized by division with F(0)

$$\mathsf{\Phi }\left(\omega \right)=\frac{{{\displaystyle \sum }}_{k=1}^N{M}_{k}exp\left(i\omega \frac{k}{N}\right)}{{{\displaystyle \sum }}_{n=1}^N{M}_n}={\displaystyle \sum }_{k=1}^N{p}_{k}exp\left(i\omega \frac{k}{N}\right)$$

(3)

where $p_k=M_k/{\displaystyle \sum }_{n=1}^N{M}_n$. The quantity p_k describes the probability to observe an earthquake event at natural time χ_k. The normalized power spectrum can then be obtained from (4), as $\mathsf{\Pi }\left(\omega \right)=|\mathsf{\Phi }\left(\omega \right){|}^2$. In the context of probability theory, and for natural frequencies of $\varphi$ less than 0.5, Π(ω) reduces to a characteristic function for the probability distribution p_k. It has been shown that the following relation holds [29,56]

$$\mathsf{\Pi }\left(\omega \right)=\frac{18}{5{\omega }^2}-\frac{6\cos \omega }{5{\omega }^2}-\frac{12\sin \omega }{5{\omega }^3}$$

(4)

According to probability theory, once the behavior of the characteristic function of the distribution is known around zero, then the moments and, hence, the distribution itself can be approximately determined. For ω → 0, (4) leads to

$$\mathsf{\Pi }\left(\mathsf{\omega }\right)\approx 1-{\kappa }_1{\omega }^2$$

(5)

where κ₁ is the variance in natural time, given as

$${\kappa }_1=\langle {\chi }^2\rangle -{\langle \chi \rangle }^2={\displaystyle \sum }_{k=1}^N{p}_k{\chi }_k^2-{\left({\displaystyle \sum }_{k=1}^N{p}_k{\chi }_k\right)}^2$$

(6)

It has been shown that the properties of Π(ω) at ω → 0, i.e., the values of κ₁ = 0.07, can signify the approach of a complex system towards some critical point [35], such as that of an impending large earthquake (see [24,31,37] and references therein). Figure 6 shows an earthquake timeseries in conventional (Figure 6a) and natural time (Figure 6b) domains. Figure 6c shows the corresponding power spectrum Π(ω) for the critical stage with κ₁ = 0.070, based on Equation (6), while the two other curves are for non-critical stages. Theoretically, it has been shown that κ₁ approaches 0.083 as N → ∞, when there are no long-ranged correlations in the system [35].

As a new event occurs, the pair (χ_k, p_k) is rescaled and κ₁ varies. It has been verified that when the parameter κ₁ converges to the value 0.070, the system enters a critical state [33,35,56].

Furthermore, the entropy in the NT domain, S_nt, is defined as [35]

$$S_{nt}=\langle \chi ln\chi \rangle -\langle \chi \rangle ln\langle \chi \rangle ={\displaystyle \sum }_{k=1}^N{p}_k{\chi }_{k}ln{\chi }_{\kappa }-({\displaystyle \sum }_{\kappa =1}^{{\rm N}}{p}_k{\chi }_k) ln({\displaystyle \sum }_{\kappa =1}^{{\rm N}}{p}_k{\chi }_k)$$

where $\langle f\left(\chi \right)\rangle ={\displaystyle \sum }_{\kappa =1}^{{\rm N}}{p}_{k}f({\chi }_k)$.

The entropy, S_nt, is a dynamic quantity that depends on the sequential order of events. Moreover, upon the time reversal T, i.e., Tp_m = p_{N − m + 1}, the entropy, S_nt−, is further defined. When the analysed seismicity approaches a “true” critical state, the following conditions should be fulfilled [35,56]:

(i).

The “average” distance D, defined by the normalized power spectra Π(ω) of the evolving seismicity and by the theoretical estimation of Π(ω) for κ₁ = 0.070, should be less than 10⁻².

(ii).

The parameter κ₁ should approach the critical value of κ₁ = 0.070 by “descending from above”.

(iii).

Both natural time entropies, S_nt and S_nt−, should be lower than the entropy of uniform noise S_u = (ln2/2) − 1/4 when κ₁ approaches 0.070.

(iv).

Since the dynamic evolution of the system is expected to be self-similar in the critical state, the time of the true coincidence should not vary upon changing (within reasonable limits) either the magnitude threshold, M_th, or the area used in the calculation.

In [28], the authors proposed the use of the time marker indicated by MRWA in the seismicity evolution before the major event as the initiation point for the NT analysis. In the frame of this approach, the two independent methods (MRWA and NT analysis) were integrated to identify the approach to the critical stage in the earthquake preparation process. In particular, the initial application of MRWA in a broader time period of the regional seismicity before the major event reveals time segments where the NT analysis is going to investigate for indicators suggesting the entrance to the critical stage.

In Figure 7, the computed Π(ϕ) curves are shown as they approach the critical Π(ϕ) in the regional seismicity, for a threshold magnitude of M_th = 2.0 and for areas of radius R = 30 km, R = 50 km and R = 80 km, respectively, around the epicenter of the main event. This analysis clearly demonstrates that, from 28 February to 2 March 2021, one to three days before the M_w6.3 earthquake of 3 March 2021, the critical Π(ϕ) was approached. In all cases, for R = 30 km, R = 50 km and R = 80 km, the NT analysis starts at approximately one to two months around the corresponding time markers indicated by MRWA (see Figure 5). It may, thus, be considered that the critical point for the regional seismicity was approached around that time. What happened during the last time period before the main event can be seen in Figure 7, which depicts the time evolution of Π(ϕ), for 0 ≤ ϕ ≤ 0.5, for M_th ≥ 2.0 and R = 30 km, 50 km and 80 km, when calculations started on 6 August 2019 for R = 30 km and R = 50 km and on 7 November 2019 for R = 80 km. It also becomes interesting that around that time, when the critical point was reached, seismicity started to occur in the epicentral region registering a few shallow, weak earthquakes prior to the mainshock. These events, which may be considered as foreshocks, do not affect the analysis, as their occurrence coincides with the approach to the critical point.

Applying the NT analysis to the seismicity that occurred prior to the M_w6.3 event in the Thessaly region, starting the analysis from approximately the time markers in the lower scales indicated by MRWA and up to the time of the mainshock occurrence, we observe that all criticality requirements are fulfilled. The latter is more clearly demonstrated by the parameters D, κ₁, S_nt and S_nt−, as they evolved event by event, and are computed and plotted in the natural time domain and in the conventional time, approximately 100 days before the mainshock (Figure 8). In Figure 8, we observe that all the requirements are fulfilled a few days before the mainshock for all three cases that we study, i.e., for R = 30 km, R = 50 km and R = 80 km around the epicenter of the mainshock. The results, thus, indicate that the regional seismicity presented criticality characteristics a few days before the main event.

4. Concluding Remarks

In the present work, we investigated the regional patterns of seismicity in the area of the Thessaly (M_w6.3) strong earthquake on 3 March 2021, by applying MRWA and NT analysis, two methods that have been used for the identification of critical stages in the preparation process of major earthquakes. The analysis was performed in the natural time domain, with an approximate starting point indicated by MRWA. The latter showed a decrease in the standard deviation of the wavelet coefficients σ_wav(m) at much lower scales, similar to the observations in [26,27,54] prior to the occurrence of major events. Within this joint approach, the initial application of MRWA in regional seismicity around the epicenter, and for a wide time period before the mainshock, indicated a time segment where the NT analysis was applied in order to explore possible indicators that suggested the entrance to a critical stage.

The results demonstrated that regional seismicity approached criticality a few days before the M_w6.3 earthquake that occurred on 3 March 2021 in the Thessaly region, in agreement with the results in [57]. In other words, the curve of the power spectrum, Π(ϕ), in the natural time domain that characterizes the evolution of the regional seismicity, coincided with the theoretical curve of critical point phenomena a few days before the M_w6.3 mainshock, in a similar way to that of non-equilibrium critical systems. Hence, the analysis of the regional seismicity in the natural time domain, initiated at approximately the time marks indicated by the results of MRWA, pointed to an approximate date of the impending large M_w6.3 earthquake of 3 March 2021, within a narrow time window in the order of a few days. These results lay further support to the methodology introduced in [28] regarding the combination of MRWA and NT analyses for the identification of critical stages of regional seismicity prior to strong earthquakes, providing a novel and promising framework for better understanding the evolution of earthquake generation processes.