Fractal Patterns in Groundwater Radon Disturbances Prior to the Great M_w=7.9 Wenchuan Earthquake, China

1. Introduction

When it comes to severe natural disasters, earthquakes stand out since they may result in great loss of lives and property. Residents of large cities may experience severe effects from the massive quantity of energy produced during an earthquake, especially if the earthquake’s core is nearby. Catastrophic earthquakes are unavoidable as natural occurrences but highly challenging to foresee [1]. Consequently, the hunt for trustworthy seismic precursors is one of science’s greatest difficulties, and significant efforts have been done for many years in this area [2,3,4,5,6,7,8,9] The issue of earthquake prediction is still open, though [10]. Earthquakes are inherently complex, thus, several techniques and multi-level strategies are required [7]. In regions where severe earthquakes are possible, the associated prediction calls for the gradual downscaling of time, location, and magnitudes [2]. Along with the electromagnetic disturbances in the ULF, LF, HF, and VHF frequencies which indicate approaching earthquakes [4,9], radon-222 (hereafter radon) is a long-established precursor of impending seismic activity [2,7,11] . Radon is an inert gas created by the radioactive decay of the ²³⁸U series with a half-life of 3.86 days. When it disintegrates, it dissolves in the soil’s pores and liquids before moving on to surface, subsurface water and atmosphere. It is because that radon can travel to short and long distances in water and soil, and may be found at detectable quantities far from its generation location, that made it very significant in earthquake related studies [1]. As the above reviews mention [2,7,11], there is a great number of papers reporting pre-seismic changes of radon in groundwater, soil gas, atmosphere, wells and thermal spas. As a result, there is a substantial body of studies examining the relationship between radon emission and earthquake activity [10]. One of the key components of the present research is the abnormal behaviour of groundwater radon as a result of the catastrophic Wenchuan ($M_w$=7.9) shallow (depth=19 $km$) earthquake occurred on May 12, 2008 (calendar’s day - hereafter $day$ - 133), along the Longmenshan fault (31.0°N, 103.4°E) at Sichuan Province, China [12,13] (Figure 1)). This quake was the most destructive one of China since 1976 and the second most devastating seismic shock of this century after the great Sumatra earthquake in 2004 [14]. The groundwater data consists of one-year measurements received from China between 1 January 2018 and 31 December, 2008 by five different stations (Figure 2) between 105.6 $km$ and 526.0 $km$ from the earthquake’s epicentre (Table 1) a possible connection between the variability of geochemical signals and the Wenchuan-induced seismicity can provide insight into the associated underlying processes since the data include information about the subsurface dynamics. The Detrended Fluctuation Analysis (DFA), the use of fractal dimension (FD) with different methodologies and the power-law fractal analysis are very robust and reliable fractal methods [1,15,16,16]. Importantly, the fractal methods can recognise underlying scaling and long-range features that are characteristic for the unavoidable occurrence of earthquakes even in noisy and non-stationary time series.

2. Materials and methods

2.1. Experimental aspects

2.1.1. Earthquake activity

The Wenchuan earthquake occurred on the Longmenshan fault (LMSF), which is parallel to eastern Tibet and the Sichuan basin in the northeast and southwest directions (Figure 1) and is about 500 $km$ long and 40-50 $km$ broad [17]. According to Liu et al.[12], the Wenchuan-Maowen, Yingxiu-Beichuan and Guanxian-Jiangyou faults, as well as a number of thrust faults produced by the compression of the eastern Tibetan Plateau and the Yangtze craton, dominate the LMSF zone. Only one significant seismic event ($M_w$=6.1) occurred in the LMSF zone before the Wenchuan earthquake, and that was in 1989 [17]. The LMSF zone was dormant until 2008 [12]. The 290 $km$ long segment of the LMSF that was ruptured by the Wenchuan earthquake propagated independently 270 $km$ in the NE and 20 $km$ in the SW direction. The co-seismic surface rupture was 80 $km$. According to Liu et al. [12], the Wenchuan earthquake was caused by a shift in the LMSF’s dip angle (30° -50°SW to 60°-70°NE) and fault motion (SW thrusting motion to NE strike-slip motion). The Wenchuan earthquake disaster resulted in 69,225 fatalities, 374,640 injuries, 17,939 unaccounted-for persons, and over 5 million displaced people.

2.2. Measurement setup

Over the past 50 years, China has developed and established a structured network to track radon levels in groundwater for seismic studies. Almost all of China’ s provinces are represented in the network, which is made up of stations with a variety of instruments (Ren et al., 2012). The network is operated, maintained, and expanded with financial support from the China Earthquake Administration. The supplied radon series may be hourly (High Sampling Rate-HSR) or daily (Low Sampling Rate-LSR), depending on the sensors deployed at each station. Table 1 lists the LSR stations as Ganzi, Mingshan, Panzhihua and Sonpan, whereas the HSR station is only Kangding. To measure the amount of radon in groundwater, two instruments are used at China’s stations. Groundwater is sampled and degassed via bubbling at LSR stations. Then, using a specialised apparatus with the code FD-125, it is directed into an ion chamber or ZnS (Ag) detector where the concentration is determined by ionisations or scintillations. For the HSR Kangding station, groundwater is forced via a de-gasser and a gas-collecting device into a ZnS (Ag) detector, where radon is detected by a specialised instrument coded SD-3A with an hourly sample interval [18]. The accuracy for HSR and LSR measurements is 0.1 $Bq/L$.

3. Mathematical Aspects

3.1. Fractal and long-memory

Numerous natural physical systems can be explained by fractals. The usual fractal behaviour is often seen when the system or a section of it is stretched, translated, or rotated in space. Depending on the mathematics explaining the changes, the related system is either self-affine or self-similar. Natural systems that are self-affine and self-similar are fractals in the sense that any component of them is a little or big copy of the total, but at various scales. As a consequence, fractal systems may be studied by focusing on its constituent parts. Additionally, the scaling characteristics of a fractal system are closely related to its long-memory [19,20,21,22] and complexity [20,23,24,25,26]. Because fractal behaviour, long-memory and complexity are interconnected concepts, examining one of them will usually result in the other. For example, examining a system’s long-memory will help define its complexity and vice versa [1,27,28]. All of these characteristics can show whether a system’s past, present, and future are strongly connected. The direct methods are effective for calculating a system’s fractal features [16,27]. The methods of Katz, Higuchi, and Sevcik are used in this paper because they offer direct estimates for fractal dimension computations [1,29]. Power-law dependencies are also present in fractal systems with long memory and complexity. DFA and spectral power-law analysis are useful tools for properly delineating these connections [1,27]. These techniques are used in this paper, as well, because of this. By use of the associated Hurst exponent, all approaches may be compared. These techniques are going to be thoroughly explained in the sections that follow. Hurst exponent is introduced first, followed by the robust DFA, the techniques for direct calculation of fractal dimension and finally, the power-law analysis.

3.2. Hurst exponent

A metric known as the Hurst exponent (H) may be used to identify long-lasting connections in both time and space [30,31]. The roughness of the associated time series may be evaluated while time-evolving fractal phenomena can be identified using the Hurst exponent [1,32,33]. The Hurst exponent has been used in a variety of research topics, including hydrological applications [30,31], astrophysical applications [34,35], processes of capital markets [36,37,38,39], noisy observations of traces in traffic [40,41,42], seizures prior to epilepsy [43,44,45] climatic dynamics [46] and precursory time series before impending earthquakes [1,7,8,47] The Hurst exponent value offers important details about the time series [1,27,29,30,31,48,49,50,51]:

(i)

If $0.5<H\le 1$, the series has a positive long-range autocorrelation. A series’ high value is followed by a series’ high value, and vice versa. High Hurst exponents suggest persistent interactions that are predicted to occur in the series’ far future;

(ii)

If $0\le H<0.5$, low values follow high values in the time series, and vice versa. There is an ongoing exchange between low and high values for low H values in the time series’ future (this is known as anti-persistency);

(iii)

If $H=0.5$ associated processes are random and the time series are totally uncorrelated..

3.3. Detrended Fluctuation Analysis (DFA)

Long-range power-law connections as well as erratic oscillations observed in time series appear prior to earthquakes [1,15,52,53]. DFA is a reliable method for spotting long-range power-law relationships in noisy, brief, non-stationary signals [1,54]. This technique has been successfully applied in several scientific domains such as the study of changes in the weather and climate [55,56,57,58], DNA sequences [59,60], heart dynamics [61,62,63,64], urban-air-pollution [27,29,51], pre-earthquake recordings of radon in soil [65,66] and electromagnetic disturbances of the ULF, kHz and MHz ranges [1,67,68,69,70,71]. Theoretically, DFA can show if a temporal signal has concealed long-range linkages that result in a self-similar process. Calculating the scaling exponent of the integrated time series allows one to discover these long-term relationships in the original time series [1,29,47,56,60,61,72,73,74,75].

3.3.1. Application of DFA

The original time signal is initially integrated. The integrated signal’s fluctuations, $F\left(n\right)$, are then identified within a window of size n. The integrated time series’scaling exponent (self-similarity parameter), $\alpha$, is then calculated by fitting the linear $log\left(F\left(n\right)\right)-log\left(n\right)$ transformation with a least-squares function. The $log\left(F\left(n\right)\right)-log\left(n\right)$ line may display one crossover at a scale n where the slope exhibits an abrupt change, two crossovers at two different scales $n_1$, $n_2$ [47], or not even show a crossover at all, depending on the system dynamics. The following procedure may be used to produce the DFA of a temporal signal in a single dimension, $y_i$,($i=1,....,N$ ) [1,29,47]:

  (i)

The time series is, first, integrated:

$$y\left(k\right)=\sum _{i=1}^k\left(y\left(i\right)-\langle y\rangle \right)$$

(1)

The entire average value of the time series is denoted in Equation (1) by the symbols the symbol <...> and k stands for the various time scales.

 (ii)

The integrated time series, $y\left(k\right)$, is then separated into equal, non-overlapping bins of length, n.

(iii)

The function $y\left(k\right)$ that represents the trend in the box is then fitted. Simple linear trends or polynomials of order 2 or higher order may be used. Here, the linear function was used. This linear function’s y coordinate is denoted by the notation $y_n\left(k\right)$ in each box n.

(iv)

The local linear trend, $y_n\left(k\right)$, is then subtracted from the integrated time series $y\left(k\right)$, which is detrended in each box of length n. The detrended time series, $y_d^n\left(k\right)$, is determined in this manner and for each bin as follows:

$$y_d^n\left(k\right)=y\left(k\right)-y_n\left(k\right)$$

(2)

 (v)

The integrated and detrended time series’ fluctuations’ root-mean-square (rms) is then computed for each bin of size n as

$$F\left(n\right)=\sqrt{\frac{1}{N}{\sum }_{k=1}^N{\left\{y\left(k\right)-y_d^n\left(k\right)\right\}}^2}$$

(3)

where, $F\left(n\right)$ are the rms fluctuations of the detrended time series $y_d^n\left(k\right)$.

(vi)

For various sizes $\left(n\right)$ of the scale boxes, the method steps (i)–(v) are repeated. This reveals the specific sort of connection between $F\left(n\right)$ and n. If there are long-term relationships in the time series, then $F\left(n\right)$ and n have an exponential relationship.

$$F\left(n\right)\sim n^{\alpha }$$

(4)

The scaling exponent $\alpha$ (DFA exponent) in Equation (4) assesses the strength of the time series’ long-term relationships.

(vii)

The linear relationship between $logF\left(n\right)$ and $log\left(n\right)$, whose slope equals $\alpha$, is found by applying a logarithmic transformation to Equation (4). A strong linear connection suggests that the accompanying variations are long-lasting and, consequently, have a long memory. The square of the Spearman’s ($r^2$) is used in this paper to measure the accuracy of the linear fit. Good linear fits were defined as having $r^2\ge$ 0.95 or above [1,15,29,47,76].

3.3.2. Sliding window DFA

The following six-step procedure was followed in order to implement the sliding window DFA:

(a)

The time series was segmented in windows of 64 samples. This segmentation yields to approximately a two-month series’ part for the PZHS, SPS, GS, MSS LSR stations, which record one measurement per day. The 64 sample window was also employed for the fractal analysis of the data from three monitoring stations of urban air pollution with precisely the same measurement recording rate, namely, 1 measurement per day [29]. In a recent paper for the PZHS station, a 256 segmentation DFA was employed [65], whereas in radon in soil measurements an approach of 128-sample window was utilised [77]. Nevertheless, since the windows are shifted 1 sample forwards (sliding window technique), the whole signal is analysed, except from a 64 sample window, which was the final one. One the other hand, the 64-sample windowing yields to to a 64-hour window for the HSR station of KDS, i.e., an analysis of about 2.5-days. Despite this differentiation it is noteworthy that for a radon station in Pakistan, with the same recording rate as the one of KDS, a 64-window analysis was also utilised [16]. DFA from the data of KDS was analysed with 64 sample windows for consistency.

(b)

Every window was fitted using a least squares fit of$logF\left(n\right)$ vs $log\left(n\right)$ in accordance with Equation (4). The data were fitted to a straight line without seeking cross-overs, as in related literature [1,29,65], with the restriction that the slope of the fit display a square of Spearman’s correlation coefficient above or equal to 0.95.

(c)

The window was advanced by one sample, and the steps (A) and (B) were repeated until the signal’s end.

(d)

DFA slopes $\alpha$ were plotted against time, and the associated plot data were exported to ASCII output files for further use.

3.4. Fractal Dimension Analysis

3.4.1. Katz’s method

To determine the fractal dimension, D, according to Katz’s method, the transpose array ${[s_1,s_2,...,s_N]}^{⊺}$ of the series $s_i$, $i=1,2,...,N$ is first determined, where $s_i=(t_i,y_i)$ and $y_i$ are the measured series values at the time instances $t_i$ [78,79]. The value pairs ($t_i,y_i$) and ($t_{i+1},y_{i+1}$) correspond to two following points of the time series ($s_i$ and $s_{i+1}$), for which the Euclidean distance is:

$$dist(s_i,s_{i+1})=\sqrt{\left(t_i^2-t_{i+1}^2\right)+\left(y_i^2-y_{i+1}^2\right)}$$

(5)

The distances in Equation (7) add up in a curve whose total length is:

$$L=\sum _{i=1}^{i=N}dist(s_i,s_{i+1})$$

(6)

This curve will stretch in the planar to d, if it does not cross itself, where d is :

$$d=max\left(dist(s_i,s_{i+1})\right),i=2,3,...,N$$

(7)

By combining equations (5), (6) and (7), the Katz fractal dimension, D, becomes

$$D=\frac{log\left(n\right)}{log\left(n\right)+log(d/L)}$$

(8)

where $n=L/\overline{\mathrm{a}}$ and $\overline{\mathrm{a}}$ is the average value of the distances of the points.

3.4.2. Higuchi’s method

To calculate the fractal dimension, D, of a time series

$$y\left(1\right),y\left(2\right),y\left(3\right),...,y\left(N\right)$$

(9)

recorded at intervals $i=1,2...N$, a new sequence , $y_m^k$, is constructed as follows [80,81,82]:

$$y_m^k:y\left(m\right),y(m+k),y(m+2k),...,y(m+\left[\frac{N-m}{k}\right]k)$$

(10)

The length of the curve associated to the time series is given by [80]:

$$L_m\left(k\right)=\frac{1}{k}\left(\sum _{i=1}^{\left[\frac{N-m}{k}\right]}y(m+ik)-y(m+(i-1)k)\right)\left(\frac{N-1}{{\left[\frac{N-m}{k}\right]}^k}\right)$$

(11)

In both equations m and k are integers that specify the time interval between the series’ samples and are connected by the formula $m=1,2...k$, where $\left[...\right]$ is the Gauss notation, namely, the bigger integer part of the included value. By inserting the normalisation factor

$$\frac{N-1}{{\left[\frac{N-m}{k}\right]}^k}$$

(12)

the average value, $〈L\left(k\right)〉$, of the lengths of Equation (13), exhibits a power law of the form:

$$〈L\left(k\right)〉\propto k^{-D}$$

(13)

The Higuchi’ s fractal dimension, D, is finally calculated by the slope of the linear regression of logarithmic transformation of $〈L\left(k\right)〉$ versus k where $k=1,2,...,k_{max}$. It must be noted that the time intervals are $k=1,..,k_{max}$ for $k_{max}\le 4$, i.e., $k=1,2,3,4$, for $k_{max}=4$ and $k=\left[2^{(j-1)/4}\right]$, $j=11,12,13...,$ for $k>4$ ($k_{max}>4$). Again $\left[...\right]$ is the Gauss notation [79].

3.4.3. Sevcik’s method

The fractal dimension of a time series according to the method of Sevcik [83], is approximated from the Hausdorff dimension, $D_h$, as [79]

$${\displaystyle D_h=\underset{ϵ\to 0}{lim}\left[-\frac{log\left(N\right(ϵ\left)\right)}{log\left(ϵ\right)}\right]}$$

(14)

where $N\left(ϵ\right)$ is the number of segments of length $ϵ$ that add up to a curve that is associated with the time series. If the length of the curve is L, then $N\left(ϵ\right)=L/2ϵ$[79] and $D_h$ is:

$${\displaystyle D_h=\underset{ϵ\to 0}{lim}\left[-\frac{log\left(L\right)-log\left(2ϵ\right)}{log\left(ϵ\right)}\right]}$$

(15)

By applying a a linear transformation twice, the N points of the curve L can be corresponded to a unit square of $N\times N$ cells of the normalised metric space. With this transformation Equation (15) provides the Sevcik’s fractal dimension [79,83]:

$${\displaystyle D_h=\underset{N\to \infty }{lim}\left[1+\frac{log\left(L\right)-log\left(2ϵ\right)}{log\left(2\right(N-1\left)\right)}\right]}$$

(16)

The calculation improves as $N\to \infty$.

3.4.4. Computational methodology of Fractal Dimension

The next methodology was followed to calculate the fractal dimensions :

(i)

As in Section 3.3.2, the time series was segmented in windows of 64 samples. As mentioned, this segmentation corresponds for the LSR stations (PZHS, SPS, GS, MSS), approximately to a two-month signal. 64 sample windowing was employed as well in the fractal dimension calculation (with the same methods) from the data of the three LSR urban pollution stations with identical rate of measurements, i.e., 1 measurement per day [29]. As also mentioned in Section 3.3.2, for the HSR station KDS, the 64-sample segmentation, corresponds to approximately, 2.5-days. In a previous fractal dimension analysis (with the same methods) for the PZHS station was implemented with a 256 window [65], however, in a very recent fractal dimension analysis with identical methodology for a HSR radon station in Pakistan of the same rate of measurements as the one of KDS (one measurement per hour), the 64-window approach was utilised [16]. Finally, as in Section 3.3.2 and for consistency with the windowing of the other stations, the 64-sample window was chosen here as well for the KDS station.

(ii)

The fractal dimensions of each method were calculated as follows:

• For the Katz’s method : Fractal dimension is D of Equation (8) for n=64 and $\overline{\mathrm{a}}$=1 collected sample per measurement interval (1 day fo PZHS, SPS, GS and MSS- stations and 1 hour for the KDS station) since $\overline{\mathrm{a}}$ corresponds the distance between the points of the series which constitute the parameter L[1,16,29] .
• Higuchi’s method : Equal to the slope D of the first order least-square fit of the logarithmic transformation of Equation (13), namely the relation of $log\left(〈L\left(k\right)〉\right)$versus $log\left(k\right)$, for $kmax=16$. In the aforementioned analysis for the urban air pollution stations [29], $kmax=4$ whereas in the analysis of radon in Pakistan [16] and of the electromagnetic disturbances of the Ileia station, Greece, the approach $kmax=16$ was used. Based on the two latter papers, the $kmax=16$ was also selected here.
• Sevcik’s method : Equal to the Hausdorff dimension of Equation (16) ($D=D_h$) for N=64, namely equal to the number of samples in each window which constitutes parameter L.

(iii)

Each window was forwarded one sample (sliding window technique) and the procedure (i)-(ii) was iterated until the end of the time series.

(iv)

Time evolution plots of the fractal dimensions in accordance to the Katz’s, Higuchi’s and Sevcik’s methods were generated and the partial data were extracted to ASCII files for further use.

3.5. Power-law analysis

The fractal power-law approach is another robust technique to identify the hidden long-lasting trends that are connected with the long-term links between space and time and are detectable before earthquakes [1,15,47,68,76,84,85,86,87,88]. As with all fractal based techniques, power-law analysis describes the main core of fractality, namely the existence of a power-law. Another reason is because the earthquake-generating systems progress gradually to self-organised critical (SOC) states exhibiting fractal evolution in space and time [87]. There have been approaches of fractal power-law analysis based on the Fourier transform [87,88]. The advantageous use of wavelets instead of the Fourier transform has been pointed in several publications (e.g., [67,89,90,91]).

A time series’ power spectral density, $S\left(f\right)$, will follow a power-law, if the series is a temporal fractal.

$$S\left(f\right)=a·f^{-\beta }$$

(17)

In Equation (17) $\beta$ is the power-law exponent which quantifies the strength of he power-law connection, a is the amplification of the spectral density and f is a frequency of a transform. According to several publications, this transform was chosen to be the wavelet one based on the Morlet bases and specifically f to be the central frequency of the Morlet wavelet [1,15,67,76,84,89,90,91]. The logarithmic transform of Equation (17) gives:

$$logS\left(f\right)=loga+\beta ·logf$$

(18)

Since Equation (18) is a straight line, the values of $\beta$ and a can be determined by fitting the corresponding data with the least squares method. As with Section 3.3.2, the goodness of fit of the least-square fit is quantified by the square of the Spearman’s ($r^2$) coefficient under the constraint $r^2\ge$0.95. The technique has been also described in other publications as spectral fractal analysis or spectral power-law fractal analysis. Hereafter the phrase "power-law analysis" will be used from now on.

3.5.1. Computational methodology of power-law analysis

Following the logic of subSection 3.3.2 and subSection 3.5.1, the next steps were followed to implement the power-law analysis:

(a)

The time series was separated into 128-sample windows. This is the double window size in comparison to the other two methods. This is done because the power-law analysis does not work well with small sized windows and for this reason the 64-sample window yielded to false runs. For the LSR stations (PZHS, SPS, GS, MSS) this segmentation corresponds to a four-month signal and to a roughly 5 day signal for the KDS station. In previous publications, 128-sample window was employed in $\beta$ parameter estimations [77] for recordings of similar recording rate, while in others, a 512 sample window [91], however, with a recording rate of one measurement per 10 minutes.

(b)

The power spectrum, $S\left(f\right)$, based on the Morlet wavelet, as well as, the central Morlet frequency, f, were calculated in each segment.

(c)

Parameters $logS\left(f\right)$ and $logf$ were fitted via least-squares. Exponents $\beta$ and power amplification $\alpha$, were computed for every window by each fit , under the constraint, Spearman’s $r^2\ge 0.95$.

(d)

The steps (A) through (C) were iterated to the end of the time series. At each iteration, the the window was shifted on sample forwards. As with the other techniques, the whole time series was covered but the last window.

(e)

The $\beta$, $loga$ data were tabulated and saved in ASCII format for further use.

3.6. Further issues for chaos analysis

3.6.1. Segmentation to Chaos analysis classes

Two further classes were formed for additional analysis:

(a)

Class I: This class includes windows that, on the one hand, showed DFA least square log-log fits with Spearman’s coefficient of $r^2\ge 0.95$ and, on the other hand, had a scaling exponent for the DFA that was in the range of $1<\alpha <2$, i.e., they could be modelled by the fBm class [76]. The Class I segments:

• with differences between anti-persistency ($1.35<\alpha <1.5$ and persistency ($1.5\le \alpha <2$) are of some pre-seismic precursory value[15,47,76,84,91].
• with persistent behaviour ($1.5\le \alpha <2$) have been stated by investigators (e.g., [7,8], and references therein), as footprints of ensuing earthquakes.

(b)

Class II: This class includes time series windows with DFA textquotesingle’s $r^2<0.95$ (i.e., they do not adhere to the prominent fBm class) or $0<\alpha <1$ (i.e., they adhere to the fractional Gaussian noise (fGn) class). It is important that the Class II segments:

• have low predictability and, hence, they are of low precursory value [1,15,47,76,77,84,91].
• are the complement of the Class I ones.

3.7. Comparisons of the fractal results

According to previous publications [1,47,76,77], the results of the fractal methods can be compared to each other, but the best approach is to compare all results through the Hurst exponent.

For the Class I segments the Hurst exponent (H) is calculated from the fractal analysis parameters as follows (e.g., [1,76],and references therein):

(1)

From DFA’ $\alpha$-exponent as:

$$H=\alpha -1$$

(19)

(2)

From fractal dimension (D) as:

$$H=2-D$$

(20)

(Berry’s equation)

(3)

From power-law $\beta$ as:

$$H=0.5·(\beta -1)$$

(21)

It should be emphasized that variations from the straightforward linear connection of equations (19)–(21) are seen in the in situ data, in accordance with comprehensive argumentation [1,15,47,76]. The relationship between the fractal analysis parameters remains linear, possibly of a slightly different form, as indicated in the aforementioned works.

3.8. Meta-Analysis

The, so called, meta-analysis [1,51], is implemented by combining the outcomes from the ASCII files of all five methods, namely, DFA, Higuchi’s, Katz’s and Sevcik’s fractal dimensions and power-law analysis. A two-step process is followed:

(a)

(Step-1): According to user-defined thresholds, each ASCII output results file, is computationally scanned for out-of-threshold values. The ASCII files carrying the fractal dimension values are searched for under threshold values, whilst the ASCII files containing the DFA’s exponents and the power law $beta$-values are looked for over threshold values. New ASCII step-1 files are generated that contain the out-of-threshold values.

(b)

(Step 2): Under the restriction that each segment’s first sample date is arbitrarily considered as the date of the whole segment, the step-1 ASCII files of are computationally filtered to find areas with common dates. The above computational process, results in the full coverage of all dates except the one of the last window. The whole procedure is iterated in the results of all possible combinations of:

• DFA versus fractal analysis or versus at least two fractal dimension calculation techniques (6 combinations);
• Fractal analysis versus at least two fractal dimension calculation techniques (4 combinations);
• One fractal dimension calculation technique versus the other two (3 combinations);

Through the above repetitive process, 13 unique combinations of techniques per five, four, three and two are produced. This is very important because it is practically equivalent to the coupling of different mono-fractal methods and this fact provides a synthetic view of the results of the fractal techniques increasing, as a result, the scientific evidence regarding the underlying nature of the identified fractal disturbances.This has been pointed in recent publications [1,15,27,29,51,76].

4. Results and discussion

Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6 present the variation of the DFA exponent $\alpha$ over time in respect to the evolution of the associated square of the Spearman’s correlation coefficient versus the measured disturbances of radon in groundwater. As can be observed form these figures, the DFA scaling exponent profile is entirely distinct from the one of time series. This has been noted in earlier works as well [1,15,27,29,47,65,76,92]. The reason relies on the fact that DFA manages to locate effectively hidden forms in time series, even in non-stationary [52,54,93]. Numerous $\alpha$ exponents are within the Class-I range (Section 3.6.1). This means that the corresponding 64-sample windows are successful fBm ones ($r^2\ge$0.95) and this has been acknowledged as a notable sign of pre-seismic activity [15,47,76,77,84,91]. The specific details of each station shows interesting information. At first, Figure 2 is a noteworthy case of completely different DFA exponent and signal profiles, for the KDS station. By inspecting this figure, a first period can be observed starting from window 1 ($day$ 1) up to window 2100 (approximately $day$ 93). During this period the DFA scaling exponents $\alpha$ are enhanced with $1.5\le \alpha <2.0$ ($0.5\le H<1.0$, equation 19) whereas, most significantly, several exponents are above 1.8 ($0.8\le H<1.0$ , equation 19). All these are successful fBm segments since the corresponding Spearman’s coefficient in each window has $r^2\ge$0.95. These successful fBm segments exhibit persistency since Hurst exponents are above 0.5, while numerous from these have great persistent behaviour ($0.8\le H<1.0$). As mentioned in Section 3.6.1, the references there and those of this section, the successful fBm segments and especially those with great persistent behaviour, correspond to radon in groundwater areas with a high potential of association with seismic activity. In addition, in Figure 3 a very interesting period is observed between windows 50 and 100 (approximately between $day$ 60 and $day$ 121) with DFA exponents $1.5\le \alpha <2.0$ ($0.5\le H<1.0$) and several above 1.8 ($0.8\le H<1.0$), for the GS station. Similar is the case for the MSS station. The corresponding high DFA period ($1.5\le \alpha <2.0$, several exponents higher than 1.8) is between window 70 ($day$ 85) and 115 ($day$ 139). What makes these three cases of great significance is that the above periods with high DFA exponents are rather synchronous. They are parallel in time and there are periods common in these three stations. Indeed the KDS period is from $day$ 1 to $day$ 93, the one from GS between $day$ 60 and $day$ 121 and the period of MSS from $day$ 85 to $day$ 115. Importantly, these periods correspond to strong persistency since several Hurst exponents are in the range $0.8\le H<1.0$. They also correspond to the Class-I category, exhibiting strong fBm behaviour. The significance of identifying strong persistent and fBm behaviour has been emphasised in several publications [1,15,29,76,89,91] as a footprint of ensuing earthquakes. The interpretations of these precursory footprints in these publications are based to the asperity model [94], according to which it is the roughness of fBm profiles, their relative motion and their association with micro-crack branching and acceleration that explains why especially the persistent fBm behaviour is a great sign of precursory activity during preparation of earthquakes. These findings are very important and have to be emphasised, especially, because the DFA results of the other two stations (PZHS, SPS) have almost all areas with low DFA exponents, mostly in the Class-II category or at the lowest part of the Class-I category. This latter fact shows very clearly that the PZHS, SPS disturbances are of low precursory value, if not totally insignificant. From the above DFA evidence it can be supported from the aspect of the time, that the important time periods identified (KDS from $day$ 1 -$day$ 93, GS from $day$ 60 - $day$ 121 and MSS frrom $day$ 97 - $day$ 130), comprise pre-earthquake signs of the great Wenchuan earthquake occurred on $day$ 133. According to the literature ([2,8], see reviews), the KDS-GS-MSS time window is within the precursory window range of radon precursors. Especially, when the big magnitude ($M_w$=7.9) is accounted, the above time window extends well before the earthquake occurrence (even from $day$ 1).

The differentiation between the KDS, GS and MSS stations versus the PZHS and SPS stations becomes more important when the differentiations in the distance and underlying geology are taken into consideration. Under the aspect of distance, a potential claim for the PZHS station could be that its long distance (526.0 $km$) from the Wenchuan’s earthquake epicentre make it difficult to show significant DFA disturbances. This claim however becomes unsupported when it is accounted that the PZHS station has shown previously DFA variations of noteworthy precursory value for earthquakes with epicentres above 533 $km$[65]. In support, the GS station shows precursory DFA variations despite being 325.5 $km$ far from the epicentre, especially when the SPS station did not show precursory DFA results even located quite nearer (182.5 $km$) and in comparable distance with the other two stations (KDS and MSS), which both showed significant DFA outcomes.This a significant observation that has to be outlined. As has been pointed out in the above references [1,15,29,76,89,91], the preparation area of earthquakes includes special preferable precursory paths. Moreover, the selectivity effect that has been proposed and utilised for precursory activity (e.g., [95,96,97,98]) suggests that, during seismic preparation, there are selective paths that the disturbed activators follow. Following certain paths has been pointed for MHz electromagnetic variations [47,99] for electromagnetic precursors, radon precursors [47,100] and has been expressed in reviews of the subject [2,5,7,8,11]. In this sense the roadmap a from station to earthquake’s epicentre, gains special geophysical meaning. In this consensus , it is very important in relation that the GS and KDS stations are in the same big fault path. Even more importantly, the MSS station is also in the prolongation of this GS-KDS path and, surely, in line with the Wenchuan’s epicentre. It may be supported from the above, that the high elevation path, in association with the great fault path of the GS-KDS-MSS stations, together with the proximity of the KDS station to the Wenchuan’s epicentre, makes these stations very sensitive in recording radon disturbances with hidden traces of significant precursory value. This is reinforced by the fact that the SPS-PHZS stations are in a completely different fault-line despite in a common fault with the KDS station.The fact (as mentioned above) that KDS is also in line with the other two stations (GS, MSS) and near to the earthquake’s epicentre, makes the recordings of this station even more important. The huge magnitude ($M_w$=7.9) and the low depth (19 $km$) which yielded to a huge energy release, in association with the geologically sensitive background and the proximity of the KDS station to the earthquake ’s epicentre, may explain why post activity is observed only in this station. Indeed, there is a great period from approximately window 3000 (approximately $day$ 133) up to window 6100 (approximately $day$ 271) where the DFA scaling exponents $\alpha$ are quite higher than 1.5 and several are above 1.8 ($0.8\le H<1.0$). These DFA variations are addressed just after to the great Wenchuan earthquake ($day$ 133) and it is the first time that such post activity is found. The reader should note here for integrity that there is also an extended period near sample 2000 (approximately $day$ 88) and a short near sample 8200 (approximately $day$ 363) of completely non-successful ($r^2<0.95$) segments. This latter period provided false DFA exponents above 2 and for this reason they were cut-off from the figure. The small period around 8200 is a Class-II one. Except from these two non-precursory cases, there are also scattered Class-II exponents of low-precursory value. These are non-successful ($r^2<0.95$) fBm segments or fGn segments. Finally, as regards to the pre-seismic signs, from the geological aspect, these are in line with the literature findings since according to the reviews for radon and electromagnetic precursors [2,4,5,6,7,8,9], the epicentre’s distance of precursory activity of KDS-GS-MSS stations, is within the reported ones. Concluding here with the visual observations from the DFA results, the reader should summarise on the following facts: a) DFA outlines hidden trends not (usually) observed in the signal; b) the time-window of the identified activity is sufficient to accept it as precursory; c) the distance and the path way provide some explanations for the identification in 3 out of 5 stations; d) most importantly, even for the robust DFA method, the visual observation is not enough and a combination of different methods is necessary in order to find the precursory time periods with enhanced importance and combined evidence [1,27,29,51]. All these will be presented later in text together with combined evidence. Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11 demonstrate the temporal evolution of the fractal dimensions calculated using Katz’s, Higuchi’s and Sevcik’s methods. There are noticeable differences found. The three methods’ computed values for the fractal dimension show variations as well. All discrepancies are due to the varied calculation techniques used by the three fractal dimension techniques. This has been acknowledged in recent publications [1,15,29]. As pointed out, the methods of Katz and Higuchi estimate higher fractal dimensions than the one of Sevcik. As general trends, the Katz’s fractal dimensions are around 1 and 1.1 ($0.9<H<1$, equation 20), the ones of Higuchi’s method are roughly between 1.5 and 2 ($0<H<0.5$, equation 20)and those of the Sevcik’s method, approximately, between 1 and 1.5 ($0.5<H<1.0$, equation 20). However, what is important is not the value range but the drop-downs of the calculated fractal dimensions especially when these have a noted duration. For that reason, the details of every sub-figure of Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11 has great importance, especially in association with the results of the DFA method.

In reference to the KDS station and in the first 200 windows (from $day$ 1 to $day$ 8) of Figure 7, a noteworthy decrease is observed in Sevcik’s fractal dimension between 1 and 1.4 ($0.6<H<1.0$) and in Higuchi’s D between 1 and 2 ($0<H<1.0$). Katz’s method estimates, unexpectedly, higher fractal dimensions between 1.0 and 1.2($0<H<0.8$) in this area. Simultaneous is an increase in groundwater radon. This is a rare, serendipitous finding that has been reported in other publications [2,8]. Abrupt drops in Higuchi’s fractal dimension are observed, next, between windows 1900 ($day$ 84) and 2200 ($day$ 97). The drops in Sevcik’s fractal dimension are between window 1800 ($day$ 79) and 2000 ($day$ 88). The Sevcik’s fractal dimensions between window 2000 and 2200 were above 2 and for this reason they were cut-off, being considered as error in calculations. Once again, the Katz’s fractal dimensions are, peculiarly, higher around window 1800. As regards the GS station, a significant decrease in Higuchi’s fractal dimension is spotted in Figure 8, between window 40 ($day$ 48) and window 110 ($day$ 133). The decrease in Sevcik’s fractal dimension is roughly synchronous and of the same profile as the one of Higuchi’s, however, milder and with smaller duration, between window 60 ($day$ 73) and 110 ($day$ 133). Fractal dimensions of the MSS station (Figure 9), exhibit also decrease profiles through the calculations with the Higuchi’s and Sevcik’s methods. These decease in D values are observed between windows 70 ($day$ 85) and 140 ($day$ 170). In summary here the key periods from the fractal dimension calculations are for the KDS station between $day$ 1 and $day$ 8 and between $days$ 84 and 97. For the GS station the period is between $day$ 48 and $day$ 110 and for the MSS station between $day$ 85 and $day$ 170. Since the Wenchuan earthquake occurred on $day$ 133, the fractal dimension variations of KDS and GS can be considered, most probably, pre-seismic. The same is valid for the $day$ 85 - $day$ 133 variations of the MSS station. These findings reinforce the claims expressed in the discussion of DFA above since the fractal dimension calculation techniques are completely different from the one of DFA. Moreover, two fractal techniques show the tendencies. and, significantly, in comparable time intervals. These facts support further the necessity of using different fractal techniques in parallel. This has been expressed with emphasis in recent publications [1,15,51,76]. However, as with the DFA outcomes of the KDS station, there is a wide period between windows 3000 ($day$ 133) and 6100 ($day$ 271) during which, the fractal dimensions from Higuchi’s and Sevcik’s methods exhibit very significant variations. These period is the same as the corresponding one discussed for the DFA results and refers to post-seismic variations. Additional post-seismic variations are addressed here via Higuchi’s and Sevcik’s fractal dimension from the MSS station between $day$ 133 and $day$ 170. The reader should note here that as with the outcomes of DFA, SPSS station Figure 11, does not show certain D patterns. Despite that DFA did not provide trends in DFA profiles for the PZHS station, the corresponding Higuchi’s fractal dimension shows a small decrease between window 200 and 250. Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16 show the results from the power law method. This method is one of the wider used techniques and has been considered as one of the most powerful to identify the hidden patterns in time series [33,67,68,69,84,87,90,91,99,101,102,103,104,105,106,107,108]. When, hence, certain trends are found with the power-law method, there is strong evidence on the underlying long-memory of the associated geo-system. At first, as with the outcomes of the DFA and fractal dimension methods, the time evolution of power-law exponent, $\beta$, differs from the one of the time series. In order, however, to discuss the interesting results of each figure, the following information should be taken into consideration for the successive fractal segments:

• If $1.0<\beta \le 3.0$, then the associated time series is a temporal fractal and follows the Class-I category;
  • If $1.0<\beta <2.0$, the time series follows anti-persistent paths;
  • If $2.0<\beta <3.0$, the time series follows persistent paths;
• If $-1.0\le \beta <1.0$, the time series is of low predictability and follows the Class-II category; Moreover:
  • If $\beta =1.0$, the fluctuations of the related processes are not growing and, hence, a stationary system describes the series;
  • If $\beta =2.0$, the underlying dynamics are random and the related system has no memory;

In Figure 12 for the KDS station, as with DFA and the fractal dimension calculation techniques, a first period is observed up to window 2200 (approximately $day$ 120) with scattered successive fractal windows with $1.7<b<2.2$. These fractal epochs correspond to predictable Class-I segments with interchange between persistency and antipersistency. According to publications (e.g., [91]) this is a sign of precursory activity. Interestingly, this epoch is almost identical to those identified with DFA and the three fractal dimension techniques and this is very important. Figure 13 for the GS station shows a first period between window 1 and 40 ($day$ 1 to $day$ 58) and window 50-90 ($day$ 73 to $day$ 131). Both these periods have areas with $b>2$ (Class I, persistent) and can be considered as precursory of the great Wenchuan earthquake. The last period at the end of the analysed windows has $b<1.6$ and, as with DFA and fractal dimension techniques is on low predictability and , hence of low precursory value. This combined finding provides more evidence. In Figure 14 for the MSS station, there are two periods, between window 30 ($day$ 43) and window 60 ($day$ 87), and, between window 70 ($day$ 102) and window 110 ($day$ 160). Although these periods are anti-persistent they are within the middle part of predictable value range of the Class-I category. Since the periods match with those of the other techniques there is probability that they might be signs (pre and post) of the great Wenchuan earthquake. As with DFA and fractal dimension techniques, there is a great post seismic region within the same period as the one identified with DFA and fractal techniques. It is very important that the findings of the techniques match even for the two other stations the PZHS and SPS. Both showed no successive fractal window. This fact reinforces the findings of the DFA and Higuchi’s and Sevcik’s for the MSS station. Evaluating the results presented so-far, there is a period up to $day$ 133 (day of occurrence of the great Wenchuan earthquake) which is systematically identified in all fractal analysis results, that is from all methods, for the KDS, GS and MSS stations, whereas the PSZH and SPS stations do not show any noteworthy such period. These fractal epochs are discussed and considered as pre-seismic signatures of the great Wenchuan earthquake. The recognition of signs in the KDS, GS and MSS stations and the lack of systematic signs of the PSZH and SPS stations is attributed to the different geological paths in association to the theory of asperities and the selectivity effects in earthquake related research. It is also extensively discussed, that finding periods of over- or under- threshold fractal values is significant, but most significant is the synthetic finding of common over- or under- threshold areas with different techniques (meta-analysis, Section 3.8). When, and if, such common locations are discovered, the scientific arguments for the existence of a seismic warnings concealed in the time series increases, making the claims stronger. It is therefore very important to apply meta-analysis to the fractal results of KDS, GS and MSS stations so as to increase the presented evidence. Meta-analysis in the PSZH and SPS station’s results is unnecessary since the corresponding signs are not enough. On the other hand, a so great earthquake as the Wenchuan, is expected to be associated with other earthquakes. Moreover, other earthquakes of the overall area might also explain warnings and there is a possibility that the presumed post-seismic sings might be also be pre-seismic of other earthquakes of the area and this is the first time to mention such a view in this paper. Towards this, the earthquakes of 2008 with $M_w\ge$5.5 were accessed from USGS for an area greater than the inserted one of Figure 1 [109]. The latter reference presents online the selected area together with data of 19 earthquakes and a map representing it in a terrain of China. For integrity the data of these 19 earthquakes are presented in Table 2, whereas Figure 17 presents the earthquakes in a Google Earth Map after creating the corresponding map kml file from USGS and importing it to Google Earth. In this way, Figure 17 presents the earthquakes alternatively in comparison to USGS map [109].

In the consensus expressed above the final step of the related analysis should include: (a) the above 19 earthquakes (Wenchuan included); (b) the fractal results and the discussion on threshold setting presented in Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16; (c) the logic of connecting fractal methods expressed in Section 3.8. Towards this Figure 18 presents combined plots for each station. This multi figure conceals all the important findings. Similar representation has been adopted in other publications as well [1,27,29,65]. The subplots are mixed and have several symbols. For this reason it is crucial to delineate the significant information given:

(a)

All over- or under-threshold results of all fractal methods (step 1, meta-analysis) for the KDS, MSS and GS stations. The threshold results of each station are combined per two, three, four and five methods (step 2, meta-analysis) (total 13 combinations) versus all 19 earthquakes of Table 2 and Figure 17.

As mentioned in the previous paragraph, it is not only to identify footprints by one or more techniques (done already here as well) but more important is to link different techniques focusing on similar aspects of the problem in hand. To achieve this:

(1)

The exact over or under threshold dates were located computationally from the fractal outputs of each station (step-1, meta analysis). These dates are Year, Month, Day, Hour for the HSR station KDS and Year, Month, Day for MSS and GS stations. This is done through a serial search.

(2)

The common threshold dates from all different techniques were found through an incremental computational search. The outputs used are from all methods and with up to 13 different combination of these. All outputs were generated through special software and are stored in computer for use.

(3)

The earthquake data from USGS [109] were transformed to an adequate ASCII file for the generation of the final plot.

(b)

Wherever the symbols of different methods coincide in time, this means that the signs of seismicity is provided by more than one method. If all 13 methods coincide this means that the evidence is maximised. More techniques pointing to similar findings, more rigid the evidence is. The reader should stress that the coincidence is done on the step 1 results, that is on the fractal outputs.

In the above sense and starting from the combined findings of lesser importance, it can be observed from sub figure a of Figure 18 that there are several windows of the combination of DFA versus Higuchi’s and Sevcik’s methods (⊡,plot needs zoom to show the details) which are in the period of all 19 earthquakes. The fact that there are three techniques that show this behaviour, makes these fractal disturbances noteworthy. However, it is not possible to discriminate whether these disturbances are post seismic of the the great Wenchuan earthquake or pre-seismic of another one of the area and this is a limitation of the present methodology. On the other hand, observing carefully the data of Table 2, it is seen that earthquakes 13-18 are practically the seismic sequence of the Wenchuan earthquake since all occurred the same day as earthquake 19. Especially 16-18 happened at the same hour. Moreover earthquakes 11 and 12 are also within the post sequence of the Wenchuan earthquake because it occurred just one (11) and two (12) days after. Moreover all earthquakes from 8-19 occurred on the same month (April 2008). This peculiar facts complicate the issue of discriminating the above common 3-technique concurrent findings from characterising as post- or pre-sesmic. This is not the case for the three co-occurrences of the combinations of the fractal dimension methods ($pentagon$-yellow, namely Sevcik’s versus Katz’s and Higuchi’s methods (3 techniques), $pentagon$-blue, viz. Katz’s versus Higuchi’s and Sevcik’s methods (3 techniques) and green ⊲, namely Higuchi’s versus Katz’s and Sevcik’s methods (3 techniques)). These occurred prior to events 7,6 and 5 but could be prior to 4,3,2,1 as well (decreasing time distance between fractal disturbance and earthquake occurrence) (all these sub figures need zoom as well).

The most significant finding of this paper is left for its last part. It is clearly observed in all sub-figures that the vast majority of coincidences are prior to the great Wenchuan earthquake (19) and its synchronous post-earthquakes (10-19). This is the most important observation and the reader should stress on it. It is not only one method, but all 13 methods which coincide. Moreover during this preparation phase of the great Wenchuan earthquake, there are far more coincidences of even 3 or less techniques as well (such as red ⊳-sub figure b, namely power law analysis versus Katz’s and Sevcik’s methods (3 techniques) and yellow ⊡, that is DFA versus Higuchi’s and Sevcik’s methods (3 techniques)). In the most emphatic manner it is here declared that all these fractal disturbances (importantly, with results from meta-analysis), are, most possibly, due to the great Wenchuan earthquake. This is the most important finding of this paper and is left as last to place emphasis. It is the great magnitude of the Wenchuan, the post-seismic activity of it in association to the novel linking-combining several different fractal techniques that managed to outline such an important finding.