Evidence of non-extensivity in the seismicity observed during the 2011–2012 unrest at the Santorini volcanic complex, Greece

During the period of October 2011–January 2012, an increase of earthquake activity has been observed in the volcanic complex of Santorini Island, Greece. Herein, the magnitude distribution of earthquakes as well as the temporal distribution of seismicity are studied. The statistics of both parameters exhibit complexity that is evident in the frequency-magnitude distribution and the inter-event time distribution, respectively. Because of this, we have used the analysis framework of non-extensive statistical physics (NESP), which seems suitable for studying complex systems. The observed inter-event time distribution for the swarm-like earthquake events, as well as the energy and the inter-event earthquake energy distributions for the observed seismicity can be successfully described with NESP, indicating the inherent complexity of the Santorini volcanic seismicity along with the applicability of the NESP concept to volcanic earthquake activity, where complex correlations exist.


Introduction
The Hellenic Volcanic arc lies in the back arc region of the Hellenic subduction zone.It is parallel to the Hellenic arc and about 150 km north of it (Sakellariou et al., 2010).The formation of the volcanic arc is the result of the Hellenic orogen (Dimitriadis et al., 2009).Sousaki, Methana, Milos, Santorini and Nisyros form volcanic centers along the South Aegean Volcanic arc (Fytikas et al., 1985;Dimitriadis et al., 2009).
Thera (Santorini) is the most active volcano of the Aegean volcanic arc (Bond and Sparks, 1976) (Fig. 1).Its volcanic activity began at approx.0.6Ma ago (Perissoratis, 1995).The present-day caldera was formed during the "Minoan" eruption, at approximately 1610-1650 BCE (Bond and Sparks, 1976).After the Minoan eruption, the volcanic activity concentrated at the central part of the caldera complex, where the emission of calc-alkaline magma has led to the formation of the Palaia andNea Kameni islets between 197 BC and1950 AD (Fytikas et al., 1990;Druitt et al., 1999;Stiros et al., 2010).Earthquake activity at and around the Santorini volcanic complex is manifested with small magnitude earthquakes and low seismicity rates.The majority of earthquake foci are concentrated around the Columbos volcano, with fairly low activity observed at the rest of the complex (Bohnhoff et al., 2006;Dimitriadis et al., 2009).At present, it is not clear which part(s) of this activity is of tectonic or of volcanic origin.
Concerning the tectonics, the northern part of the island complex lies in a graben, probably the continuation of the Anhydros Basin (Perissoratis, 1995).Most of the effusive activity since 530 ka, which includes the Peristeria Volcano, the Simandiri shield, the Skaros shield, the Therassia dome complex, and the Kammeni Volcano has been associated with this feature.The contemporary volcanism is manifested by two major volcanic centers, the Nea Kameni volcano rising at the center of the caldera and the submarine Columbos volcano, located 7 km NE of Cape Coloumbo.Their evolution has been affected by two distinct NE-SW tectonic lineaments, the Kameni and Columbos fault zones, respectively (Druitt et al., 1989(Druitt et al., , 1999)).These mark the alignment of several eruptive vents and have been interpreted to comprise major normal faults (Pe-Piper and Piper, 2005).Six Plinian eruptions were aligned along Kameni Fault Zone.Independent volcanic centers at North Thera, as is the Megalo Vouno cinder cone, the Kokkino Vouno cinder cone and the  Cape Columbos tuff ring define the Columbos Fault (for details see Fouqué, 1879 andReck, 1936).In addition, several dykes located at northern Thera, have a NE-SW orientation, e.g. the dyke between Mikros Profitis Ilias and Megalo Vouno (Heiken and McCoy, 1984;Mountrakis et al., 1998).
The southern half of the island is situated at the northern flank of a NE-SW-trending basement horst, the Santorini-Amorgos Ridge (Perissoratis, 1995).
During the period of October 2011-January 2012, an increase of intra-caldera earthquake activity has been observed (Chouliaras et al., 2012; ISMOSAV, http://ismosav.santorini.net -also see Figs. 1 and 2).This activity occurred simultaneously with ground deformation, probably the most significant since the 1950 eruption (Newman et al., 2012).Such observations could indicate resurgence of intra-caldera volcanism after a quiescence of about 6 decades, albeit without assurance that an eruption would be forthcoming, as suggested by Newman et al. (2012).
Recently non extensive statistical physics (Tsallis, 2009) becoming a challenging framework for geophysical complex phenomena (Vallianatos and Telesca, 2012).The main motivation of our work is to establish a solid pathway for the analysis of seismo-volcanic effects using concepts of-non extensive statistical mechanics.Our aim is not to develop a precise model but rather to give a simple argument of physical plausibility.In the present work, we study some physical characteristics of the observed seismicity by analyzing the statistical physics of its magnitude and temporal (interevent times) distributions.The Frequency-Magnitude distribution is investigated by means of the Gutenberg-Richter law.Considering the complex pattern of seismicity in volcanic areas, we use the non-extensive statistical physics concept (NESP; Tsallis, 2009), which is a generalization of Boltzmann-Gibbs statistical physics, to study the energy and temporal patterns of the observed seismicity.Non-extensive statistical physics seems to be a suitable framework for the study of non-equilibrium systems, exhibiting scale invariance, multi-fractality and long-range interactions (e.g.Vallianatos et al., 2012a).Fracture related phenomena, such as seismo-volcanic, present such characteristics (Vallianatos, 2009).It is precisely such phenomena that constitute the scope of non-extensive statistical physics.

Frequency-magnitude distribution
The analysis of the temporal evolution of the volcanic seismicity that occurred in Santorini volcanic complex during the period of October 2011-January 2012 we based on the earthquake catalogue of the Geodynamic Institute of the National Observatory of Athens reporting earthquakes in the M L magnitude scale (http://www.gein.noa.gr/).The seismicity rate at the Santorini area exhibits a systematic increase during this period, as it is depicted in Fig. 2a.The highest rate is observed on 23-24 January 2012 and is followed by a rather extended period of lower seismicity rates (Fig. 2b).
The earthquake size distribution can be described by a power-law relationship over a broad range of magnitudes, which is commonly referred to as the Gutenberg-Richter law (Gutenberg and Richter, 1944): where N is the cumulative number of earthquakes with magnitude ≥ M and a, b are constants.The so-called "b-value" is the slope of the frequency-magnitude distribution (FMD) and describes the size distribution of the earthquake events.While for tectonic earthquakes b has typical values close to 1, higher b-values (up to 3) have been reported for volcano related earthquakes (McNutt, 1996).High b-values associated with volcanic areas have been attributed to crustal heterogeneities (Mogi, 1963), high thermal gradients (Warren and Latham, 1970) and high pore-pressures at the vicinity of the magma chambers, all of which can cause a decrease of the effective normal stress (Wyss, 1973).
The FMD for the given data set has a bimodal character (Fig. 3).The maximum likelihood solution of the Gutenberg-Richter law according to Aki (1965) as was later on revised by Utsu (1978), gives a b-value of 0.76 ± 0.04 that holds for the magnitude range 1.2-2.6.It is also evident that this approximation holds down to a magnitude of 1.2, which is thus adopted to be the magnitude of completeness (M c ) of the given data set.As we can see from Fig. 3, the slope of the FMD increases to b = 1.7 ± 0.36 for magnitudes greater than or equal to 2.6, thus rendering the distribution bimodal.Bimodal distributions in the FMD of volcanic areas have been attributed to earthquake swarms that are related to magma movements (Wiemer and Wyss, 2002).This may be the case for the observed intra-caldera seismicity during the period considered in this study.On the other hand, the SW-NE hypocenter distribution of the earthquakes along an almost vertical fault beneath the Kameni islets (see Newman et al., 2012) indicate that seismicity may also have a strong tectonic component and that the causative magma and fluid movements may actually affect and excite instability on a tectonic fault as well as in the pre-existing fractured surfaces and cracks.

A Non-extensive statistical physics approach
In this section we use the non-extensive statistical physics concept (NESP) to analyze the temporal and magnitude distributions of the observed seismicity for the period 1 October 2011-27 January 2012.The NESP concept refers to the non-additive entropy S q (Tsallis, 2009), which is a generalization of Boltzmann-Gibbs entropy S BG and has been frequently used to characterize complex dynamical systems that exhibit scale-invariance, (multi)fractality and long-range interactions (e.g.Gell-Mann and Tsallis, 2004).The Tsallis' entropy S q is non-additive in the sense that it is not proportional to the number of the elements of the system as S BG does.The Tsallis' entropy S q reads as where k B is Boltzmann's constant, p i is a set of probabilities, W is the total number of microscopic configurations and q the entropic index.The latter is a measure of the nonextensivity of the system and for the particular case q = 1 the Boltzmann-Gibbs entropy S BG is obtained.The cases q > 1 and q < 1 correspond to subadditivity and superadditivity, respectively.Recent applications of the NESP concept to solid earth physics (in regional or planetary scale) are mainly focused on seismology (Abe and Suzuki, 2003;2005;Telesca, 2010aTelesca, ,b, 2011;;Telesca and Chen, 2010) using earthquake catalogs from different seismic zones, fault lengths distribution (Vallianatos et al., 2011a;Vallianatos and Sammonds, 2011) and very recently to natural hazards (Vallianatos, 2009), plate tectonics (Vallianatos and Sammonds, 2010), rock physics (Vallianatos et al., 2011b) and geomagnetic reversals (Vallianatos, 2011).
Considering the frequency-magnitude distribution of seismicity, Sotolongo-Costa and Posadas (2004) proposed a model for the earthquake generation mechanism considering the interaction of two rough fault planes and the fragments filling the space between them, where the fragments are produced by the local breakup of the material comprising the fault planes.This interaction is supposed to modulate earthquake triggering.By using the non-extensive formalism, these authors demonstrated the influence of the size distribution of fragments on the energy distribution of earthquakes and introduced an energy-distribution function, which reduces to the Gutenberg-Richter law as a particular case.Silva et al. (2006) revisited the fragment asperity model of earthquakes and derived an energy distribution function which allows determination of the relative cumulative number of earthquakes as a function of magnitude.Moreover, they proposed a scale law between the released relative energy ε and the volume of the fragments r 3 (ε ∼ r 3 ), in agreement with the standard theory of seismic moment (Lay and Wallace, 1995).This model has been recently applied to regional seismicity, covering diverse tectonic regions (Telesca, 2010a,b, Telesca andChen, 2010).Finally, by considering that the earthquake energy ε is related to the magnitude m as m = 2 3 log(ε), Telesca (2011Telesca ( , 2012) ) has introduced a new expression for the cumulative number of earthquakes as a function of the magnitude: This relationship describes from the first principles and in NESP formalism, the cumulative distribution of the number of earthquakes N with magnitude m greater than the threshold M in a seismic region, normalized by the total number of earthquakes.The constant a expresses the proportionality between the released energy ε and the fragment size r and q M is the entropic index.
We have applied the NESP model of Eq. ( 4) to the normalized cumulative magnitude distribution of the Santorini seismicity for the entire magnitude range above the threshold and the result is presented in Fig. 4. The model describes rather well the observed magnitude distribution for q M = 1.39 ± 0.035 and α = 286.6 ± 78.This fine agreement between the model of Eq. (4) and the observed seismicity, is a good indicator of the suitability and the success of the NESP model in recovering the main characteristics of the earthquake dynamics in the case of volcano related seismicity.The latter has been also pointed out by Telesca (2010a), who has successfully applied the NESP model to the earthquake  4).The values for the best fit regression to the data are q M = 1.39 ± 0.035 and a = 286.6 ± 78.The 95 % confidence intervals for q M and α are also plotted (dashed lines).
activity of the two well-known volcanoes of Southern Italy, Vesuvius and Mt.Etna.
Next, we used the cumulative distribution of the lapse time between successive earthquakes (inter-event times), which can provide useful insights on the complexity that characterizes the physical mechanism of the seismicity process, in order to investigate the temporal properties of the observed seismicity.The inter-event time T is defined as T = t (i+1) − t (i) .Then P (> T ) is the cumulative distribution of the probability finding an event greater than the inter-event time T .In Fig. 5, P (> T ) is plotted for earthquake events with M ≥ M c .The observed P (> T ) exhibits an exponential tail that holds for any T > T c = 5.8 h (Fig. 5).In the statistical physics concept, the exponential function represents a solution of the S BG entropy (Eq. 3) that is used to characterize uncorrelated, random processes.The latter implies that P (> T ) exhibits two domains: an uncorrelated one for inter-event times greater than T c and a more complex one for shorter inter-event times that mainly occur at the jerks of increasing seismicity rate (number of events per day greater than 4 -see Fig. 2a).An uncorrelated domain for the larger inter-event times has been also observed by Bak et al. (2002) for Californian seismicity, in which the correlated earthquake events have been interpreted as a sequence of avalanches in a self-organized critical process.
In order to investigate the complexity that characterizes the regime of P (> T ) where T < T c , we used the NESP concept.Along this vein, Abe and Suzuki (2005) investigated the temporal properties of the seismicity in California and Japan and more recently Vallianatos et al. (2012b) investigated the spatiotemporal properties of the 1996 Aigion (Greece) earthquake aftershock sequence.In these studies, the inter-event time distributions P (> T ) were all nicely fitted with the qexponential distribution where the q-exponential function is defined as, q) , for 1 + (1 − q) x ≥ 0 and e q (x) = 0, for 1 + (1 − q) x < 0. (6) The inverse of Eq. ( 5) is the q-logarithmic distribution ln q (P ( T )).These results imply the existence of complex correlations in the temporal evolution of seismicity.After the estimation of the appropriate q that describes the observed distribution P (> T ), the q-logarithmic distribution ln q (P (> T )), where is linear with T (Abe and Suzuki, 2005;Vallianatos, 2011 and references therein).
In order to test the applicability of this approach to the shorter inter-event times, we selected the events that occurred on 23 and 24 January 2012.During these two days, a sudden increase of seismicity was observed (Fig. 2a), with swarmlike characteristics.The earthquake activity for this period also fulfills two other requirements.First, it contains a sufficient number of events (> 50) for M ≥ M c , in order to have as much statistically significant results as possible and secondly comprises, in their majority, inter-event times shorter than T c .In Fig. 6, the P (> T ) for this period and for the events with M ≥ M c is plotted.As evident, this cumulative distribution can be nicely fitted with the q-exponential distribution when q T = 1.52.This result implies that non-extensive thermodynamics apply when the volcano-related seismicity takes on  5).The value of q for the best fit regression is qT = 1.52.Inset: the q-5 logarithmic distribution for T ≤ Tc, exhibiting a correlation coefficient of ρ = 0.9614.The solid line represents the q-exponential distribution of Eq. ( 5).The value of q for the best fit regression is q T = 1.52.Inset: the q-logarithmic distribution for T ≤ T c , exhibiting a correlation coefficient of ρ = 0.9614.a swarm-like character and that complex correlations exist among the earthquake events during such seismicity bursts.
Moreover, we have used the NESP concept to investigate the inter-event energy distribution, using the approach of Caruso et al. (2007).Earthquakes have often been considered an example of self-organized critical (SOC) phenomenon (e.g.Bak and Tang, 1989;Sornette and Sornette, 1989;Bak et al., 2002), where in a critically stressed crust avalanches of earthquake events can occur on all scales (Main, 1996).Caruso et al. (2007), by considering the avalanche size differences of the dissipative Olami-Feder-Christensen model (OFC -Olami et al., 1992) on a small world topology and the energy differences between actual tectonic earthquakes, have shown that their probability density functions deviate from the normal Gaussian distribution and exhibit a rather generalized Gaussian, namely the q-Gaussian distribution that has been derived in the NESP concept (e.g.Gell-Mann and Tsallis, 2004).According to Caruso et al. (2007), these findings support the hypothesis that a self-organized critical mechanism is the origin of the seismic events.
Following this approach, we investigated the probability distribution function of the energy differences between the successive earthquakes for the volcano-related seismicity considered in this study.We considered the quantity S = exp(M) to be a measure of energy in accordance with the work of Caruso et al. (2007), in order to relate the earthquake magnitude M with the energy released during the earthquake (the latter being an exponential function of the magnitude M).Considering this quantity, we defined the function R(t) = S(t + ) − S(t), where is a discrete time interval and we used it to describe the differences in the energy released between successive earthquakes.Caruso et al. (2007) www.nat-hazards-earth-syst-sci.net/13/177/2013/Nat.Hazards Earth Syst.Sci., 13, 177-185, 2013  The data is well fitted by a q-Gaussian curve (solid line) for the value of the entropic index q R = 2.24 ± 0.09 (95 % confidence intervals -dotted curves).
In Fig. 7 we plot the probability density function (PDF) of R(t) for all earthquakes with M ≥ M c , normalized to zero mean and unit variance in accordance with the normalization condition p(R)dR = 1. Figure 7 indicates that the PDF of R(t) deviates from the standard Gaussian shape due to the existence of heavy tails and can rather be described by the q-Gaussian function of the form p(x) = A 1 − (1 − q) x 2 /B 1/(1−q) , where x = (R − < R >)/σ R , with σ R being the standard deviation of R. For q = 1, p(x) reduces to the Gaussian (normal) distribution (e.g.Gell-Mann and Tsallis, 2004).For q R = 2.24 ± 0.09, we can see that the q-Gaussian function fits the observed PDF of the normalized R(t) very well, confirming the results of Caruso et al. (2007) and implying that a self-organized critical mechanism may be at the origin of the observed volcano-related seismicity.Furthermore, the probability distribution of the avalanche and earthquake sizes (S) can be described by a power law that includes the term S −τ (Caruso et al., 2007).Based on the probability distribution function of S = exp(M) of Santorini seismicity with M ≥ M c , a power-law exponent τ S = 3.0 ± 0.37 is obtained (Fig. 8).Caruso et al. (2007) suggested that an interrelation between the power-law exponent τ S and the q-value of the q-Gaussian function that best describes the PDF of R(t) exists.On the hypothesis that there is no correlation between the sizes of two events, the probability of observing R(t) is given by (Caruso et al., 2007):

CE1
The correct is , we can state ….. where K is a normalization factor, ε is a small positive value and 2 F 1 is the hypergeometric function.As pointed out in Caruso et al. (2007), the probability distribution P (R), which depends on τ S , can be approached by means of q-Gaussians, whose values of q R do not depend on ε.Then, by applying Eq. ( 8) for various τ S values, Caruso et al. (2007) found the following stretched exponential interrelation between q R and τ S : According to Eq. ( 9), the value of the power-law exponent τ S = 3.0 ± 0.37 obtained for the observed seismicity, corresponds to q R ≈ 1.64 if no correlations are taken into account, which is dissimilar to the value of q R found by fitting the q-Gaussian to the data set at hand.We would like to note that Eq. ( 9) is valid for uncorrelated events, which is not the case of seismicity as has been also pointed out by Caruso et al. (2007) who stated that "temporal correlations among avalanches (earthquakes) do surely exist and a certain degree of statistical predictability is likely possible", and thus considering the applicability of Eq. ( 9) with caution.
On the other hand, Bakar and Tirnakli (2009) analyzed numerically the Ehrenfest's dog-flea SOC model and obtained, with high accuracy, the value of τ S = 1.517 for the powerlaw exponent, along with q R = 2.35 for the q-Gaussian function, which is in agreement with Eq. ( 9) and in good agreement, taking into account experimental errors, to the value of q R obtained for the Santorini seismicity.The result of Bakar and Tirnakli (2009) has been achieved using a simple prototype SOC model (different from that used by Caruso et al., 2007) and as noted by Sarlis et al. (2010), who suggested that q R = 2.35 is the value that one should use in the q-Gaussian to check whether P (R) resulting from earthquake catalogs can be approached, emphasizing to the fact that the result of Bakar and Tirnakli (2009) is much more general and the q R value is not a fitting parameter anymore.
In addition, in a recent work Celikoglou et al. (2010) propose instead of self-organized criticality mechanism (Bak and Tang, 1989) the application of the coherent noise model (CNM) (Newman, 1996) as a simple and robust mechanism to describe seismicity characteristics.This approach validated very recently using the Olami-Feder-Christensen model (Zhang et al.,201l).The CNM model is based on the notion of an external stress acting coherently onto all agents of the system without having any direct interaction with them, leading to a power-law distribution of event sizes without exhibiting criticality.Applying such a model Celikoglu et al. (2010) propose that on CNM principles that q R = τ S + 2 τ S , which for τ S ≈ 3 gives q R ≈ 5/3 a value slightly different from the approximated relation presented by Caruso et al. (2005) (see Eq. 9).Furthermore, in the frame of CNM presented by Celikoglu et al. (2010) a possible source of the dissimilarity between the observed and the estimated values of q R could be the limited number of events (i.e. a size effect).In fact, the coherent noise model has tested by Celikoglu et al. (2010) for data that span 6 decades providing very clear q-Gaussian curves for the return distributions (∼ 8 decades) in contrast to the limited range of present data.This could be one of the possible reason for the disagreement between the experimentally estimated q R and the calculated one using τ S .The later implies that when the number of events is small, the size distribution has a very short powerlaw region, which causes the distribution of R(t) to deviate from the q-Gaussian with the correct q-values, but keeping its shape.We note that since it is generically extremely difficult (if not impossible) to achieve a size-independent case (which theoretically demands an infinity number of events), the size effect, we can state that for seismo-volcanic systems, even in the size dependent case, the q-Gaussian could still describe the distribution in a satisfactory way.

Conclusions
In this study we have investigated the magnitude and temporal distribution of seismicity during the recent unrest at the intra-caldera region of the Santorini volcanic complex.The frequency-magnitude distribution departs from the powerlaw distribution at higher magnitudes, exhibiting a bimodal character.This complexity in the earthquake energy distribution is well reproduced by the fragment asperity NESP model for the value of the entropic index q M = 1.39.By considering the differences between the energy released by successive earthquakes, we have seen that they follow a q-Gaussian distribution with an entropic index of q R = 2.24, implying that non-linear dynamics control the evolution of the volcano-related seismicity.Furthermore, we have investigated the temporal evolution of the observed seismicity by considering the inter-event time distribution between the successive earthquakes.The observed distribution exhibits two domains: an uncorrelated one for larger inter-event times (> 5.8 h) and a more complex one for shorter inter-event times.In order to investigate the latter domain, we have considered the inter-event time distribution of the 23-24 January 2012 earthquake activity, when a sudden increase in the seismicity rate has occurred.The observed distribution can be described very well with the q-exponential function for the value of q T = 1.52.The results obtained from the time series analysis indicate that Boltzmann-Gibbs thermodynamics apply in the uncorrelated domain of larger interevent times, whereas for the shorter inter-event times, when a swarm-like sequence appears, complex correlations exist and non-extensive thermodynamics are more appropriate descriptors of the seismicity process.
All the results discussed in this study indicate the inherent complexity of the Santorini volcanic seismicity, the applicability of the NESP concept to volcanic earthquake activity and the usefulness of NESP in investigating complex Earth dynamic systems, and could be in principle used to define organization of volcanic seismicity contributing to hazard estimation.
Figure 1.a) The modern active trenches (thick dark lines with solid barbs) for the 3 Fig. 1.(a) The modern active trenches (thick dark lines with solid barbs) for the Hellenic subduction zone as they are indicated by Royden and Papanikolaou (2011) and the hellenic volcanic arc (black dashed line) with the major volcanoes indicated as red triangles.(b) The seismic events (open circles) of Santorini volcanic complex area for the period 1 October 2011-31 January 2012.

Figure 6 .
Figure 6.Log-log plot of the inter-event times cumulative distribution for the period 3

Fig. 6 .
Fig. 6.Log-log plot of the inter-event times cumulative distribution for the period 23 January 2012-24 January 2012 for M ≥ M c .The solid line represents the q-exponential distribution of Eq. (5).The value of q for the best fit regression is q T = 1.52.Inset: the q-logarithmic distribution for T ≤ T c , exhibiting a correlation coefficient of ρ = 0.9614.

Figure 7 .
Figure 7. Probability density function of p(x) (solid circles) for the events with M ≥

Fig. 7 .
Fig. 7. Probability density function of p(x) (solid circles) for the events with M ≥ M c on a semi-log plot, where x = (R − < R >)/σ R (see text).The dashed curve represents the standard Gaussian shape.The data is well fitted by a q-Gaussian curve (solid line) for the value of the entropic index q R = 2.24 ± 0.09 (95 % confidence intervals -dotted curves).
left column, 9 th line from the top This approach validated very recently using the Olami-Feder-Christensen model (Zhang et al., 201l).10 Page 7, left column, 15 line from the top Applying such a model Celikoglu et al., (2010) propose that

Fig. 8 .
Fig. 8. Probability distribution function of the earthquake energies with M ≥ M c scaled as S = exp(M) in a log-log plot.The solid line represents the power law fit with an exponent τ S = 3 ± 0.37 (95 % confidence intervals -dashed lines).