Investigation of the temporal fluctuations of the 1960 – 2010 seismicity of Caucasus

The time-clustering behaviour of the seismicity of the Caucasus spanning from 1960 to 2010 was investigated. The analysis was performed on the whole and aftershockdepleted catalogues by means of the method of Allan Factor, which permits the identification and quantification of timeclustering in point processes. The whole sequence is featured by two scaling regimes with the scaling exponent at intermediate timescales lower than that at high timescales, and a crossover that could be probably linked with aftershock time activiation. The aftershock-depleted sequence is characterized by higher time-clustering degree and the presence of a periodicity probably correlated with the cyclic earth surface load variations on regional and local scales, e.g. with snow melting in Caucasian mountains and large Enguri dam operations. The obtained results were corroborated by the application of two surrogate methods: the random shuffling and the generation of Poissonian sequences.


Introduction
The investigation of the temporal fluctuations of seismicity represents an important step to perform in order to gain a better knowledge of the dynamics of the seismic process under study.For instance, this issue can be considered as a crucial aspect to consider in the general framework of studies devoted to the seismic hazard, to the feasibility and reliability of the probability estimation of occurrence of future earthquakes, to the knowledge of the statistical distribution of earthquake occurrence.The Poissonian distribution, characterized by an exponential decreasing function of the seismic interevent time, models seismic sequences that appear to be homogeneously distributed on time, and explains some characteristics like absence of memory, absence of correlation and independence of all the earthquakes.But, at both short and long timescales, earthquakes correlate with each other (Kagan and Jackson, 1991).The timecorrelation in seismic sequences is manifested by the property of time-clusterization: seismicity evolves through time changing from phases of high temporal density of events to phases of low temporal density, giving to the whole sequence the feature of ensemble of clusters (Telesca et al., 2000).Time-clustering structures in earthquake sequences were identified by using the coefficient of variation, which is defined as the ratio between the standard deviation and the mean value of the interevent times, or by means of the histogram interevent times, which estimates the probability density function of the interevent time.Both these measures are not sufficient to fully characterize the time-clusterization of a sequence, because the first is not informative of the temporal ranges over which the sequence is time-clusterized, the second is not informative of the correlation properties of the sequence that are linked with statistics of the second order.The identification of clusterized behaviour in seismicity was performed at a regional and local level (Telesca et al., 2001(Telesca et al., , 2009)), revealing also a non-trivial relationship with space (Telesca et al., 2003), time (Telesca et al., 2007) and magnitude (Telesca and Macchiato, 2004).

Study area
The Caucasian segment of the Mediterranean Alpine Belt is located between the still converging Eurasian and Africa-Arabian lithosphere plates and is a typical wide zone of continent-continent collision.The main tectonic units of Caucasus are platform and fold-thrust units, which from north to south are the Scythian (pre-Caucasus) young platform; the fold-thrust mountain belt of the Great Caucasus including zones of the Fore Range, Main Range, and Southern Slope; the Transcaucasian intermountain depression superimposed mainly on the rigid platform zone (Georgian massif); the Achara-Trialeti and the Talysh fold-thrust mountain belts; the Artvin-Bolnisi rigid massif; the Loki (Bayburt)-Karabagh-Kaphan fold-thrust mountain belt; the Lesser Caucasus ophiolitic suture; the Lesser Caucasian part of the Taurus-Anatolia-Central Iranian platform; and the Aras inter-mountain depression at the extreme south of the Caucasus (Fig. 1).According to GPS data, the rate of the convergence is ∼20-30 mm y −1 , of which 2/3 is accommodated in the Sevan-Akera ophiolitic suture and the rest chiefly by crustal shortening inside South Caucasus.The tectonic movement of plates formed a net of active faults.Though the two largest events in Caucasus, Spitak (1988, M = 6.9) and Racha (1991, M = 6.9) earthquakes, occurred within a 14 month interval, the recurrence time of strong earthquakes is large enough.Comparison of model probabilistic seismic hazard (PSH) map with real maps of seismic shaking for the last century shows that the optimal PSH map, which does not miss strong events, is a map for 2 % probability of being exceeded in 50 yr, i.e. for average recurrence time 2500 yr (Chelidze and Javakhishvili, 2003;Matcharashvili et al., 2009).
In this study we present the analysis of the time fluctuations of the 1960-2010 seismicity of Caucasian area, fo-cusing on the time-clustering behavior of the earthquakes.We investigated the seismicity that occurred from 1 January 1960 to 31 December 2010 with depth h ≤50 km.The data were extracted from earthquake catalogue for Caucasus and the adjacent territories of Northern Turkey and Northern Iran available at the Institute of Geophysics and Seismic Monitoring Center, Tbilisi, Georgia.The total number of the events in the catalogue is 63 254 for magnitude ranging between 0 and 7.The completeness of the catalogue was checked by the Gutenberg-Richter law and the number of earthquakes versus magnitude was computed to get the completeness magnitude M c , which in this study is defined as the magnitude at which a power law can model 90 %, or more, of the frequency-magnitude distribution (Wiemer and Wyss, 2000).Figure 2 shows that the completeness magnitude M c of the whole Caucasian catalogue (black circles) is 2.0 and the b-value, estimated by using the maximum likelihood estimation method (Aki, 1965), is 0.677.Since large events activate aftershocks, which increase the time-clusterization of the seismicity, in order to avoid bias in the time-clustering behavior due to the presence of aftershocks, we declustered the catalog by means of the Reasenberg's algorithm (1985) producing a data set of 46 802 events, whose completeness magnitude M c , defined as above, is 2.2 and the b-value is 0.725 (Fig. 2, downward red triangles).From Fig. 2 we can also observe that the magnitude-frequency distributions of the whole and the aftershock-depleted catalogues are almost identical at high magnitudes.The plot shows also the binned frequency-magnitude distribution of the whole (upward black triangles) and aftershock-depleted (upward red triangles) catalogues.

Allan Factor method and data analysis
A seismic sequence is a realization of a marked stochastic point process that describes earthquakes occurring randomly in time and marked by their magnitude.Such representation was also utilized to model diverse point processes like lightning (Telesca et al., 2005), starquakes (Telesca, 2005), solar flares (Telesca, 2007), fires (Telesca and Song, 2011) or anthropic events (Telesca and Lovallo, 2006, 2007, 2008).One of the main properties to investigate in a point process is the time-clustering, which, if present, leads to power-law shaped statistics (Thurner et al., 1997), whose exponent quantifies the degree of clusterization.The method used to detect and quantify the time-clustering in the seismic sequences is the Allan Factor (AF).This method provides the estimation of the scaling exponent in fractal point processes, permitting to discriminate between Poissonian sequences (whose fractal exponent is nearly zero) and clusterized processes (whose fractal exponent is significantly larger than zero).The AF method represents the counterpart in the point process domain of the power spectral density, which is the well-known method to identify and quantify correlation structure in continuous time domain.Therefore, like the power spectral density, if the AF is almost flat for all the timescales, then the point process is memoryless, uncorrelated, indicating that a Poissonian model is well suited to fit the temporal distribution of the events; otherwise, if it behaves as a power law of the timescale, the process is time-clusterized and the numer- ical value of the power-law exponent, which is the fractal exponent, quantifies the strength of the inner time-correlations.
method is based on a counting process representation of the earthquake sequence, in which each event is a Dirac's delta centered on the occurrence time, with amplitude corresponding to the magnitude.The observation period is divided into equally spaced contiguous counting windows of duration T , and a sequence of counts {N k (T)} is produced, where N k (T ) is the number of events falling into the k-th window of duration T .The duration T is called timescale.The AF is given by the variance of successive counts for a specified timescale T divided by twice the count average Since the AF is defined through the difference of successive counts, the nonstationarities, at least of the first order, are removed (Viswanathan et al., 1997).Varying the timescale T , a relationship AF(T)∼T is obtained and the time-clustering of the sequence can be revealed by the power-law, with 0<α<3 (Lowen and Teich, 1995).In this case the seismic sequence is clusterized with degree α.For Poissonian sequences, the AF is approximately flat, with value close to 1 for almost any available timescale.The timescale T 1 is the fractal onset time and indicates the crossover between Poissonian behavior and clusterization.
Figure 3 shows the AF of the whole Caucasian seismic sequence (M ≥ 2.2) for timescales T varying from 10 min to about 5 yr; the upper timescale, which is nearly 1/10 of the whole period, is due to poorer statistics if considering higher timescales.The AF curve shows several features: (i) The process is not Poissonian because the AF, plotted in bilogarithmic scales, is not constant for any timescale, but shows different dynamical regimes with scaling behavior; (ii) the lower timescale region up to approximately 3 h represents the Poissonian regime and the cutoff timescale of about 3 h can be considered as the fractal onset timescale; (iii) a scaling region is visible at intermediate timescales from 3 h to about 2.5 days with scaling exponent ∼0.64; (iv) a second scaling region involves the high timescales from 2.5 days up to the highest available timescales (around 5 yr), with scaling exponent ∼0.5; (v) the presence of the two scaling regions indicates the co-existence of two different time dynamics in the process generating the seismic sequence; it was generally observed that a two-fold scaling regime with the scaling exponent in the first region (intermediate timescales) higher than that of the second region (high timescales) is indicative of aftershock effects (Currenti et al., 2005).
Two surrogate methods have been applied to check the significance of the results generating (i) one hundred Poissonian sequences with the same number and the same mean intervent time of the original sequence, and (ii) one hundred shuffled sequences with the same interevent interval probability density function as the original one.To each surrogate (Poissonian or shuffled) the AF method was applied, and for each timescale the 95th percentile among the AF values of the surrogates for that timescale were computed.The 95 % confidence AF curve is given by the set of the 95th percentiles.Figure 3 shows also the 95 % confidence AF curves for Poissonian (red) and shuffled (green) surrogates.The AF curve of the original sequence is significantly above both the 95 % confidence AF curves, indicating that the Caucasian sequence is not a realization of a Poisson process and that the time-clustering depends on the ordering of the interevent intervals.
Figure 4 shows the results of the AF method applied to the declustered catalogue (M ≥ 2.2).Similarly to the whole sequence, several features can be evidenced: (i) The declustered sequence is not a realization of a Poisson process, in fact, not only the AF is not flat for any timescale, but is also significantly different from the 95 % confidence curve obtained for 100 Poissonian surrogates (red); (ii) two scaling regions are visible: one involving the timescales from about 14 days to 6 months with scaling exponent ∼1.53 and the other with scaling exponent ∼2.45 involving the timescales higher than 1.25 yr; (iii) a clear drop of the AF curve is around 1.25 yr, indicating the presence of a cycle in the temporal fluctuations of the seismicity with period around 1.25 yr; such possible cycle could be linked with the time fluctuations of the seismicity that could be induced by different causes, e.g. by the influence of Enguri high dam water reservoir (Telesca, 2011).Possibility of such influence of quasi periodic reservoir water level variation on the local seismicity in western Caucasus was already described (Matcharashvili et al., 2008).On regional scales presence of cycles in seismicity with around one year period can also be explained by seasonal variations of surface loading (snow melting and precipitation) in mountainous Caucasian region.Possibility of such yearly variations were documented for seismicity in Himalayas (see Fig. 3, Bettinelli et al., 2008); (iv) the high similarity between the AF of the original declustered sequence and 95 % confidence AF curve for the shuffles (green) indicates that the time-clustering of the seismicity is due to the shape of the probability density function of the interevent times and does not depend on their ordering, contrary to the whole sequence; (v) the appearance of the two scaling regions could be due to the drop at about 1.25 yr, indicative of the reservoir-induced seismicity.

Conclusions
The time-scaling properties of the 1960-2010 seismicity of the Caucasus was investigated by using the statistical method of Allan Factor that permits the detection and quantification of time-clustering in point processes.The whole and aftershock-depleted catalogues show very different scaling behaviour: The whole sequence is characterized by a twofold scaling behaviour with the scaling exponent at intermediate timescales lower than that at high timescales, with a crossover timescale probably linked with the aftershock time activation; the aftershock-depleted sequence reveals a striking higher time-clustering degree with the presence of a possible periodicity, which could be linked with the cyclic dam operations of Enguri.The surrogate analysis indicate that the whole catalogue is significantly non-Poissonian and that the time-clustering behaviour depends on the orderings of the interevent times.Concerning the aftershock-depleted catalogue, the surrogate method suggest that it is not modelled by a Poisson process, but the time-clustering behaviour depends on the probability density function of the interevent times.

Fig. 1 .
Fig. 1.Schematic map of the seismic sources of Caucasus.

Fig. 2 .
Fig. 2. Cumulative frequency-magnitude of the whole (black circles) and aftershock-depleted (downward red triangles) catalogues.The plot shows also the binned frequency-magnitude distribution of the whole (upward black triangles) and aftershock-depleted (upward red triangles) catalogues.