Vol49,4_5,2006 961 ANNALS OF GEOPHYSICS, VOL. 49, N. 4/5, August/October 2006 Key words seismic regime – strong electrical dis- charges – non-linear dynamics 1. Introduction The dynamics of seismic process is far from being clearly understood and modeled. It pres- ents recent several aspects showing that seismic- ity is certainly not a pure random process under a multidisciplinary approach. Magnitude, time and spatial distribution of earthquakes present aspects of self similarity or fractal character evi- denced by several authors (De Rubeis et al., 1993; Turcotte, 1997). On the other hand seis- micity cannot be deterministically explained al- though efforts have been made to show its quasi periodic character. A direct consequence of this situation is the almost complete impossibility to precisely predict earthquakes (Geller et al., 1997; Main, 1999). In recent years, non-linear dynamics has of- fered several tools to analyze and characterize seismicity. Concepts like phase space recon- Influence of strong electromagnetic discharges on the dynamics of earthquakes time distribution in the Bishkek test area (Central Asia) Tamaz Chelidze (1), Valerio De Rubeis (2), Teimuraz Matcharashvili (1) and Patrizia Tosi (2) (1) Institute of Geophysics, Georgian Academy of Sciences, Tbilisi, Georgia (3) Istituto Nazionale di Geofisica e Vulcanologia, Roma, Italy Abstract From 08/01/1983 to 28/03/1990, at the Bishkek ElectroMagnetic (EM) test site (Northern Tien Shan and Chu Val- ley area, Central Asia), strong currents, up to 2.5 kA, were released at a 4.5 km long electrical (grounded) dipole. This area is seismically active and a catalogue with about 14100 events from 1975 to 1996 has been analyzed. The seismic catalogue was divided into three parts: 1975-1983 first part with no EM experiments, 1983-1990 second part during EM experiments and 1988-1996 after experiments part. Qualitative and quantitative time series non- lin- ear analysis was applied to waiting times of earthquakes to the above three sub catalogue periods. The qualitative approach includes visual inspection of reconstructed phase space, Iterated Function Systems (IFS) and Recurrence Quantification Analysis (RQA). The quantitative approach followed correlation integral calculation of reconstruct- ed phase space of waiting time distribution, with noise reduction and surrogate testing methods. Moreover the Lem- pel-Ziv algorithmic complexity measure (LZC) was calculated. General dynamics of earthquakes’ temporal distri- bution around the test area, reveals properties of low dimensional non linearity. Strong EM discharges lead to the increase in extent of regularity in earthquakes temporal distribution. After cessation of EM experiments the earth- quakes’ temporal distribution becomes much more random than before experiments. To avoid non valid conclusions several tests were applied to our data set: differentiation of the time series was applied to check results not affected by non stationarity; the surrogate data approach was followed to reject the hypothesis that dynamics belongs to the colored noise type. Small earthquakes, below completeness threshold, were added to the analysis to check results robustness. Mailing address: Prof. Tamaz Chelidze, Institute of Geo- physics, Georgian Academy of Sciences, Alexidze Street 1, 380093 Tbilisi, Georgia; e-mail: chelidze@ig.acnet.ge 962 Tamaz Chelidze, Valerio De Rubeis, Teimuraz Matcharashvili and Patrizia Tosi struction and fractal dimension of the attractor help to distinguish between a purely random process and a complicated process driven by a finite, limited set of rules. This was allowed by the recognition that non linearity can produce complex dynamic behavior in systems driven by a finite number of factors. The enormous gap between «simple» lin- ear deterministic models and random, compli- cated and strongly unpredictable processes seems to be filled with these new analytical tools. The aim is to render tractable, in a cer- tain way, phenomena and data otherwise not clearly depicted. The present work investigated the influence of strong EM discharges on earthquakes tempo- ral distribution. Experiments on the triggering effect of MHD (magneto-hydro-dynamic) soundings on the microseismic activity of the region were performed in 1975-1996 by IVTAN (Institute of High Temperatures of Russian Academy of Sciences) in the Central Asia test area (Tarasov, 1997; Tarasov et al., 1999; Jones, 2001). Dur- ing these experiments deep electrical sounding of the crust was carried out at the Bishkek test site in 1983 to 1989. The source of electrical energy was the MHD generator, and a load was an electrical dipole of 0.4 Ω resistance with electrodes located at a distance of 4.5 km from each other. When the generator was fired, the load current was 0.28-2.8 kA, the sounding pulses had durations of 1.7 to 12.1 s, and the en- ergy generated was mostly in the range of 1.2- 23.1 MJ (Volykhin, 1993). Evidence of some relationships between EM discharges and seismic activity have been pointed out under a statistical aspect and in a time range of days after EM experiments (Tarasov, 1997). Here the general dynamical as- pect is considered. A good seismic catalogue of the area is available well before, during and well after this period. A simple causal relation- ship between the two processes is not strongly evident. Relations appear to be present but data noise is also relevant. It is useful to investigate if the seismic dynamics, in periods before, dur- ing and after EM experiments, is influenced by the introduction of strong electric current into the ground. 2. Methods The investigation was performed according to the general scheme of time series non-linear analysis (Abarbanel et al., 1993; Sprott and Row- lands, 1995; Kantz and Schreiber, 1997; Goltz, 1998; Hegger and Kantz, 1999). In general data analysis can be performed firstly under a more qualitative and visual approach and subsequently a more quantitative methodology can be applied. The qualitative approach includes visual in- spection of reconstructed phase space. Namely, p-dimensional phase space from the scalar time sequences was reconstructed by the method of time delay (Packard et al., 1980; Takens, 1981). According to Takens’ theorem it is possible to catch the essential dynamic properties of a sys- tem by a reconstruction of its phase space by only one variable. Two- and three-dimensional phase space portraits encapsulating essential dynamic properties of the analyzed complex process were used as qualitative tests. Also oth- er qualitative tools have been used such as Iter- ated Function Systems (IFS) (Jeffrey, 1992) and Recurrence Plots (RP) (Eckman et al., 1987). Generally recurrence analysis is a graphical method designed to locate hidden recurring pat- terns and structure in time series. The RP is de- fined as (2.1) here εi is a cut-off distance (often ε = 0.1σ, σ- standard deviation), Θ(x) is the Heaviside func- tion, Θ(x) = 0 if x < 0 and Θ(x) = 1 if x ≥ 0. Ac- cording to Eckman the values one and zero in this matrix commonly are visualized as black and white (Eckman et al., 1987). The black points indicate the recurrences of the investigat- ed dynamic system revealing their hidden regu- lar and clustering properties. By the definition, RP has a black main diagonal (line of identity) formed by distances in the matrix compared with themselves. In order to understand RP it should be stressed that it visualizes the distance matrix which represents autocorrelation in the series at all possible time (distance) scales. As far as distances are computed for all possible pairs, on the RP plots elements near the diagonal correspond to short range correlation, whereas R x x,i j i i jεΘ= − −_ i Influence of strong electromagnetic discharges on the dynamics of earthquakes time distribution ... the long range correlations are revealed by the points distant from the diagonal. Hence if the an- alyzed dynamics (time series) is deterministic (ordered, regular), then the recurrence plot shows short line segments parallel to the main diagonal. Qualitative patterns of unknown dynamics presented as fine structure of RP are often are too difficult to be considered in detail. Therefore one uses a modern quantitative method of analysis of complex dynamics for RP approach (Recurrence Quantitative Analysis or RQA) (Webber and Zbilut, 1992). RQA is especially useful to over- come the difficulties often encountered dealing with non-stationary and rather short real data sets. As a quantitative tool of complex dynamics analysis RQA defines several measures mostly based on diagonally oriented lines in the recur- rence plots: recurrence rate, determinism, maxi- mal length of diagonal structures, entropy, trend, etc. In the present work recurrence rate − RR(t) and determinism − DET(t) measures based on analysis of diagonal oriented lines in the recur- rence plot have been calculated (Weber and Zbi- lut, 1994; Marwan et al., 2002). Generally the Recurrence Rate RR(t) is the ratio of all recurrent states (recurrence points) to all possible states and is therefore the proba- bility of the recurrence of a certain state. Sto- chastic behavior causes very short diagonals, whereas deterministic behavior causes longer diagonals. The ratio of recurrence points forming diag- onal structures to all recurrence points is called the determinism DET(t). DET(t) is the propor- tion of recurrence points forming long diagonal structures of all recurrence points. Again, sto- chastic and heavily fluctuating data cause none or only short diagonals, whereas deterministic systems cause longer diagonals. An Iterated Function System (IFS) is a Hutchinson operator on Kk (finite set of functions moving points around in some space), which maps a set of points to another set of points. If a Hutchinson operator is repeatedly applied to a compact set, in the limit it will render the unique attractor of the IFS (Peitgen et al., 1992). For time series analysis purpose IFS attractors can be used as qualitative measures of self similarity of analyzed dynamics (more order is in time series more regular are structures in IFS attractor). We use IFS as an additional qualitative tool for de- tection of hidden structure in analyzed time se- ries (Sprott and Rowlands, 1995). These tests enables a first qualitative visual inspection of unknown dynamics and help to un- cover general properties of analyzed process. Qualitative analysis reveals the possible existence of specific attractors, e.g., strange ones which point to the deterministic chaotic behavior. As the main tool for quantitative analysis of earthquakes dynamics, correlation integral cal- culation of reconstructed phase space of tempo- ral distribution has been performed (Abarbanel et al., 1993; Kantz and Schreiber, 1997; Hegger and Kantz, 1999). This approach is based on idea of correlation sum. Correlation sum C(r) of a set of points in the vector space is defined as the fraction of all possible pairs of points which are closer to each other than a given distance r. The basic formula useful for practical application is (2.2) where again Θ(x) is the Heaviside step func- tion, ; x i − x j; is the Euclidian norm, i = j are ex- cluded. For fractal systems for enough long time series and for a small r, C(r)∝r ν relation- ship is correct. Commonly such dependence is correct only for the restricted range of r values, so called scaling region. Correlation dimension ν or d2 is defined as . (2.3) In practice d2 value is found from the slopes of logC(r)/log(r) curves for different phase space dimensions. The correlation dimension of an unknown process is the saturation value of d2, which does not change by increasing of phase space dimension. In order to reduce possible spurious conclu- sions on considered dynamics, noise reduction and surrogate testing methodologies were used (Kantz and Schreiber, 1997; Hegger and Kantz, 1999). Besides, as additional quantitative test for rel- atively short time series, the Lempel-Ziv algorith- mic complexity measure (LZC) was calculated ( ) ( ) lim og log l d r C r r 2 0 ν = = " ( ) ( ) C r N N r x x 1 2 i j j i N i N 11 Θ= − − − = += _ i// 963 964 Tamaz Chelidze, Valerio De Rubeis, Teimuraz Matcharashvili and Patrizia Tosi (Lempel and Ziv, 1976). It is based on the trans- formation of the original one dimensional time series into a finite symbol sequence and is de- fined as (2.4) where N is the length of the original time series, L (N) ∼ Nw (N)(logbNW (N )+1) is the total length of the encoded sequence, where Nw (N ) ≤ N is the total number of code words. Being one of the tools for time series non-linear analysis, LZC is ( ) lim supC N L n LZ N = "3 especially suitable for relatively short real data sets because it is not so demanding to the time se- ries length as other methods (Zhang and Thakor, 1999; Matcharashvili and Janiashvili 2001). 3. Data and analysis In the present study non-linear analysis was performed on about 14100 time intervals (in seconds) between earthquakes from 1975 to 1996. In the original catalogue, the energy of Fig. 1a,b. Results of the completeness analysis of the Ivtan Site seismic catalogue. a) Cumulative number of events versus time for magnitude class step = 0.5 for 0.< m > 3.0. Note that cumulative number of events is rescaled among magnitude classes. b) Log cumulative number of earthquake versus magnitude (Gutenberg- Richter relation) for the whole catalogue; values of regression fit equation: Y = − 0.83∗X + 5.40. Coefficient of de- termination, R-squared = 0.995. a b Influence of strong electromagnetic discharges on the dynamics of earthquakes time distribution ... the events was expressed as energy class, which we converted to magnitude using the following relation (Riznichenko, 1985) (3.1) where m is magnitude and E is the energy class. The completeness of the catalogue was in- vestigated first by considering the realization of the Gutenberg-Richter relationship at low mag- nitudes: departure from a straight line was in- terpreted as a lack of completeness at low mag- nitudes. As a result the catalogue was consid- ered complete, under the sole magnitude aspect, for m ≥ 1.7. The Gutenberg-Richter b-value was found to be equal to 0.83 with a reasonably good fit. Earthquakes with magnitude higher than 6 seem to show behavior typical of characteristic events. A second test was oriented to check the time completeness. As is well known, a catalogue’s completeness changes with time, usually as a re- sult of improving seismic-network performance (e.g., increasing number of stations), leading to greater magnitude sensitivity. The completeness analysis was performed by employing the method of Mulargia et al. (1987). The method consists in separating all events into magnitude classes and plotting separately the cumulative number of events versus time. Assuming that during the considered time interval seismicity had a constant rate, the flat behavior at the begin- ning of the time period may be due to a lack of data; this is normal for low-magnitude ranges. Only for magnitudes higher than 2.0 is our catalogue complete over the entire time period (number of earthquakes n = 5297). If a lower magnitude limit is desired, the time period from year 1980 is more appropriate (fig. 1a,b). As a result of the analysis performed, a relatively complete catalogue was obtained with a lower magnitude threshold of 1.7 from the year 1980. For the present study the catalogue was an- alyzed under the time aspect, specifically on in- ter-event (waiting) times. A catalogue subset of waiting times was used according to complete- ness analysis, i.e. all time spans and m > 2.0. Then all data used were selected by the same procedure for confirmation of results and to test their robustness. m E 18 4 = − 4. Results and discussion Figures 2a-f presents the results of qualitative analysis of waiting times sequences above the mentioned threshold. Namely in fig. 2a,c,e, IFS clumpiness test (Jeffrey, 1992; Sprott and Row- lands 1995) and in fig. 2b,d,f, recurrence plot analysis (Zbilut and Webber, 1992) results reveal that after the beginning of the experiments some structure in plots is visible, pointing to the in- creased amount of functional interdependence in earthquake temporal distribution. As to the quantitative approach, the variation of correlation dimension versus dimension of phase space where reconstructed dynamics is embedded (embedding dimension) is presented in fig. 3. The integral time series (5297 time in- tervals) for the whole period of observation (1975-1996) containing time intervals sequences between all events above the threshold reveals clear low correlation dimension (d2 = 1.22 ± 0.43, asterisks). Shorter time series also were consid- ered. Namely 1760 waiting times data before (1975-1983), 1953 waiting times during MHD experiments (1983-1988) and 1584 waiting times of the period after experiments (1988-1992). Time series before and especially during MHD experiments also have low correlation dimension (d2 < 5). Namely d2=3.83± 0.80 before and d2 = = 1.04 ± 0.35 during experiments. On the other hand, opposite to what was mentioned above, af- ter cessation of experiments (fig. 3, triangles) cor- relation dimension of waiting times sequences noticeable increases (d2 > 5.0), exceeding low di- mensional threshold (d2=5.0). That means that af- ter termination of experiments the extent of regu- larity or extent of determinism in process of earth- quake temporal distribution decreases. The con- sidered process becomes much more random both qualitatively (fig. 2e,f) and quantitatively (fig. 3. triangles). For clarity in fig. 3, results for random number sequence are shown too (diamonds). The found low correlation dimension of waiting interval time series is in good accor- dance with earlier published results for Cauca- sus (Matcharashvili et al., 2000) as well as with Goltz’s results (Goltz, 1998) for other seismic active regions. This result together with qualitative analysis results shown in fig. 2a-f provide evidence that 965 966 Tamaz Chelidze, Valerio De Rubeis, Teimuraz Matcharashvili and Patrizia Tosi Fig. 2a-f. Qualitative analysis of temporal distribution of earthquakes (complete catalogue, M ≥ 1.7) before the be- ginning of EM experiments (1975-1983), during experiments (1983-1988) and after accomplishing of experiments (1988-1992). IFS-clumpiness test for inter-event time interval sequences: a) before experiments; c) during experi- ments; e) after experiments. Recurrence plots analysis of waiting times sequences: b) before experiments; d) during experiments; f) after experiments. Note diagonal lines in IFS plot and compact structure in RP during experiments. after the beginning of EM discharges the tempo- ral distribution of earthquakes around IVTAN test area becomes more regular, or events of cor- responding time series become functionally much more interdependent. At the same time the absence of typical phase space structure (not shown here), IFS and recurrence plot attractors (fig. 2a-f) do not al- low us to consider the process deterministical- ly chaotic. a c e b d f Influence of strong electromagnetic discharges on the dynamics of earthquakes time distribution ... In order to the reduce effects of possible noise we analyzed waiting time series after a noise reduction procedure (Shreiber, 1993; Kantz and Schreiber, 1997). Namely, we used the meth- odology of non-linear noise reduction (which in fact is phase space non-linear filtering) instead of common linear filtering procedures. The latter, as is well known, may lead to destroying the origi- nal non-linear structure of analyzed complex processes (Hegger and Kantz 1999; Schreiber and Schmitz, 2000). Non-linear noise reduction relies on the exploration of a reconstructed phase space of considered dynamic process instead of frequency information of linear filters (Schreiber, 1993; Kantz and Schreiber, 1997; Hegger and Kantz, 1999). Correlation dimension versus embedding space dimension of noise-reduced time series is presented in fig. 4. As follows from the obtained results, correlation dimensions are not essential- ly affected by unavoidable noises. Therefore ob- tained results assure that the differences found in d2-phase space dimension (P) dependence be- fore, during and after experiments (fig. 3) are in- deed related to dynamic changes in temporal distribution of earthquakes after beginning of MHD discharges experiments. When describing unknown dynamics of wait- ing times fluctuation, the differentiation of origi- nal time series can be useful to avoid improper conclusions related to the effects of trends or non- stationarity in data sets, even when those are not clearly visible (as in the case of considered time series) (Goltz, 1998). As shown in fig. 5, differen- tiation of our time series according to Goltz (1998) does not lead to significant changes in ob- tained results (see fig. 3). So our results could not be affected by trends or non-stationarity in used data sets. 967 Fig. 3. Correlation dimension versus embedding di- mension of waiting times sequences (complete cata- logue): asterisks – integral time series (1975-1996); circles – before beginning of experiment (1975-1983); squares – during experiments (1983-1988); triangles – after experiments (1988-1992); diamonds correspond to random number sequence. Fig. 4. Correlation dimension versus embedding dimension of waiting times sequences (complete cat- alogue) after noise reduction: diamonds – before ex- periments; squares – during experiments; triangles – after experiments. Fig. 5. Correlation dimension versus embedding dimension of differenced waiting times sequences (complete catalogue): diamonds – before experi- ments; squares – during experiments (1983-1988); triangles – after experiments. 968968 Analysis of different time series may be im- portant also in the sense of inherent dynamic structure testing (Prichard et al., 1994). Name- ly, the test is based on finding that estimated non-linear measure (correlation dimension in our case) for the differentiated series is larger than that estimated for original data, if the structure of their dynamics is caused by the lin- ear stochasticity. At the same time for chaotic (low dimensional) processes these measures are the same. From this point of view analysis of differentiated time series before detailed surro- gate testing provides the first additional evi- dence that variation of waiting times indeed has inherent non-linear structure and that their dy- namic properties are not caused by linear rela- tionships between data. Indeed, the curves of figs. 3 and 5 are characterized by comparable values of correlation dimension. Moreover, in order to have a basis for more reasonable rejection of spurious conclusions caused by possible linear correlations in consid- ered data sets, we have used the surrogate data approach to test the null hypothesis that our time series are generated by a linear stochastic process (Theiler et al., 1992; Rapp et al., 1993, 1994; Kantz and Schreiber, 1997). In other words we would like to reject correctly the pos- sibility that revealed dynamics belongs to the colored noise type. Namely, Random Phase (RP) and Gaussian Scaled Random Phase (GSRP) surrogates sets for waiting times series were used (Matcharashvili et al., 2000). RP- surrogate sets are obtained by destroying the non-linear structure through randomization of phases of Fourier transform of original time se- ries and then performing a backward transfor- mation. GSRP surrogate sets were generated in a three-step procedure. At first a Gaussian set of random numbers was generated, which has the same rank structure as the original time series. After this phase randomised surrogates of these Gaussian sets were constructed. Finally the rank structure of the original time series was re- ordered according to the rank structure of the phase randomized Gaussian set (Theiler, 1992). Figure 6a,b shows the results of surrogate testing of waiting time sequences before and dur- ing experiments using d2 as a discriminating met- ric. For each of our data sequences, we generated 75 of RP and GSRP surrogates. There are sever- al ways to measure the difference between the discriminating metric measure of original (given by Morig) and surrogate (given by Msurr), time se- ries (Rapp, 1994). Investigators often use the significance criterion S=⎪〈 Msurr〉 − M orig⎪/σsurr, where σsurr is standard deviation of Msurr (Theil- er et al., 1992). The significance criterion S, according to Theiler et al. (1992) for analyzed time series before experiments is significant: 22.4 ± 0.2 for RP and 5.1 ± 0.7 for GSRP surrogates. Conse- quently after the beginning of experiments the null hypothesis that original time series is the linearly correlated noise was rejected with sig- nificant value of S criterion: 39.7 ± 0.8 for RP and 6.0 ± 0.5 for GSRP surrogates. Fig. 6a,b. Correlation dimension versus embedding dimension of original (diamonds) and surrogate (squares – GSRP, triangles – RP) waiting time se- quences: a) before beginning of experiments, b) dur- ing experiments. a b Tamaz Chelidze, Valerio De Rubeis, Teimuraz Matcharashvili and Patrizia Tosi 969 Influence of strong electromagnetic discharges on the dynamics of earthquakes time distribution ... 969 These results can be considered strong enough evidence that analyzed time series are not linear stochastic noise and that correspon- ding processes of earthquakes temporal distri- bution before and especially during experi- ments are characterized by inherent low dimen- sional non-linear structure. According to the IVTAN catalogue, each considered time period contains one large (M ≈ 6.1-6.3) earthquake. Therefore in order to refine whether the above results are caused by special properties of a separate large earthquake or reflect total changes in dynamics caused by EM discharges, we analyzed waiting time se- Fig. 7a-f. Qualitative analysis of 1000 data waiting times sequences (complete catalogue), after largest events before beginning of EM experiments (1975-1983), during experiments (1983-1988) and after accomplishing of experiments (1988-1992). IFS-clumpiness test for inter-event time interval sequences: a) before experiments; c) during experiments; e) after experiments. Recurrence plots analysis of waiting times sequences: b) before exper- iments; d) during experiments; f) after experiments. a b d f c e 970 quences (above appropriate threshold) after each largest event. Namely, 1000 consecutive waiting time intervals after 24/03/1978 M=6.1 (K = 15.0), 24/01/1987 M=6.3 (K = 15.3) and 798 time intervals after 30/12/1993 M=6.1 (K = 15.0) events were analyzed. It is important to say, that each of these relatively short time series are localized in corresponding time peri- ods named above as «before», «during» and «after» experiments. It becomes clear from IFS-clumpiness and RQA analysis results (fig. 7a-f) that qualitative- ly this situation is like that shown in fig. 2a-f, i.e. after the beginning of experiments, consid- ered dynamics become more regular and, after accomplishing experiments, dynamics is most random-like. Quantitatively, it is shown in fig. 8 that these short time series generally reveal that after ex- periments analyzed dynamics also become more random than before. Some differences are noticeable in saturation values of correlation di- mension (in fig. 8) before (circles, d2 = 3.1 ± 0.4) and during (squares, d2 = 2.1 ± 0.7) experiments. The latter may be caused by too limited data length for proper non-linear analysis of these time series (untypical shape of curve at high Fig. 8. Correlation dimension versus embedding dimension of 1000 data waiting times sequences (complete catalogue) after largest events: circles – time period before beginning of experiments (1975- 1983); squares – time period during experiments (1983-1988); triangles – time period after accom- plishing of experiments (1988-1992). Tamaz Chelidze, Valerio De Rubeis, Teimuraz Matcharashvili and Patrizia Tosi embedding dimensions) as well as by artificial- ly increased fraction of aftershocks in short time series which contains events only after the largest earthquakes. In any case, our main conclusion about low dimensional dynamical structure of earthquake temporal distribution during experiments and increasing randomness after termination re- mains valid even for periods of separate large earthquakes. The above conclusion on the increase in reg- ularity in earthquakes’ temporal distribution after the beginning of experiments in some degree is confirmed also by the results of Lempel-Ziv algo- rithmic complexity (CLZ) measure calculation (Lempel and Ziv, 1976). Indeed CLZ is larger when necessary code words are longer i.e. when regular patterns of analyzed time series are minor. Indeed, measured values of Lempel-Ziv complexity before, during and after experi- ments for original time series above threshold consequently are CLZ = 0.99 ± 0.07; CLZ = 0.87± ± 0.05; CLZ = 1.00 ± 0.08. Quantitative RQA results also lead to the same conclusion: namely RR(t) = 9.6, DET(t) = = 3.9 before experiments, RR(t) = 25, DET(t) = = 18 during and RR(t) = 3, DET(t) = 1.5 after ex- periments. The increase in order in earthquake temporal distribution under EM influence is confirmed for short time interval sequences above threshold af- ter the largest earthquakes. Indeed Lempel-Ziv complexity measure values were: CLZ = 0.98 ± ± 0.08; CLZ = 0.74 ± 0.05; CLZ = 1.00 ± 0.09 be- fore, during and after MHD runs consequently (note that CLZ = 0.04 for periodic and CLZ = 1 for random processes). Also, the increase in order in temporal distribution is documented by RQA results for mentioned short time series; namely RR(t) = 9.8, DET(t) = 6.5 before experiments, RR(t) = 19.5, DET(t) = 49.3 during and RR(t) = = 7.1, DET(t) = 1.6 after experiments. In other words for situations where the shape of d2 (fig. 8), is not informative for find- ing changes in dynamics possibly due to too short time series, Lempel-Ziv and RQA analy- sis reveals the increase in regularity. On the basis of results of previous research it is known that, small earthquakes play a very important role in the general dynamics of earth- 971 Influence of strong electromagnetic discharges on the dynamics of earthquakes time distribution ... quake temporal distribution (Matcharashvili et al., 2000). Therefore, additionally we carried out an analysis of time series containing all available from the whole catalogue waiting time sequences, including those between small earthquakes below magnitude threshold. This test is also valid to check results robustness in case of adding a new not necessarily complete set of data to our original set. The total number of events in the whole catalogue increased up to Fig. 9a-f. Qualitative analysis of temporal distribution of earthquakes including small events (whole catalogue, all events) before beginning of EM experiments (1975-1983), during experiments (1983-1988) and after accom- plishing of experiments (1988-1992). IFS-clumpiness test for waiting times sequences: a) before experiments; c) during experiments; e) after experiments. Recurrence plots analysis of inter-event time interval sequences: b) before experiments; d) during experiments; f) after experiments. a c e b d f 972 Tamaz Chelidze, Valerio De Rubeis, Teimuraz Matcharashvili and Patrizia Tosi 14100, while in the complete catalogue for the three above-mentioned periods (before, during and after MHD experiments) there were about 4000 data in each one. Results of IFS and recurrence plots tests analysis of these time series are shown in fig. 9a-f. Noticeable qualitative differences in wait- ing time intervals distribution dynamics during as well as after accomplishment of experiments is obvious. Figure 10 presents the results of correlation dimension calculation for these time series. Prac- tically there are no differences from results ob- tained for the case with m > 2.0 (fig. 3). Namely, according to fig. 10, the integral time series (14100 time intervals) for whole period of ob- servation (1975-1996) reveals clear low corre- lation dimension (d2 = 2.40 ± 0.71) (diamonds). The time series before beginning of experiment (squares) is characterized by correlation dimen- sion (d2 = 3.50 ± 0.63) which is still below ac- cepted low dimensional threshold (d2 = 5.0). During experiments (fig. 10, triangles) correla- tion dimension of time interval sequence, in comparison with the situation before, notice- ably decreases (d2 = 1.71 ± 0.09). After termina- tion of experiments the correlation dimension of waiting time interval sequences noticeably increases (d2 > 5.0), exceeding low dimensional threshold (d2 = 5.0). As in the case of the com- plete catalogue, it means that after termination of experiments the extent of determinism in the process of earthquake temporal distribution de- creases. The considered process becomes much more random both qualitatively (fig. 9c,f) and quantitatively (fig. 10 circles). Both the complete and whole catalogues of waiting time sequences reveal low dimensional non-linear structure in temporal distribution of earthquakes before and especially during exper- iments, which was confirmed by 70 surrogate testing analyses (fig. 11a,b). The significance criterion S for analyzed time series before ex- periments gives: 32.3 ± 0.2 for RP and 5.3 ± 0.6 for GSRP surrogates; consequently after the be- ginning of experiments the null hypothesis that Fig. 10. Correlation dimension versus embedding dimension of waiting times sequences of the whole catalogue: diamonds – integral time series (1975- 1996); squares – before beginning of experiment (1975-1983); triangles – during experiments (1983- 1988); circles – after experiments (1988-1992). Fig. 11a,b. Correlation dimension versus embed- ding dimension of original waiting time sequences of whole catalogue (triangles) and their surrogate (cir- cles-GSRP, squares-RP): a) before beginning of ex- periments; b) during experiments. a b 973 Influence of strong electromagnetic discharges on the dynamics of earthquakes time distribution ... the original time series is the linearly correlated noise was rejected with significant value of S criterion: 46.2 ± 0.5 for RP and 6.5 ± 0.7 for GSRP surrogates. The correlation dimension versus embed- ding space dimension of noise-reduced time se- ries of the whole catalogue is presented in fig. 12. It is clear from this picture that calculated values of correlation dimension are not affected by noises as for the complete catalogue. Ob- tained results show that differences found in the d2-phase space dimension (P) relationship be- fore and during experiments in both catalogues are indeed caused by dynamic changes in tem- poral distribution of earthquakes during EM ex- periments. We also analyzed waiting time sequences after each of the largest (M ≈ 6.1-6.3) events for the whole catalogue, namely, 1000 consecutive waiting time sequences after 24/03/1978 M = = 6.1 (K=15.0), 24/01/1987 M = 6.3 (K=15.3) and 30/12/1993 M = 6.1 (K=15.0) event. As is shown in fig. 13, these short time series gener- ally reveal dynamic characteristics similar to those of the time series obtained from the com- plete catalogue. Found differences which are noticeable in saturation values of correlation di- mension before (circles, d2 = 2.0 ± 1.1 in fig. 13) and during (squares, d2 = 3.2 ± 0.8, fig. 13) ex- periments may be caused both by shortness of these time series or by influence of increased fraction of aftershocks. Thus, conclusions concerning the influence of hot and cold EM runs on the general charac- teristics of earthquakes temporal distribution dy- namics remain valid for small earthquakes too. It is interesting to note that on the laborato- ry scale the effect of triggering and synchro- nization of acoustic emission during slip im- posed by strong EM field is well documented in numerous experiments (Chelidze et al., 2002, 2005; Chelidze and Lursmanashvili, 2003). 5. Conclusions The question whether electromagnetic ex- periments on a specific site can influence the dy- namics of a seismic region is a complex argu- ment. A complete answer to it, if any could be given, would involve a repeated set of analyses for different seismic regions over large periods of time with and without EM experiments. A theoretical explanation showing the cause and effect relationships between the two phenomena is also fundamental. This paper addressed the question under statistical aspect involving non- Fig. 12. Correlation dimension versus embedding dimension of inter-event time interval sequences of whole catalogue after noise reduction: diamonds – before experiments; squares – during experiments; triangles – after experiments. Fig. 13. Correlation dimension versus embedding dimension of 1000 data waiting times sequences of the whole catalogue after largest events: circles – time period before beginning of experiments (1975- 1983); squares – time period during experiments (1983-1988); triangles – time period after accom- plishing of experiments (1988-1992). 974 Tamaz Chelidze, Valerio De Rubeis, Teimuraz Matcharashvili and Patrizia Tosi linear dynamics methods. These methods were chosen because there are not trivial, simple and direct relations between the two phenomena: this means that relations are of a complicated or- der. Moreover seismicity is very probably a crit- ical process with a per se complicated evolution: under given conditions possible relations must not be direct and simple. With non-linear meth- ods the time evolution of seismicity was investi- gated looking at relations with EM experiments. Waiting times constitute the aspect analyzed. The whole time period was divided into three parts, the middle being the one when EM exper- iments took place. Phase space attractor, reconstructed with de- lay time technique, for the whole time period shows low correlation dimension values; this in- dicates, at least, the presence of a few processes driving seismicity. The same analysis on the three sub catalogues confirms the result, with the exception for the period after the EM experi- ments: strong EM discharges lead to the increase in the extent of regularity in earthquakes’ tempo- ral distribution while, after cessation of EM in- fluence, earthquakes’ temporal distribution be- comes much more random than before experi- ments. This is the main result of the analysis and it was confirmed changing the conditions of the analysis itself. Non-linear noise reduced time se- ries confirmed such results as did surrogate test- ing. The middle period contains a large seismic event (January 24, 1987 m = 6.3 derived from energy class K=15.3) this event has certainly a well identified aftershocks activity and this can be a strong factor influencing time dynamics. The root question is: is this event with its related sequence responsible for the change in the dy- namics of analyzed data? If the answer is yes we are forced to answer immediately the new ques- tion if this earthquake is related to EM experi- ments? But it must be noted that inside the other two periods there are also important events of comparable magnitudes and the analysis was conducted on the three sequences of catalogue after each strong event separately. General re- sults confirmation was shown. The same results were revealed with the use of the whole cata- logue regardless of the completeness criteria. This analysis is not certainly exhaustive: the seismic catalogue covers a broad area and all complete data were used, no distinction was made for space location of seismic events. The energy aspect was not fully considered: all events were considered equal regardless of their magnitude. These are strong simplifications and results must be considered under these con- straints. However results appear to be consistent: EM experiments influence seismic time dynam- ics to some extent, increasing the regularity of waiting times. After EM experiments seismic waiting times increase their random character to a level higher than before experiments. REFERENCES ABARBANEL, H.D., R. BROWN, J.J. SIDOROVICH and L.SH. TSIMRING (1993): The analysis of observed chaotic data in physical systems, Rev. Mod. Phys., 65 (4), 1331- 1392. CHELIDZE, T. and O. LURSMANASHVILI (2003): Electromag- netic and mechanical control of slip: laboratory exper- iments with slider system, Non-linear Processes Geo- phys., 20, 1-8. CHELIDZE, T., N. VARAMASHVILI, M. DEVIDZE, Z. CHELIDZE, V. CHIKLADZE and T. MATCHARASHVILI (2002): Labora- tory study of electromagnetic initiation of slip, Ann. Geophysics, 45 (5), 587-598. CHELIDZE, T., T. MATCHARASHVILI, J. GOGIASHVILI, O. LURSMANASHVILI and M. DEVIDZE (2005): Phase syn- chronization of slip in laboratory slider system, Non- linear Processes Geophys., 12, 1-8. DE RUBEIS, V., P. DIMITRIU, E. PAPADIMITRIOU and P. TOSI (1993): Recurrent patterns in the spatial behaviour of Italian seismicity revealed by the fractal approach, Geophys. Res. Lett., 20, 1911-1914. ECKMANN, J.P., S. KAMPHORST and D. RUELLE (1987): Re- currence plots of dynamical systems, Europhys. Lett., 4 (9), 973-977. GELLER, R.J., D.D. JACKSON, Y.Y. KAGAN and F. MULARGIA (1997): Earthquakes cannot be predicted, Sciences, 275, 1616-1617. GOLTZ, C. (1998): Fractal and Chaotic Properties of Earth- quakes (Springer, Berlin). HEGGER, R. and H. KANTZ (1999): Practical implementation of non-linear time series methods: the TISEAN pack- age, Chaos, 9, 413-440. JEFFREY, J.H. (1992): Chaos game vizualization of se- quences, Comput. Graphics, 16 (1), 25-33. JONES, N. (2001): The quake machine, New Scientist, June 30, 34-37. KANTZ, H. and T. SCHREIBER (1997): Non-linear Time Se- ries Analysis (Cambridge University Press). LEMPEL, A. and J. ZIV (1976): On the complexity of finite sequences, IEEE Trans. Infor. Theory, IT-22, 75-81. MAIN, I. et al. (1999): Is the reliable prediction of individ- ual earthquakes a realistic scientific goal?, in Nature De- bates (available on line at http://www.nature.com/nature/ debates/earthquake/equake_frameset.html). MARWAN, N., N. WESSEL, U. MEYERFELDT, A. SCHIRDEWAN 975 Influence of strong electromagnetic discharges on the dynamics of earthquakes time distribution ... and J. KURTHS (2002): Recurrence-plot-based meas- ures of complexity and their application to heart rate variability data, Phys. Rev. E, 66, 026702.1-026702.8. MATCHARASHVILI, T. and M. JANIASHVILI (2001): Investiga- tion of variability of indexes of myocardial contractili- ty by complexity measure in patients with arterial hy- pertension, in Non-linear Dynamics in Life and Social Sciences, edited by W. SULIS and I. TROFIMOVA (IOS Press, Amsterdam), 204-214. MATCHARASHVILI, T., T. CHELIDZE and Z. JAVAKHISHVILI (2000): Non-linear analysis of magnitude and waiting time interval sequences for earthquakes of Caucasian region, Non-linear Processes Geophys., 7, 9-19. MULARGIA, F., P. GASPERINI and S. TINTI (1987): Contour mapping of Italian seismicity, Tectonophysics, 142, 203-216. PACKARD, N.H., J.P. CRUTCHFIELD, J.D. FARMER and R.S. SHAW (1980): Geometry from a time series, Phys. Rev. Lett., 45, 712-716. PEITGEN, H.O., H. JURGENS and D. SAUPE (1992): Chaos and Fractals: New Frontiers of Science (Springer, NY). PRICHARD, D. and J. THEILER (1994): Generating surrogate data time series with several simultaneously measured variables, Phys. Rev. Lett., 73 (7), 951-1018. RAPP, P.E., A.M. ALBANO, T.I. SCHMAH and L.A. FARWELL (1993): Filtered noise can mimic low-dimensional chaotic attractors, Phys. Rev. E, 47 (4), 2289-2297. RAPP, P.E., A.M. ALBANO, I.D. ZIMMERMAN and M.A. JUMENEZ-MONTERO (1994): Phase-randomized surro- gates can produce spurious identification of non-ran- dom structure, Phys. Lett. A, 192 (1), 27-33. RIZNICHENKO, YU.V. (1985): Problems of Seismology (Nau- ka, Moscow), p. 24 (in Russian). SCHREIBER, T. (1993): Extremely simple non-linear noise- reduction method, Phys. Rev. E, 47 (4), 2401-2404. SCHREIBER, T. and A. SCHMITZ (2000): Surrogate time se- ries, Physica D, 142, 346-352. SPROTT, J.C. and G. ROWLANDS (1995): Chaos Data Analyz- er; the Professional Version (AIP, NY). TAKENS, F. (1981): Detecting strange attractors in turbu- lence, in Dynamical Systems and Turbulence, edited by D.A. RAND and L.S. YOUNG, Springer Lecture Notes in Mathematics, 898, 366-381. TARASOV, N.T. (1997): Crustal seismicity variation under electric action, Trans. Russ. Acad. Sci., 353A (3), 445- 448. TARASOV, N.G., N.V. TARASOVA, A.A. AVAGIMOV and V.A. ZEIGARNIK (1999): The effect of high-power electro- magnetic pulses on the seismicity of the Central Asia and Kazakhstan, Vulkanol. Seismol., 4-5, 152-160 (in Russian). THEILER, J., S. EUBANK, A. LONGTIN, B. GALDRIKIAN and J.D. FARMER (1992): Testing for non-linearity in time series: the method of surrogate data, Physica D, 58 ,77-94. TURCOTTE, D.L. (1997): Fractals and Chaos in Geology and Geophysics (Cambridge University Press), 2nd edition. VOLYKHIN, A.M., V.D. BRAGIN and A.P. ZUBOVICH (1993): Geodynamic Processes in Geophysical Fields (Nauka, Moscow). WEBBER, C.L. JR. and J.P. ZBILUT (1994): Dynamical assessment of physiological systems and states using recurrence plot strategies, J. Appl. Physiol., 76, 965- 973. ZBILUT, J.P. and C.L. WEBBER JR. (1992): Embeddings and delays as derived from quantification of recurrence plots, Phys. Lett. A, 171, 199-203. ZHANG, X. and N.V. THAKOR (1999): Detecting ventricular tachicardia and fibrillation by complexity measure, IEEE Trans. Biomed. Eng., 46 (5), 548-555. (received December 10, 2005; accepted March 2, 2006)