Evolution of the multifractal parameters along different steps of a seismic activity. The example of Canterbury 2000–2018 (New Zealand)

Since the publication of Mandelbrot’s book (Mandelbrot, 1983), numerous studies have been devoted to characterize the self-similar (or scale-invariant) properties of a wide variety of natural phenomena by using the concepts of fractal geometry and fractal dimension (Enescu et al., 2005). The framework to characterize such non-linear dynamic systems seems to be better described by multifractals (Hirabayashi et al., 1992). In the case of the seismic activity of Earth, several studies ascertain that the temporal and spatial (epicentral and hypocentral distribution) behavior, as well as the energy distribution of earthquakes, has multifractal properties. For instance, the multifractal structure of the spatial distributions of earthquakes has been characterized in different areas of Earth, such as the Kanto region in Japan (Hirata and Imoto, 1991), Italy (Godano et al., 1996), the Himalayan region (Roy and Mondal, 2012), Iran (Zamani and Agh-Atabai, 2011), California, and Greece (Hirabayashi et al., 1992). Time distribution of events has also been analyzed by means of multifractal techniques in regions such as Romania (Enescu et al., 2005), central Italy (Telesca and Lapenna, 2006), and west India (Aggarwal et al., 2015), among others. In the first case (Enescu et al., 2005), whereas the time distribution of seismicity at short scale (minutes–hours) depicted a heterogeneous multifractal pattern, possibly caused by the aftershocks, at large scale, the time distribution of magnitudes changed toward monofractality. Conversely, in the second case (Telesca and Lapenna, 2006), the analysis of the elapsed times between consecutive seismic events permitted detection of a loss of multifractality along the aftershock sequences of every mainshock. In the third case (Aggarwal et al., 2015), the authors of the research observed an increase in the seismic magnitude multifractality along the aftershock sequences. More recently, Lana et al. (2020) have analyzed the mono/multifractal behavior of distances between consecutive aftershocks of a high magnitude mainshock, taking the seismic crisis of Landers, Northridge, and Hector mines in Southern California as the database.

Nowadays, it is absolutely accepted by researchers on Earth sciences that the tectonic stresses, generated by displacements of plate tectonics due to large-scale convective phenomena in Earth’s mantle, are the cause of seismic activity. Very detailed analyses should be necessary for a better understanding of these complex non-linear mechanisms, being relevant for the application of the fractal theory (Turcotte, 1997), to improve the knowledge of this physical problem. The possibility of detecting changes on fractal parameters, assimilated to warnings of high magnitude earthquakes, is also another question to be carefully considered.

The multifractal detrended fluctuation (MF-DF) algorithm is widely used to characterize multifractal scaling properties of non-stationary time series (Kantelhardt et al., 2002). Nevertheless, depending on the analyzed problem, other algorithms and theories should be convenient or necessary. Two examples could illustrate these facts. From the point of view of seismology, when the spatial distribution of hypocenters is analyzed, processes proposed by Turcotte (1997) and Goltz (1997) would be appropriate. Other examples could be the analyses of turbulent data and Brownian signals, which should be better analyzed by means of wavelet theory (Muzy et al., 1991; 1994). In this paper, we apply the MF-DF algorithm to analyze the degree of complexity associated with the seismicity generated at the Canterbury region (New Zealand) along the years 2000–2018. In particular, we focus our analysis on two time series: the interevent hypocentral distance, ∆δ(t), and the interevent time, ∆τ(t), series, defined, respectively, as distances and time intervals between two consecutive seismic events. Whereas ∆δ(t) describes the spatial trajectory of consecutive earthquake hypocenters, ∆τ(t) measures the waiting time between two consecutive earthquakes. In both cases, the evolution of these series would depend on the tectonic stress release in previous earthquakes and the specific state of stress on the rupture points (hypocenters). At the same time, this stress release could notably depend on the magnitude of recent earthquakes. Consequently, the quantities of both time series could be notably different along a time interval of very moderate earthquake magnitude (background seismicity), close in time to a high magnitude earthquake (main event) with a high increase in tectonic stress, and after a main event (aftershock interval) with a restructuration of crustal Earth strain and tectonic stress close to the mainshock domain.

An illustrative example of the non-linear complexity of the seismic activity could be the analyses based on the self-organized criticality (SOC) algorithms proposed by Bak et al. (1988; 1990) and applied by Scholz (1991), Correig et al. (1997), Turcotte et al. (2003), and Monterrubio et al. (2015), among others, to, respectively, analyze the mechanism of earthquakes and faults, rock fractures at micro- and macro-scales, and aftershock series.

As mentioned before, one the main objectives of this paper is to quantify the complexity of the seismic activity in a tectonic region of New Zealand by analyzing changes in the multifractal spectrum parameters along the series ∆δ(t) and ∆τ(t). Additionally, bearing in mind that the state of tectonic stresses and seismicity rates could notably change before an imminent mainshock of outstanding magnitude, ∆δ(t) and ∆τ(t) series are analyzed by means of the MF-DF algorithm with the aim of detecting possible warnings of forthcoming mainshocks. Concepts of multifractal theory were also applied years ago with the aim of analyzing the spatial distribution of hypocenters within tectonic domains of high magnitude mainshocks (Turcotte, 1997; Goltz, 1997). The algorithms used by these authors were different from the MF-DF used in this paper, given that the spatial density of hypocenters (a three dimensional problem) was analyzed instead of the ∆δ(t) and ∆τ(t) series (a one dimensional problem). Nevertheless, the same kind of multifractal parameters was successfully determined, being verified as the multifractal structure of the hypocenter’s spatial density.

From 2000 to 2018, a total of 15 889 events (Mw > 1) were recorded in the Canterbury region (New Zealand). This is a region with high seismic activity because plate boundary deformation creates a south-easterly advancing one, with the repetitive structural pattern governed by the propagation of northeast-striking thrust assemblages (Campbell et al., 2012). In 2010, the Canterbury earthquake, also known as the Darfield earthquake (Quigley et al., 2012), struck the south island of New Zealand with a moment magnitude of Mw = 7.2 at 4:35 a.m. local time on September 4. This event is the largest registered during the period 2000–2018. Additionally, three large events of magnitudes Mw = 6.2, Mw = 6.0, and Mw = 5.9 occurred in 2011 (February 22, July 13, and December 23, respectively). During 2010 and 2012, 28 events of magnitude 5 < Mw < 5.9 also shook this region. Figure 1 depicts the spatial distribution of seismicity in agreement with the seismic catalog obtained from GeoNet (GNS Science–Earthquake Commission, New Zealand; https://www.geonet.org.nz). The first event of the database was recorded on February 1, 2000, UTC 09:31:31.009, while the last event was detected on May 19, 2018, UTC 16:23:54.674, thus encompassing 6682 days of seismic activity. At the present study, the seismic activity analysis has been constrained from February 1, 2000 to February 22, 2016, covering [Figs. 1(a) and 1(b)] the square perimeter defined by coordinate intervals 171.8° E−172.9° E and 42.0° S−43.6° S.

Before the multifractal analysis, events equaling or exceeding Mc, the magnitude of catalog completeness (Popandopoulos et al., 2016) has been selected from the database. The frequency–magnitude distribution and Mc are computed by applying the Gutenberg–Richter law (Gutenberg and Richter, 1944),

with N being the cumulative number of earthquakes greater than a given magnitude Mw and parameters a and b depending on regional characteristics such as seismicity and stress field (Evernden, 1970; Ozturk, 2012). Definitively, the dataset used in this analysis consists of 9889 events with magnitudes larger than Mc = 2.5, fulfilling Eq. (1).

The evolution of the elapsed time and distances between consecutive earthquakes (exceeding Mc = 2.5), posteriorly analyzed in terms of multifractality, is described in Figs. 2(a) and 2(b). These evolutions are notably conditioned by the aftershock process after a mainshock. The elapsed time evolution is characterized by an evident increasing trend, up to the activation of a mainshock, and a sharp decrease immediately after. This sudden decrease has to be generated by the aftershock process following a mainshock, which is characterized by a high concentration of moderate and low magnitude earthquakes along the tectonic fault associated with the mentioned mainshock. From physics and geologic points of view, the spatial rearrangement and dissipation of tectonic stresses cumulated before the mainshock would be the cause of this aftershock activity. A relatively similar evolution is detected for the distances between consecutive events. Whereas the spatial distribution of consecutive earthquakes before a mainshock (basically background seismicity) depends on the spatial distribution of tectonic stress and geological faults, the seismic activity is notably concentrated on aftershocks, much more numerous along several months that the background earthquakes of the whole tectonic region are expected to be much more extended than the aftershock area.

The multifractal properties of nonstationary series are analyzed by means of the multifractal detrended fluctuation, MF-DF, analysis (Koscielny-Bunde et al., 1998; Talkner and Weber, 2000; Kantelhardt et al., 2002; and Shadkhoo and Jafari, 2009). The MF-DF methodology is described in detail by Kantelhardt et al. (2002). Assuming that {x_k} is a time series of length N, a brief overview of the algorithm steps is given below:

1. First step: Computation of the “profile” of the series:

   $Yi=∑k=1ixk−x,i=1,…,N$

   (2)

   〈x〉 is the average value of the series.

2. Second step: Division of the profile Y(i) into Ns = int(N/s) non-overlapping segments of equal length s. Given that the length N of the series is often not a multiple of the considered segment lengths, a short part at the end of the profile may remain. To not disregard this part of the series, the same procedure is repeated starting from the opposite end. Consequently, 2 Ns segments are obtained.

3. Third step: Computation of the local variance F²(s, ν) for every segment ν of length s by using a four order least-square polynomial fitting to obtain the differences between “profile” segments (first step) and the corresponding polynomial.

4. Fourth step: Computation of the q-order fluctuation function:

   $Fsq=12Ns∑12NslnF2s,νq/21/qν,q≠0,q∈R,$

   (3)

   $Fs0=14Ns∑12NslnF2s,ν,q=0.$

   (4)

   Given that it is necessary to repeat steps 2, 3, and 4 for several scales s, in agreement with (Kantelhardt et al., 2002), it is convenient that these scales vary within (m + 2, N/4), with m = 4 the chosen polynomial order (third step).

5. Fifth step: The q-order fluctuation function is expected to fit a power-law dependence on the segment length s:

   $Fsq≈shq,$

   (5)

   and h(q), the generalized Hurst exponent, can be well determined by a linear regression of ln{F(S)_q} vs ln(S).

In the case of non-stationary series, such as fractal Brownian signals, the exponent h(q = 2) will be larger than 1.0 and will satisfy h(2) = H + 1, where H is the well-known Hurst exponent (Movahed and Hermanis, 2008). For stationary time series, the value h(q = 2) is identical to the Hurst exponent. H > 0.5 indicates persistence in long-range correlation; H = 0.5 manifests the random character of the series, while H < 0.5 reflects anti-persistence. In the case of multifractal series, if positive values of q are considered, the segments ν with large variance (i.e., large deviations from the corresponding polynomial fit) will dominate the Fq(s) average. Thus, for positive values of q, h(q) corresponds to the scaling behavior of the segments with large fluctuations. For negative values of q, the segments ν with small variance F²(s, ν) will dominate the Fq(s) average, and h(q) then describes the scaling behavior of the segments with small fluctuations (Movahed and Hermanis, 2008; Burgueño et al., 2014).

The singularity spectrum, f(α), is related to the generalized Hurst exponent h(q) by means of the Legendre transform (Kantelhardt et al., 2002) as follows:

where α is the singularity strength or Hölder exponent and f(α) denotes the dimension of the subset of the series. The multifractal scaling exponent, also known as the mass exponent, is

and the Hölder exponent is defined as

The function f(α) describes the subset dimension of the series characterized by the same singularity strength α, with the singularity strength with the maximum spectrum designed as α₀. Small values of α₀ mean that the underlying process loses the fine-structure, that is, it becomes more regular in appearance; conversely, a large value of α₀ ensures higher complexity. The shape of f(α) may be fitted to a quadratic function around the position α₀,

The coefficient B manifests the asymmetry of the spectrum, being null for a symmetric spectrum. A right-skewed spectrum, B > 0, indicates a fine structure, while left-skewed shapes, B < 0, point to a smooth structure. The width of the spectrum, W, can be obtained by extrapolating the fitted curve f(α) to zero or, in other words, extrapolating the multifractal spectrum to q → ±∞. The spectral amplitude is defined as

with f(α_max) = f(α_min) = 0 and α_max (q → −∞) being larger than α_min (q → +∞). Given that q has been chosen ranging within the (−15, +15) interval, α_max and α_min have been obtained by numerically extrapolating Eq. (9) to f(α) = 0. The wider the range of possible fractal exponents α is, the stronger the multifractality would be. In other words, the wider the range of α is, the more complex the structure of the physical process (Burgueño et al., 2014) will be. A signal with a high value of α₀, a wide range of fractal exponents, and a right-skewed shape is more complex than a signal with opposite characteristics (Shimizu et al., 1982). The spectral width will be null for pure monofractal series, and h(q) will be independent of q. Therefore, there will be a unique value of α and f(α), with α being the Hurst exponent and f(α) being equal to 1.

The complexity index (CI) was proposed by (Shimizu et al., 1982) to quantify the degree of complexity, considering the parameters α₀, B, and W. The definition of this complexity index is based on the standardized values of α₀, B, and W,

where z(α₀), z(B), and z(W) are the corresponding standardized quantities. This standardized process is carried out taking into account the α₀, B, and W samples obtained for every moving window used to analyze the time evolution of the multifractality. Consequently, a sample of α₀, B and W is obtained for every one of the three multifractal parameters, being quantified for their corresponding averages and standard deviations. Finally, CI is defined as

where 〈Z〉 is the mean value and σ(Z) is its standard deviation. Values of Z ≥ 0 indicate high complexity, and those of Z < 0 indicate low complexity.

Given that the empiric samples of the multifractal spectrum f(α) sometimes do not fit well the quadratic function given by Eq. (9), the parameter B could be submitted to some uncertainties. With the aim of avoiding these computational uncertainties, the asymmetry of the multifractal spectra in terms of B has been substituted in this paper by the quotient (α_max − α₀)/(α₀− α_min), with this mathematical expression being an alternative version of the asymmetry. In agreement with the new definition, a quotient close to 1.0 implies a high symmetry. Conversely, quotients notably higher than 1.0 and lesser than 1.0 will imply right and left asymmetry, respectively.

The MF-DF algorithm is applied to quantify the complexity of the seismicity detected in the Canterbury region (New Zealand) during the period 2000–2018 and, more specifically, to analyze the multifractal behavior when a large event occurs. In this way, if changes in the evolution of some multifractal parameters occur close to a mainshock, they could be assumed as possible warnings of a forthcoming event. The multifractal spectrum f(α) is computed for the series of interevent hypocenter distances ∆δ(t) and interevent times ∆τ(t), for pairs of consecutive earthquakes. Both series describe the spatial and time behavior of seismicity since they give a measure of the degree of spatial and time earthquake correlation. A series of 9884 events (equaling or exceeding Mc = 2.5) is analyzed, by considering moving windows with a length of 1000 events (a long enough segment data for a right application of the MF-DF algorithm) and a shift of two events between consecutive windows. The total number of windows, Nw, required to cover all the events can be reduced to 4000, and the process described in Sec. III is used to compute f(α) and all the multifractal parameters for every moving window.

The window length is defined as an equal number of events, instead of equal time length, because the seismic activity, as expected, is not homogeneously distributed in time. It is worth mentioning that in the case of the first event (Darfield mainshock), given that a large enough number of events before the mainshock are not available, the accuracy of multifractal parameter determination is not sufficient to derive reliable results. Something similar occurs with the last event of Mw 5.7, given that the number of aftershocks with a magnitude larger than 2.5 suddenly decreases. In order to quantify the statistical significance of possible changes in the values of H(q = 2), α₀, W, α_max, and CI for the interevent ∆δ(t) and ∆τ(t) series, Buishand’s test (Buishand, 1982) is applied. This homogeneity test is able to statistically determine with a certain percentage of confidence if a series may be considered homogeneous along the time or, conversely, if a sharp change (a maximum for instance) or a time trend is detected. An example of successful application to climatic temperature series can be found in (Martínez et al., 2010).

The computational process of Buishand’s test can be described very briefly as follows:

Residuals of the original data series {x} are obtained by computing

If the original series is homogeneous, {X^*} will depict fluctuations around zero. On the contrary, a minimum or a maximum value of X^*(k) will be associated with positive or negative changes on {x} close to the critical point k of the series. The parameter R of Buishand’s test is defined as

for 1 ≤ k ≤ N, with σ(x) the sample standard deviation of the original series. Tables of the parameter R/N^1/2 for different percentages of the statistical significance of the detected sharp changes on {x} can be found in (Buishand, 1982; Wijngaard et al., 2003).

The analysis of the multifractal behavior of the seismicity is focused on the evolution of parameters H(q = 2), α₀, W, α_max, and CI, the last one including the contribution of α₀, W, and the asymmetry quantified by the quotient (α_max − α₀)/(α₀ − α_min), instead of parameter B. The evolution of these parameters along the 2000–2018 period is investigated for ∆δ(t) and ∆τ(t) series. Figure 3 shows some examples of multifractal spectra for three mowing windows including segments of consecutive 1000 elapsed times. The first example is characterized by a relevant symmetry, with signs of a left shift; and the second one again depicts quite similar characteristics but at the same time a notable increase on the spectral amplitude. The third example is characterized by notable symmetry with signs of a moderate right shift and remarkable small spectral amplitude. This last example also shows shortcomings sometimes appearing in the quantification of the multifractal parameters. In this last example, spectral amplitudes for the right side depart from a second order polynomial evolution, and in these cases, the asymmetry and spectral amplitudes are better quantified by taking into account extreme (α_max − α₀) and (α_min − α₀) values extrapolated to q tending to ±∞. These shortcomings could generate uncertainties or computational errors on the values of several multifractal parameters. Consequently, some of the relatively small oscillations in the parameters along the moving windows could not only depend on the complexities of the physical mechanism but would be generated by the just mentioned shortcomings. The numerous multifractal computations applied to the elapsed time and interevent distance series have also shown that sometimes, parameter q [Eqs. (3)–(5)] has not been possibly expanded up to ±15, being then again necessary for the just mentioned extrapolation to ±∞. Despite these shortcomings, the multifractal parameter maxima are well established and detected.

Figures 4–7 show the evolution of H, α₀, W, and α_max, respectively, for the set of moving windows. In every one of these figures, the upper plot shows the results for ∆δ(t) series and the lower plot, those for ∆τ(t) series. Vertical dashed lines indicate the first window including each one of the largest events detected (mainshocks of Mw = 6.2, 6.0, and 5.9).

Bearing in mind that the windows length is 1000 elements and the window shift is two elements, for a better interpretation of every parameter evolution, it has to be remembered that a mainshock is included along 500 moving windows. An example, available for the other figures, of the set of moving windows including a mainshock is shown in Fig. 4. Then, the maximum of the multifractal parameter is detected when the last earthquakes previous to the mainshock and a good part of aftershocks are included in the moving window. Consequently, these maxima would be correlated with the activation of the mainshock and its aftershock series, but they would not constitute a direct warning of a forthcoming high magnitude earthquake. These peaks are systematically more outstanding for ∆τ(t) than for ∆δ(t) series. This fact would reveal that the interevent time series have a more complex behavior than the interevent distance series close to and just after a large earthquake. Arrows point to peaks not related to any of the three largest events occurred in the analyzed period, thus indicating possible false detections of mainshocks. It is worthy mentioning that, as it can be observed in Fig. 4, the evolution of the Hurst exponent, H(q = 2), indicates a notable increase in persistence (long memory), especially for the ∆τ(t) series, when the activation of a mainshock is approaching. It is also relevant that only a few short periods of anti-persistence (short memory) are detected. A similar behavior to that observed for these four parameters, is also obtained for the complexity index, CI (Fig. 8). Forthcoming mainshocks of notable magnitude are accompanied by clear increases in the complexity index for both series and especially for ∆τ(t). Given that the parameter CI would represent a synthesis of the complexity degree of the evolution of elapsed time and distances between consecutive shocks, it could probably constitute the most relevant warning element to detect immediate earthquakes.

As mentioned before, some of the maxima detected for the multifractal parameters (indicated by arrows in the Figures) are not related to any of the four largest events included in the analyzed period and could then be interpreted as possible false detections of mainshocks. These peaks of multifractal parameters could be attributable to earthquake swarms, defined by space-temporal clusters of earthquakes with moderate magnitude, affecting a reduced domain and activated along a relatively short time. The three possible swarms are detected between the consecutive mainshocks of magnitudes Mw 6.2 and 6.0, 6.0, and 5.9 and after this mainshock magnitude are characterized by an almost equal magnitude range (Mw = 4.0–5.6, 4.0–5.8, and 4.4–5.5, respectively), areas of ∼12 × 23 km², 16 × 40 km², and 32 × 72 km² and depths ranging from 2.5 km–11.0 km, 4.5 km–10.5 km, and 4.5 km–18.5 km. Their lengths vary from 1 day to 11 days–36 days. To compare the main characteristics of the three mainshocks and swarms, both detected by the multifractal parameters, their details are summarized in Table I. Swarms used to be characterized by a short interval of moderate magnitudes along days, weeks, or even a few months and spatially distributed in areas with more extension than those corresponding to mainshocks and their aftershock series. Consequently, swarms and mainshocks are easily distinguished when seismic catalogs are analyzed a posteriori. A different question is the possibility of distinguishing mainshocks and swarms by means of the multifractal analysis, especially if a future objective is the predictability of mainshocks.

The assumed wrong detections appear at least one time for every multifractal parameter, being remarkable that the effects of the three swarms are observed more times for the parameters W, α_max, and CI corresponding to Δδ(t) than for those related to Δτ(t) series. Consequently, it could be concluded that interevent time series, Δτ(t), would be less sensitive to swarm interference and it would be then the most convenient series to be analyzed, in order to reduce the possible detection of wrong earthquakes. Nevertheless, two peaks appearing in Δτ(t) plots for parameters H and W, before the mainshock of magnitude Mw 6.2, must be also considered. Given that these peaks are not observed in the rest of the plots and possible associated swarms have not been detected in the seismic catalog, they could be the consequence of some computational irregularity concerning Eq. (3) for high positive or negative values of q.

To assess whether a multifractal peak is associated with a mainshock or with a wrong detection due to the effects of a swarm, Buishand’s homogeneity test becomes a very useful discriminating factor. According to Buishand’s test, significant peaks must be characterized by sharp minima of the Buishand parameter R/N^1/2, with N being the number of samples analyzed (number of moving windows) and R being the parameter quantified by Eq. (14). For all the multifractal parameters, peaks linked to the triggering of a mainshock have been found to be characterized by statistical significances equaling or exceeding 95% probability. Conversely, peaks associated with swarms have not reached such a high significance. As an example, Fig. 9 illustrates the evolution of the Buishand parameter for the complexity index, CI, and both time series. For the interevent distance series, only the swarm detected between mainshocks of Mw 6.2 and 6.0 is associated with a high significance (notably the negative Buishand parameter), which then leads to a detection of a false mainshock. By contrast, for the interevent time series, swarms are not linked to statistically significant minima of R/N^1/2. Consequently, the results derived for the CI would permit detection of, with notable certainty, the mainshocks of magnitude Mw 6.2 and 6.0 since they are associated with deep minima of the Buishand parameter, and, with less certainty, the Mw 5.9 mainshock (lower statistical significance). Beside this, a comparison of Buishand’s test results for both time series would permit removal of the possibility of false mainshock detections. Two swarms would not be linked to statistically significant minima of the Buishand parameter for any of the time series. Although the analysis of interevent distance series would lead to an erroneous detection of the third swarm as a mainshock, the lack of the corresponding minimum of the Buishand parameter for the interevent time series would permit discarding this wrong detection.

As a summary, the MF-DF analysis confirming multifractal behavior of ∆δ(t) and ∆τ(t) series is an appropriate methodology to quantify the complexity characterizing the seismicity of this New Zealand seismic active region. The results indicate that events with high magnitude are related to increases in H, α₀, α_max, W, and CI, with the multifractal parameters reaching peaks with statistical significance close to or exceeding 95%, in agreement with Buishand’s test. According to the evolution of multifractal parameters shown in Figs. 4–8, ∆τ(t) would be the most suitable series to detect high magnitude events, given that higher peaks of the parameters are more frequently reached for ∆τ(t) than for ∆δ(t). This different behavior for both series could be explained bearing in mind that close to the activation of a mainshock, specially at the beginning of the aftershock series, the elapsed time ∆τ(t) between consecutive events is notably shortened. Conversely, although the spatial density of hypocenters also increases, the distance between two consecutive events is not necessarily shortened, given that the previous (precursors) and posterior (aftershocks) seismic activity related to the mainshock event of high magnitude remains confined all along the area delimited by the high magnitude earthquake. Additionally, a lesser possibility of misleading detections of swarms as mainshocks has been found for ∆τ(t) than for ∆δ(t) series.

The MF-DF algorithm offers the possibility of analyzing the evolution of several multifractal parameters along periods of background seismicity, close to a mainshock and along the subsequent aftershock period. Consequently, the obtained results permit to confirm the multifractal structure of the tectonic processes generating the seismic activity on Earth, characterized by a notable degree of complexity and non-linearity, with some of these questions also quantified years ago by geologists and geophysicists from other points of view to those proposed in this paper.

With respect to the detection of large peaks related to swarms instead of mainshocks, the results obtained suggest that a comparison of statistically significant peaks derived from Δτ(t) and Δδ(t) series, together with the application of Buishand’s test, can distinguish between these two different kinds of earthquakes. In fact, the assumed detected swarms and mainshocks by means of multifractal analyses are absolutely coincident with the real episodes when the seismic catalog is revised a posteriori.

Another question to bear in mind is the high number of secondary maximum and minimum multifractal parameters along the moving window process. Despite the comparison between Δτ(t) and Δδ(t), as well as the application of Buishand’s test, permitted to detect and confirm swarms and mainshocks characterized by the main peaks, it should not be forgotten that these oscillations on the parameters could be consequence of a not absolute accomplishment of two relevant properties concerning multifractal structures of time series. One of them is the different long-range correlations of short and long time; the other is a broad density of probability for the empiric data analyzed. Additionally, the fact that the partial accomplishment of these two constraints could contribute to some incomplete multifractal spectra, with some of them shown in Fig. 3, should not be discarded, where sometimes the empiric multifractal spectrum cannot be expanded within the same range of positive and negative parameters q.

The detection of statistically significant peaks of these parameters close to a mainshock could also be a starting point for further research on warnings of high magnitude earthquakes. Given that the multifractal parameter peaks are detected after the mainshocks, except for some swarms, two questions should be considered for future analyses addressed to detect useful warnings of mainshocks. First, peaks achieved by the multifractal parameters of Δτ(t) and Δδ(t) series have to be generated by the activation of a mainshock (a sliding on a previous tectonic fault or Earth’s crust rupture, generating a new fault) and the high seismic activity (redistribution of tectonic stresses) at the beginning of the aftershock process. Second, strategies for detecting forthcoming mainshocks of notable magnitude would be based on the analysis of increasing statistically significant trends on the multifractal parameters before the detection of a maximum. In this way, the detected fluctuations of these parameters, with a relevant number of small secondary maxima and minima, would not mask the evolution toward a high magnitude event. The authors of this paper recognize that the evolution of some multifractal parameters before a mainshock, despite confirming the multifractal structure, does not facilitate this strategy. Some examples could be, for instance, the sudden increase detected on the multifractal amplitude, W, the maximum Hölder exponent, α_max, or the complexity index, CI. Nevertheless, the evolution of the Hurst exponent and the Hölder central exponent, α₀, would facilitate the detection of time trends associated with a forthcoming mainshock.

A first overview of Figs. 2(a) and 2(b) could lead to the wrong decision that plotting the time evolution of these two parameters would be sufficient to detect the imminent triggering of a mainshock of high magnitude. Despite notable sharp changes being quite coincident with a mainshock, two shortcomings have to be underlined. First, these changes are more probably due to the mentioned aftershock episodes than the main event. Second, the almost flat evolution or linear increase of these parameters before the mainshock is very different when the evolution is analyzed step by step. At a more detailed time scale, it is characterized by strong oscillations, with flat or linear increasing tendencies, making the algorithm for short-time prevention of triggering of high magnitude earthquakes very complex.

In short, the multifractal analysis permits achieving two results. First, a better description of the complexity, nonlinearity, and persistence/randomness characteristics of the seismic activity along background seismicity (low magnitude earthquakes), swarms of earthquakes with moderate magnitude, and mainshocks of high magnitude and their associated aftershocks can be achieved. Second, advances on the detection of multifractal parameters characteristics, which would permit establishment of warnings of forthcoming mainshocks of high magnitude, can be achieved.

The data that support the findings of this study are available from the corresponding author upon reasonable request and also from GeoNet (GNS Science–Earthquake Commission, New Zealand) (https://www.geonet.org.nz).

The research leading to these results has received funding from the European Union’s Horizon 2020 research and innovation program, under Grant Agreement No. 823844, the ChEESE CoE Project.