Multifractal fluctuations in earthquake-related geoelectrical signals

Multifractal fluctuations in the time dynamics of geoelectrical data, recorded in a seismic area of southern Italy, have been analysed. We find that the multifractality of the signal depends mostly on the different long-range properties for small and large fluctuations. Furthermore, we quantitatively characterized the multifractality of the geoelectrical time series, on the basis of the characteristic parameters derived from the multifractal spectrum (maximum α0, asymmetry B and width W) (Shimizu Y, Thurner S and Ehrenberger K 2002 Fractals 10 103). We then analysed the time evolution of the multifractal parameters and we found that the multifractal degree of the signal, revealed by the variation of the width of the multifractal spectrum, is enhanced in association with the occurrence of the largest seismic event. This study aims to suggest another approach to investigate the complex dynamics of geoelectrical signals.


Introduction
The search for correlations between variability of geophysical signals and seismicity represents one of the most pervasive research fields in the context of studies devoted to tectonic processes. Indirect information on the dynamics underlying the tectonic processes can be derived from the recording of the variations of geophysical and geochemical parameters. Among these, geoelectrical signals might be very useful to monitor and understand a collection of seemingly complex phenomena related to earthquakes [1]. As an example, variations in the stress and fluidflow fields can cause changes in the geoelectrical field, in resistivity, and in other electrical variables [2], so that the investigation of these induced fluctuations can provide information on the driving mechanisms during seismic activity.
The use of monofractal methods to extract quantitative information from earthquake-related geophysical signals is well known [3]- [6]. Monofractals are homogeneous objects, in the sense that they have the same scaling properties, characterized by a single singularity exponent [7]. Generally, there exist many observational signals which do not present a simple monofractal scaling behaviour. The need for more than one scaling exponent can derive from the existence of a crossover timescale, which separates regimes with different scaling behaviours [8]. Different scaling exponents could be required for different segments of the same time series, indicating a time variation of the scaling behaviour. Furthermore, different scaling exponents can be revealed for many interwoven fractal subsets of the time series [9]; in this case the process is not a monofractal but multifractal. Thus, multifractals are intrinsically more complex and inhomogeneous than monofractals [10] and characterize systems featured by very irregular dynamics, with sudden and intense bursts of high-frequency fluctuations [11]. The simplest type of multifractal analysis is given by the standard partition function multifractal formalism, developed to characterize multifractality in stationary measures [12]. This method does not correctly estimate the multifractal behaviour of signal affected by trends or non-stationarities. Another multifractal method was the so-called wavelet transform modulus maxima (WTMM) method [13], based on the wavelet analysis and it was able to characterize the multifractality of non-stationary signals. This method involves tracing the maxima lines in the continuous wavelet transform over all scales, implying a big effort in programming. Alternatively, a method based on the generalization of the detrended fluctuation analysis (DFA) has been developed [14]. This method, called multifractal-DFA (MF-DFA) is based on the DFA [15], and is able to reliably determine the multifractal scaling behaviour of non-stationary signals, without requiring a big effort in programming. Furthermore, the MF-DFA allows a global detection of multifractal behaviour, while the WTMM method is suited for the local characterization of the scaling properties of signals [16]. In a recent paper [17], a comparative study of applicability of the MF-DFA and the WTMM methods in proper detection of mono-and multifractal character of data has shown that WTMM overestimates the random multifractal component of the data (i.e. related to the probability distribution of the signal fluctuations), which is especially strikingly evident in the case of monofractal signals, while the MF-DFA performance is correlation-dependent; and the MF-DFA is recommended as a more reliable tool of the multifractal analysis in many situations.

Experimental data
We study a geoelectrical data set recorded at Giuliano station (40.688 • N, 15.789 • E), located in one of the most seismically active areas of southern Italy. The signal consists of voltage difference between two unpolarizable electrodes inserted at a depth of 1 m in the ground to avoid external temperature effects. The distance between the electrodes is 100 m. Figure 1 shows the time series, which consists of minute-sampled geoelectrical values, recorded from 1 March 2001 to 31 May 2003. It also shows the the earthquakes (with magnitude M 3.0) occurring during the observation period and satisfying Dobrovolskiy's rule [18]. This law, which is a theoretical relation between earthquake magnitude, distance from the epicentre and volumetric strain, states that detectable seismically induced strain exceeds 10 −8 . From this relation the maximum distance from the epicentre in which the effects of the earthquake are detectable is r = 10 0.43M , where r is measured in km. We can observe that the signal evidences large fluctuations in association with the occurrence of the selected earthquakes, more sharply in correspondence with the largest event (M = 4.1) of the seismic sequence.

Multifractal framework
Variations in the stress distribution around a focal area can lead to a complex structure for the geoelectrical data time series [4,5]. To analyse different aspects of the variability of its fluctuations and correlations at all experimentally observed temporal scales, one needs to perform diverse techniques of statistical analysis of the recorded data. Standard statistical analyses characterize variability in a very limited way. As an example, the power spectrum analysis is able to reveal scaling properties of only stationary signals, giving spurious results for non-stationary signals. DFA is a well-established method for determining data scaling behaviour in the presence of possible trends without knowing their origin and shape [19]. Its advantages over conventional methods (e.g. the power spectrum) are that it permits the detection of intrinsic dynamical features, like scaling and long-range correlations, embedded in a seemingly non-stationary time series, thus avoiding the spurious detection of apparent scaling, which may be an artefact of extrinsic trends [20].
The previous statistical methods are based on the assumption that the signal is monofractal, meaning that one single-scaling exponent is sufficient to describe its scaling properties. But observational signals could be multifractal, therefore they require more than one scaling exponent to be completely described.
Multifractals can be considered as a generalization of fractal geometry, essentially developed for the description of geometrical patterns [21]. Indeed, fractal geometry has been introduced to describe the scaling relationship between patterns and the measurement scale: the 'size' of the fractal object varies as the scale raised to a scaling exponent, given by the fractal dimension. The concept of multifractal leads to the existence of an infinite hierarchy of sets, each with its own fractal dimension. Therefore, multifractals require an infinite family of different exponents.
To perform a multifractal analysis for all values of q, Kantelhardt et al [14] have developed the MF-DFA, used to analyse the data presented in this paper and described in the next section.

MF-DFA
The main feature of multifractals is that they are characterized by high variability on a wide range of temporal or spatial scales, associated to intermittent fluctuations and long-range power-law correlations.
The MF-DFA operates on the time series x(i), where i = 1, 2, . . . , N and N is the length of the series. With x ave we indicate the mean value

5
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT We assume that x(i) are increments of a random walk process around the average x ave , thus the 'trajectory' or 'profile' is given by the integration of the signal Furthermore, the integration will reduce the level of measurement noise present in observational and finite records. Next, the integrated time series is divided into N S = int(N/s) non-overlapping segments of equal length s. Since the length N of the series is often not a multiple of the considered timescale s, a short part at the end of the profile y(i) may remain. In order not to disregard this part of the series, the same procedure is repeated starting from the opposite end. Thereby, 2N S segments are obtained altogether. Then we calculate the local trend for each of the 2N S segments by a least-squares fit of the series. Then we determine the variance for each segment ν, ν = 1, . . . , N S and for ν = N S + 1, . . . , 2N S . Here, y ν (i) is the fitting line in segment ν. Then, after detrending the series, we average over all segments to obtain the qth order fluctuation function where, in general, the index variable q can take any real value except zero. Repeating the procedure described above, for several timescales s, F q (s) will increase with increasing s. Then analysing log-log plots of F q (s) versus s for each value of q, we determine the scaling behaviour of the fluctuation functions. If the series x i is long-range power-law correlated, F q (s) increases for large values of s as a power law In general the exponent h(q) will depend on q. For stationary time series, h(2) is the well-defined Hurst exponent H [22]. Thus, we call h(q) the generalized Hurst exponent. Monofractal time series with compact support are characterized by h(q) independent of q. The different scaling of small and large fluctuations will yield a significant dependence of h(q) on q. For positive q, the segments ν with large variance (i.e. large deviation from the corresponding fit) will dominate the average F q (s). Therefore, if q is positive, h(q) describes the scaling behaviour of the segments with large fluctuations; and generally, large fluctuations are characterized by a smaller scaling exponent h(q) for multifractal time series. For negative q, the segments ν with small variance will dominate the average F q (s). Thus, for negative q values, the scaling exponent h(q) describes the 6 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT scaling behaviour of segments with small fluctuations, usually characterized by a larger scaling exponent.
The value h(0) corresponds to the limit h(q) for q → 0, and cannot be determined directly using the averaging procedure of equation (5) because of the diverging exponent. Instead, a logarithmic averaging procedure has to be employed It is worth noting that h(0) cannot be defined for time series with non-compact support, where h(q) diverges for q → 0. Concerning the orders q, it has shown the existence of a linearization effect in the behaviour of the estimators as a function of q: the estimated exponents account for scaling exponents only for values of q within an interval q ∈ [q − , q + ] [23,24]. If q is outside that interval, the estimates of the exponents systematically behave as a linear function of q. This linearization effect is neither related to infiniteness of moments (first-order phase transitions), nor to a finite size effect (firstand second-order phase transitions) nor to an estimation difficulty effect. The linearization effect is a limitation that is intrinsic to the nature of the multi-exponent multifractal processes [24].

The singularity spectrum
The multifractal scaling exponents h(q) defined in equation (7) are directly related to the scaling exponents τ(q) defined by the standard partition function multifractal formalism [14]. Suppose that the series x k is a stationary and normalized sequence. Then the detrending procedure of the MF-DFA is not required. Thus, the DFA can be replaced by the fluctuation analysis (FA), for which the variance is defined as Inserting this definition in equation (5) and using equation (7), we obtain For sake of simplicity, we assume that the length N of the sequence is an integer multiple of the scale s, obtaining N S = N/s and therefore The term [y(νs) − y((ν − 1)s)] is the sum of the numbers x k within each segment ν of size s. This sum is known as the box probability p s (ν) in the standard multifractal formalism for normalized series x k . The scaling exponent τ(q) is usually defined via the partition function Z q (s) where q is a real parameter as in the MF-DFA. Equation (11) is identical to equation (10), therefore In this equation h(q) is different from the generalized multifractal dimensions D(q) = (τ(q)/q − 1); in fact, while h(q) is independent of q for a monofractal series with compact support, D(q) depends on q in that case [14]. The assumption of compact support of the series leads to the fractal dimension of the support D(0) = −τ(0) = 1. Therefore, monofractal series with long-range correlations are characterized by linearly dependent q-order exponent τ(q), i.e. the exponents τ(q) of different moments q are linearly dependent on q with a single Hurst exponent Long-range correlated multifractal signals have a multiple Hurst exponent, i.e. the generalized Hurst exponent h(q), where τ(q) depends nonlinearly on q [25].
The singularity spectrum f(α) is related to τ(q) by means of the Legendre transform [26] where α is the Hölder exponent and f(α) indicates the dimension of the subset of the series that is characterized by α. The singularity spectrum quantifies in detail the long-range correlation properties of a time series [27]. The multifractal spectrum gives information about the relative importance of various fractal exponents present in the series. In particular, the width of the spectrum indicates the range of present exponents. To be able to make quantitative characterization of multifractal spectra, we fitted the spectrum to a quadratic function [28] around the position of its maximum at α 0 , i.e. f(a) = A(α − α 0 ) 2 + B(α − α 0 ) + C: the coefficients are obtained by an ordinary leastsquares procedure. In this fitting the additive constant C = f(α 0 ) = 1. Parameter B serves as an asymmetry parameter, which is zero for symmetric shapes, positive or negative for a leftor right-skewed (centred) shape, respectively. To obtain an estimate of the range of possible fractal exponents, we measured the width of the spectrum, extrapolating the fitted curve to zero [28]. The width of the spectrum was then defined as W = α 1 − α 2 with f(α 1 ) = f(α 2 ) = 0. Obviously, negative values for f(α) are permitted, and the proposed measure of the width of the multifractal spectrum is equivalent to the calculation of the range (maximum-minimum) of the h(q) scaling exponents. The measure of the width of multifractal spectra as a measure of degree of 8 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT multifractality has also been suggested in relation to the investigation of past climate variations [29]. Another measure of multifractality degree has been proposed [9], with the estimate of the Lévy index α in the power-law fitting function of the scaling exponent function K(q) for universal conservative multifractals. Both the measures are not in opposition, but the measure of the multifractality degree by means of the width of the singularity spectrum takes into account also also the negative moment orders q, and, thus, allows the capture of the information deriving from large as well as small fluctuations, whose scaling properties in the MF-DFA approach are investigated with positive and negative q values respectively. With low α 0 , the process becomes correlated; for example if the process had the tendency to move upward in the past, it will move upward with a probability larger than 1/2 in the next time step. Roughly speaking, a small value of α 0 means that the underlying process is more regular in appearance. The width of the spectrum W is a measure of how wide the range of fractal exponents found in the signal is; and, thus, it measures the degree of multifractality of the series. The wider the range of possible fractal exponents, the 'richer' is the process in structure. Finally, the asymmetry parameter B captures the dominance of low-or high-fractal exponents with respect to the other. A rightskewed spectrum indicates relatively strongly weighted low-fractal exponents, and high ones for left-skewed shapes.

Results
The data examined in this paper present clear irregular dynamics (figure 1), characterized by sudden bursts of high-frequency fluctuations, which suggest performing a multifractal analysis, thus evidencing the presence of different scaling behaviours for different intensities of fluctuations. Furthermore, the signal appears non-stationary, and, for this reason, we applied the MF-DFA.
We calculated the fluctuation functions F q (s) for timescales s ranging from 5 × 10 2 min to N/4, where N is the total length of the series. The length of the series (N ∼ 1.2 × 10 6 ) allows us to consider the estimated exponents reliable. We calculated the fluctuation functions F q (s) for q ∈ [−10, 10], with 0.5 step. Figure 2 shows the MF-DFA fluctuation functions F q (s) only for q = −10, 0 and +10. The lower and higher value of the q-interval are chosen in accordance with the values estimated in the framework of first-and second-order multifractal phase transition [9], and with the values of the orders of moments generally investigated in power-law correlated time series [29]. The different slopes of the fluctuation curves indicate that small and large geoelectrical fluctuations scale differently. Figure 3(a) shows the q-dependence of the generalized Hurst exponent h(q) determined by fits in the regime 5 × 10 2 min < s < N/4. A measure of the degree of multifractality can be given by the range (maximum-minimum) ρ, which is 0.95 ± 0.03. This measure can reveal the degree of the strength of the multifractality of the series, if compared with the range of the generalized Hurst exponents h shuf (q) of the shuffled time series (r shuf ≈ 0.05). Figure 4 shows the nonlinear dependence of the τ function on q and the multifractal spectrum f(α).
In order to quantitatively characterize the multifractality of the signal, we calculated the three multifractal parameters as shown in section 3.3. In the case examined in this paper the signal is correlated (α 0 = 0), and has a 'rich' time structure denoted by an enough wide range of fractal exponents (W ∼ 1.22) and is characterized by a relative dominance of low-fractal exponents (B ∼ −0.66).  We investigated the type of multifractality that underlies the q-dependence of the generalized Hurst exponent. Generally, two different types of multifractality in time series can be discriminated: (i) due to a broad probability density function for the values of the time series and (ii) due to different long-range correlations for small and large fluctuations, that is the sets of the large and small signal values scale differently. Both of them need a multitude of scaling exponents for small and large fluctuations. The easiest way to discriminate between these two types of multifractality is by analysing the corresponding randomly shuffled series [14]. In the shuffling procedure, the values are put into random order, and although all correlations are destroyed, the probability density function remains unchanged. Hence, the shuffled series coming from multifractals of type (ii) will exhibit simple random behaviour with h shuf (q) = 0.5. While those coming from multifractals of type (i) will show h(q) = h shuf (q), since the multifractality depends on the probability density. If both types of multifractality characterize the time series, thus the shuffled series will show weaker multifractality than the original one. Figure 3(b) shows the results of the generalized Hurst exponents versus q, averaged over ten randomly shuffled versions of the original time series. The error bars delimit the 1 − σ range around the mean values, where σ indicates the standard deviation. The h shuf (q)-values range around 0.5, but with a slight q-dependence; this indicates that the most multifractality of the geoelectrical data is due to different long-range correlations for small and large fluctuations. A measure of the degree of multifractality can be given by the range (maximum-minimum) of the generalized Hurst exponents ρ; the range of the shuffled series is ρ shuf ≈ 0.05, much smaller than the range of the original series, which can be considered as characterized by a strong multifractal degree.
We computed the singularity spectrum by the calculation of the partition function exponent τ(q), which is approximately linear versus h(q) (figure 4(a)). The singularity spectrum of the shuffled series is shown in figure 4(b); fitting with a quadratic function, we estimate the following values of the multifractal parameters: α 0,shuf ≈ 0.49 suggests the absence of correlations or the presence of very weak correlations; W shuf ≈ 0.19, as ρ shuf , indicates a very low degree of multifractality; the asymmetry B shuf ≈ −0.57 suggests that the shuffled series are characterized by a relative dominance of low-fractal exponents with respect to high ones. We investigated the variation of the multifractal behaviour of the series, analysing the variation with time of the three multifractal parameters: the maximum α 0 (t), the asymmetry B(t) and the width W(t), introduced in section 3.3. We have also calculated the time variation of the range ρ(t), that is the difference between the maximum and the minimum generalized h q exponent. We calculated the set of the generalized Hurst exponent {h q (t) : −10 q +10} in overlapping time windows 10 5 min long; the time shift between two successive windows was set to 5×10 3 min. By means of the Legendre transform, from the set {h q (t)}, we calculated the singularity spectrum and the multifractal parameters; each value was associated with the time of the last sample in the window. The window duration allows us to obtain reliable values of the exponents, estimated over enough long data sets; the time shift permits sufficient smoothing among the values. Figure 5 shows the time variation of the three parameters, derived from the singularity spectrum, and the range ρ, along with the average among the values obtained from ten shuffled series. The maximum α 0 ( figure 5(a)), the width W (figure 5(c)) and the range ρ ( figure 5(d)) of the original series are well above the corresponding parameters obtained from the shuffled series; in contrast, the asymmetry B ranges as well as the shuffled-B value ( figure 5(b)). Figure 6 shows the variation of the multifractal width W and the range ρ along with the earthquakes occurred in the area during the observation period; we observe a clear increase of both parameters in association with the largest earthquake of the sequence, indicating an enhancement of the multifractal degree during the preparation stage of the earthquake. However, since there are few rather large earthquakes (M > 3), the found correlation between the seismic events and the change in the multifractality of the electrical signal cannot be stated in an objective and statistically significant manner. But the peculiar variation of the degree of multifractality of geoelectrical signals suggest that also this feature should also be investigated in order to have a deep understanding of the geophysical process generating such signals.

Discussion
The analysis presented in the previous sections shows the multifractal nature of geoelectrical signals measured in seismic areas. This multifractality is due to the different long-range correlation behaviour characterizing small and large signal fluctuations. In particular, it has been shown that the degree of multifractality is characterized by a clear increase just prior to the largest earthquake of the sequence of seismic events occurred during the observation period. Nevertheless, there are other smaller enhancements that are not correlated with earthquake occurrence. Although earthquake precursory phenomena of a geoelectrical nature have been reported (in particular, transient variations of the Earth's electric field have been observed prior to seismic activities [30]; variations in the stress and fluid-flow fields can produce changes in the geoelectrical field [31]), a definite statement about the possibility to forecast future large earthquakes analysing the shape of the multifractal spectrum or the time patterns of its parameters cannot be given on the basis of a single data set. It should be necessary to analyse longer data sets, recorded in different tectonic settings.
However, the multifractality observed in the geoelectrical signals recorded in seismic areas can reflect the irregularity and heterogeneity of the crust, within which a phenomena generating geoelectrical fields occur. In fact, it is well known that the most relevant phenomenon that could originate the geoelectrical field is known as electrofiltration or streaming potential: an electrical signal is produced, when a fluid flows in a porous rock due to a pore pressure gradient. The phenomenon is generated by the formation within the porous ducts of a double electrical layer between the bounds of the solid, that absorbs electrolytic anions and cations distributed in a diffused layer near the boards. Due to the pressure gradient, the fluid flows and transports a part of the cations, giving on one side of the layer an excess of positive charges. This produces an induced electric field along the length of the duct and the associated potential differences (streaming potential). The streaming potential can be responsible for the voltage measures on the ground surface preceding an earthquake [32]. In a seismic focal region, the increasing accumulation of strain can cause dilatancy of rocks [33]. The phenomenon of dilatancy consists in the formation and propagation of cracks inside a rock as stress reaches a critical value [34]. If the rocks in the focal region and surrounding volumes are saturated with fluids, the voids generate pressure gradients. Hence, fluids invade the newly opened voids and flow until the pressure balances inside the whole system of interconnected pores. During the fluid invasion, the condition of rock hardening can be reached: the rock suddenly weakens and the earthquake is triggered. Therefore, the structure of the geoelectrical signal is linked to the structure of the seismic focal zone.
Furthermore, the enhancement of the multifractal degree before the occurrence of an earthquake could be related with the geometry and the structure of individual fault zones, which can be represented by a network with an anisotropic distribution of fracture orientations, and consisting of fault-related structures including small faults, fractures, veins and folds. This is a consequence of the roughness of the boundaries between each component and the interaction between the distinct components within the fault zone [35]. In fact, earthquake faulting is characterized by irregular rupture propagation and non-uniform distributions of rupture velocity, stress drop and co-seismic slip [36]. These observations indicate a non-uniform distribution of strength in the fault zone, whose geometry and mechanical heterogeneities are important factors to be considered in the prediction of strong motion. Experimental studies on the hierarchical nature of the processes underlying fault rupture, leading to the possibility of recognizing the final preparation stage before a large earthquake, have been performed [37].

Conclusions
The geophysical phenomenon underlying the earthquake-related geoelectrical variability is complex and is governed by physical laws that are not completely known. The multifractal analysis has led to a better understanding of such complexity, by means of the multifractal parameters, maximum, asymmetry, width and range. Using the surrogate series, obtained by a random shuffling procedure, the multifractality of the series is shown to be mainly due to different long-range correlations for small and large fluctuations. The singularity spectrum has led to a better description of the signal, revealing a clear enhancement of its multifractal degree (given by the variation of the multifractal width or the range) in association with the occurrence of the largest earthquake. Nevertheless, there are other smaller enhancements with no earthquake occurrence.
Of course, in order to be able to assess significant correlations between large earthquakes and patterns of multifractal parameters, and, in particular, the use of such patterns to perform feasible earthquake prediction, we need to investigate data sets longer than that examined in this paper.