Natural Hazards and Earth System Sciences Recent seismicity detection increase in the Santorini volcanic island complex

Santorini is the most active volcanic complex in the South Aegean Volcanic Arc. To improve the seismological network detectability of the seismicity in this region, the Institute of Geodynamics of the National Observatory of Athens (NOA) recently installed 4 portable seismological stations supplementary to the 3 permanent stations operating in the region. The addition of these stations has significantly improved the detectability and reporting of the local seismic activity in the NOA instrumental seismicity catalogue. In this study we analyze quantitatively the seismicity of the Santorini volcanic complex. The results indicate a recent significant reporting increase mainly for events of small magnitude and an increase in the seismicity rate by more than 100 %. The mapping of the statistical significance of the rate change with the z-value method reveals that the rate increase exists primarily in the active fault zone perpendicular to the extensional tectonic stress regime that characterizes this region. The spatial distribution of the b-value around the volcanic complex indicates a low b-value distribution parallel to the extensional stress field, while the b-value cross section of the volcanic complex indicates relatively high b-values under the caldera and a significant b-value decrease with depth. These results are found to be in general agreement with the results from other volcanic regions and they encourage further investigations concerning the seismic and volcanic hazard and risk estimates for the Santorini volcanic complex using the NOA earthquake catalogue.

In order to improve the detection of the local seismicity in Santorini, the Institute of Geodynamics of the National Observatory of Athens (NOA) in cooperation with the Geophysical Laboratories of the University of Athens and Thessaloniki installed a network of 8 telemetric seismological stations in 1994.The seismological data from this network identified two clusters of epicenters, one in the Kameni caldera and one in the Columbo submarine volcano to the northeast (Panagiotopoulos et al., 1996).An increase in the seismic activity was observed in 1995 and this activity reached a maximum on 27 May 1996 with a swarm of shocks with magnitudes M = 3.8−4.2.Stavrakakis et al. (1996) determined the source spectra of the strongest shocks in this sequence and the seismic source parameters indicate a constant source radius irrespective of the source energy.This result is consistent with similar results from other volcanic areas due to the volcanic processes.
The first permanent seismological station on Santorini was installed by NOA in 1997 and in 2010 and 2011 two permanent seismological stations were also installed on the nearby volcanic island of Milos (http://bbnet.gein.noa.gr/NOA HL/index.php/real-time-plotting/hl-network).In addition to the permanent seismological stations, NOA installed a portable seismological network comprising of four broadband stations on Santorini in June 2011.This portable network transmits in real time waveform data that are incorporated in the everyday analysis and the earthquake catalogue of NOA thus significantly increasing the reporting of small seismic events (http://bbnet.gein.noa.gr/NOAHL/ index.php/real-time-plotting/santorini-network).Earthquake catalogues are not simple lists, but rather complex data sets containing both natural and man-made changes.During any period, the seismologists responsible for a seismic network may improve the ability to detect or to locate events more accurately, by adding more stations, improving signal to noise conditions, improving the signal processing and the magnitude calculation and they may also proceed to uninstall certain stations.These manmade artifacts in earthquake catalogues are apparent changes in the seismicity rate and mask the true estimation of tectonic seismicity patterns and the seismic hazard and risk estimates (Habermann, 1982(Habermann, , 1987(Habermann, , 1991;;Habermann and Craig, 1988;Zuninga and Wiemer, 1999;Tormann et al., 2010).Habermann (1982) was first to introduce a quantitative technique for identifying man-made artifacts in cata- logues with the magnitude signature method, which proved to be a valuable discriminator between natural and artificial variations in seismicity.The method efficiency to identify detection changes, reporting changes and systematic changes in magnitudes has been tested for several earthquake catalogues (Habermann, 1987(Habermann, , 1991)).More recent techniques for detecting artificial seismicity in earthquake catalogues use the spatial mapping of seismicity rate and frequency magnitude distribution (Wiemer and Wyss, 1994;Zuninga and Wyss, 1995;Zuninga and Wiemer, 1999;Zuninga et al., 2005).
In this publication, the recent increase in the seismicity of Santorini is analyzed using the NOA earthquake catalogue to evaluate whether the detection increase is due to natural (earthquake preparation) or due to artificial (man-made) causes.

Data analysis and results
In this study we investigate quantitatively the seismicity catalogue of NOA (http:www.gein.noa.gr/services/cat.html)for the Santorini region using the ZMAP/IASPEI software package, in order to present reproducible results (Wiemer, 2001).
The seismicity map of the region surrounding the volcanic complex (Fig. 2) reveals two separate clusters of seismic events: one in the caldera and one just to the north with a NE-SW epicentral distribution, at the Kameni and Columbo faults, respectively.This data set contains 911 events from 1966 until the end of 2011, with a magnitude range M L (local) = 1-5 and hypocentral depths from 1 to 167 km.
The blue curve in Fig. 3 shows the time evolution of the cumulative seismicity for this region.A power law increase in the reported seismicity is observed in general with the NOA seismicity catalogue, due to the expansion of the seismological network and the improvement in detecting seismic events in recent years (Chouliaras, 2009).The sudden increases in the seismicity for the Santorini region are due to the occurrence of seismicity swarms and also due to the seismic sequence of the largest recent earthquake on 26 June 2009 with a magnitude M L = 5 to the NE of Thera island along the Columbo fault line.
In order to study seismicity rates, one must separate dependant events such as swarms and aftershocks from the independent events.For this reason we apply Reasenberg's (1985) declustering algorithm to the NOA earthquake catalogue for the Santorini region, in order to determine the declustered cumulative curve (marked in red) in Fig. 3.This curve contains 532 events and appears much smoother, thus, indicating that the earthquake clusters have been removed effectively.
It is apparent from the clustered catalog cumulative curve in Fig. 3 that more than 30 % of the events were detected in the period after the year 2010 and this rate change is investigated with the method of Zuninga and Wiemer (1999).The seismicity rate of the declustered catalogue for a background period from 2005 until 2010 is compared to the rate for the recent period from 2010 until the end of 2011 (14 December 2011) in Fig. 4a and b.By comparing the cumulative and the noncumulative frequency-magnitude distribution curves (normalized by the duration of the time periods), we determine a rate increase of more than 100 % for the recent period and this increase is noticed mainly for events with M L < 3. The magnitude shift between the two periods for magnitudes with M L < 3 (Fig. 4b) may be due to the improvements in NOA magnitude reporting procedures as of February 2011.Since the rate increase for the recent period is found to be magnitude-dependant, the magnitude signature method (Habermann, 1987) is employed in order to analyze the results in the magnitude domain.Magnitude signatures determine the significance of the rate change between two time periods as a function of magnitude.In order to determine the difference between the means of each period, the most general of the parametric statistical tests is applied to the data, namely the z-test.The two rates are compared using the formula Eq. ( 1) that defines the normal deviate z as (1)  where, M1, M2 are the mean rates during the two periods, S1, S2 the respective statistical deviations and N1, N2 the respective number of samples.The resulting z-value is interpreted in terms of significance, as the number of standard deviations from the mean of a standard distribution (i.e.z = 1.64 is 90 % significance, z = 1.96 is 95 % significance, z = 2.57 is 99 % significance).
Figure 4c shows the magnitude signature for the two periods of investigation.The z-value in the vertical axis has positive values in the upper part of the plot, thus, indicating a rate decrease, whereas the opposite holds for the lower part of the plot with negative values indicating a rate increase.The horizontal axis indicates the magnitude band for which the calculation is made.Magnitude bands on the left side include events with magnitudes less than those shown on the axis and magnitude bands on the right side include events with magnitudes greater than those shown on the axis.
There are different characteristic magnitude signatures for each type of seismicity rate change.Detection changes, reporting changes and magnitude shifts have different and distinct characteristic appearances on a magnitude signature and for this reason, our result in Fig. 4c is compared with actual and theoretical magnitude signatures.This comparison indicates that the Santorini magnitude signature is in good agreement with the magnitude signatures determined from periods of detection increases due to the additional installation of seismological network stations (see Figs. 2, 3 in Habermann, 1987).Characteristic for this increase are the negative z-values for smaller events on the left side of the plot and less significant change (smaller z-values) for the larger events on the right side of the plot.Another interesting result in Fig. 4c is the occurrence of z-values of opposite signs on the opposite sides of the magnitude signature, indicative of a magnitude shift for events of small magnitudes (M < 3.5) (see Figs. 11, 16, 17 in Habermann, 1987).
In order to investigate the spatial dependence of the rate change shown in Figs.4a-c, we mapped the z-value for the Santorini region using the gridding technique of Wiemer and Wyss (2002) in Fig. 5.A significant rate increase (negative z-values) is observed over Kameni and Thera islands and this increase is elongated in a NE-SW direction from the Columbo fault, almost perpendicular to the NNW-SSE extensional stress regime of the region.The pattern resembles a Coulomb stress pattern and this agrees with the results of more detailed studies that show that the pattern of increased or decreased seismicity accordingly follows the advancement or retardation of faulting by the Coulomb fracture criteria (Stein et al., 1992;Wyss and Wiemer, 2000).
The inverse relationship between the b-value and the stress regime has been validated by laboratory tests on rock specimens (Mogi, 1962;Scholz, 1968), as well as by quantitative studies of earthquake catalogues from different regions of the world with different stress regimes (Smith, 1981;Hauksson, 1990;Frohlich and Davis, 1993;Wiemer and Wyss, 1997;Kagan, 2002;Schorlemmer et al., 2004;Schorlemmer and Wiemer, 2005a, b).In volcanic regions a characteristic high b-value has been reported and this is attributed to the phenomenon of creep in the magma chamber due to the migration of fluid or magma (Wyss et al., 1997(Wyss et al., , 2001;;Wiemer et al., 1998).
Assuming that earthquake magnitudes obey the Gutenberg-Richter law (Gutenberg and Richter, 1944), we determined magnitude of completeness (M c ) and the b-value (b) for the declustered seismicity catalogue of Santorini by the use of the weighted least squares (WLS) fitting method and the FMD curve as shown in Fig. 6 (Wiemer, 2001).The magnitude of completeness (M c ) is defined as the lowest magnitude, at which all earthquakes in a space -time volume are reliably detected by a seismological network (Rydelek and Sacks, 1989;Taylor et al., 1990;Wiemer and Wyss, 2000;Woessner and Wiemer, 2005).M c is an important parameter for statistical investigations in seismology and it is traditionally determined from the maximum point of the derivative of the FMD (marked in green) curve in Fig. 6 (Wiemer, 2001).This procedure yields a value of M c = 3.4 and b = 1.47 for the declustered earthquake catalogue of the Santorini region.
Changes in the detection and reporting of a seismological network may cause the M c and the b-value to vary, so it is important to study the temporal dependence of M c for the Santorini area before and after the observed rate increase shown in Fig. 4. To determine the M c as a function of time for the NOA declustered catalogue, we employed the method of Wiemer (2001) and in Fig. 7 we observeed a M c variation from 3.4 to 3.6 until the year 1994.
In 1994, the NOA seismological network of 18 analog short period stations was complemented with 9 additional permanent digital short period stations and after 1997 the analog seismological network began to be upgraded with new installations of digital broadband stations.Presently, the NOA Hellenic Broadband Seismological Network (NOA HL) consists of 56 stations that provide real time waveform and parametric data (http://bbnet.gein.noa.gr/NOA HL/index.php/real-time-plototing/hl-network).The seismological network station increase was complemented with improvements in the analysis software and in the reporting procedures as of February 2011 and this has resulted in a significant improvement in the detection of seismic events in Greece.This improvement may be seen in Fig. 7, in which we observe a decrease in the M c to a value of M c = 3 before the year 2010, and after that a very sharp decrease to M c = 2 towards the end of 2011.
The b-value exhibits heterogeneities in time and space, depending on the stress patterns and in this way it acts as a stress meter (Schorlemmer and Wiemer, 2005b).For the Santorini volcanic complex, we mapped the spatial distribution of the b-value using the declustered NOA catalogue with M c > 3.5 from 1966 to the end of 2011 and the gridding technique of Wiemer and Wyss (2002).In Fig. 8, we note a bvalue variation from 0.85 to 1.35 and a zone of low b-values (appearing in blue) that trends NNW-SSE, crossing the volcanic islands and separating the two zones of higher b-values (appearing in red) to the East and to the West.This alignment of the low b-value zone is found to follow the NNW-SSE extensional stress field orientation for this region (Fytikas et al., 1990;Vougioukalakis et al., 1995;Mountrakis et al., 1996).
Concerning the Santorini region, an increase in the reported events in the NOA seismicity catalogue occurred after the installation of two permanent seismological stations on the nearby island of Milos in 2010 and 2011, and also after the installation of the portable network of four broad band seismological stations in June 2011.The significant detection improvement results in a change of the M c and this change will result in a different b-value distribution map for this region for the recent period.
To investigate the recent b-value distribution, we used the NOA catalogue to determine the spatial dependence of the b-value in the cross section of the main seismic activity in Fig. 9a.The cross section of the main seismic activity contains 160 events with M c >2 between 2010 and the end of 2011 in a NNE-SSW distribution, parallel to the faulting orientation and perpendicular to the extensional stress field.We found the main clusters of seismic activity in the Kameni and Columbo active faults, while on either side away from the volcanic complex, we observed a more diffused seismicity pattern.
Figure 9b shows the b-value cross section to a 40 km depth range and this reveals a significant lateral and vertical heterogeneity around and beneath the volcanic complex, with b-values ranging from 0.7 to 1.3.Relatively higher b-values (0.9-1.3) are observed under the volcanic complex to a 15 km depth and lower b-values (0.7-0.9) are observed laterally on either side.The high b-values near the surface show a decrease with depth in an inverted cone distribution.This bvalue distribution is generally in agreement with results from other volcanic areas (Mt.St. Hellens, Mt.Spurr, Tokohuku, Kilauea) that report higher b-values (up to b = 2) near the surface, due to the low ambient stress and the highly frac-tured environment around the magma chamber.The decrease of the b-value with depth is mainly due to the increases in stress, pore pressure and the thermal gradient (Wiemer and McNutt, 1997;Wiemer and Wyss, 1997;Wyss et al., 2001).

Discussion and conclusions
In this study we used the NOA seismicity catalogue to quantify valuable parameters such as the magnitude of completeness and the b-value from the frequency-magnitude relation, as well as the recent seismicity rate change (z-value) for the Santorini volcanic complex.
Many studies have reported anomalies in the seismicity rate during earthquake and volcanic preparation processes.To reveal possible spatial and temporal rate anomalies, one must first carefully analyze the employed earthquake catalogue and identify any possible artificiality due to man-made changes.This artificiality exists in all catalogues and it depends on the increases or decreases in the detection, resulting from the installation or removal of seismic stations as well as to possible changes in the operative procedures and reporting of the seismological network operator, which may change with time.
As part of the detection improvements of the NOA seismological network in the vicinity of the Santorini volcanic island complex, a local network of 4 portable broad band seismic stations was installed in June 2011 to supplement the first permanent station that was installed on Thera island in 1997.In 2010 and 2011 two other permanent stations were also installed on the nearby volcanic island of Milos and as a result the detectability and reporting of the local seismicity in the NOA seismicity catalogue improved significantly.
The NOA seismicity catalogue contains more than 900 events for the region surrounding the Santorini islands, from 1966 till the end of 2011 and more than 30 % of these events were detected after the installation of additional seismological stations in 2010 and 2011.
Quantitative analysis of the normalized cumulative and noncumulative FMD curves, for the period before and after the new seismological station installations, reveals a seismicity rate increase of more than 100 % for the investigated region.This increase concerns mainly small magnitude (M < 3) earthquakes and is also confirmed by the negative zvalue determined by the magnitude signature method.
The significance of the rate change is mapped using the zvalue method and we found the rate increase to mainly exist in the Kameni and Columbo faults in a NNE-SSW direction parallel to the active faults and perpendicular to the NNW-SSE extensional stress direction.
The frequency-magnitude distribution of the NOA seismicity catalogue for the region surrounding the volcanic complex determined bulk values for M c = 3.4 and b = 1.47 for the entire declustered catalogue.Improvements in the NOA seismic network detectability, by the increase of local seismic stations and the improved magnitude reporting, caused the M c and the b-value to vary with time, from a value of M c = 3.5 in 1966 to M c ≈ 2 at the end of 2011.
Using the declustered catalogue data from 1966 till the end of 2011, we determined a b-value map that indicates a low b-value zone across the volcanic complex, in a NNW-SSE direction, parallel to the dominant extensional stress, and higher b-values on either side.
The recent seismological data after the additional station installations in 2010 and 2011 have significantly lower M c (>2) and 160 of these events were used to determine a bvalue cross section along the direction of the principal seismic activity of the region (NNE-SSW).The b-value cross section across the volcanic complex indicates prominent vertical and lateral heterogeneities.The distribution pattern of b-values map the Santorini volcanic complex with relatively high b-values in the top fifteen kilometers and decreasing bvalues with depth and lateral distance (on either side of the complex).
The recent b-value distribution on Santorini is in qualitative agreement with reported b-value maps determined in other volcanic regions, in which mapping of the magma chamber was attempted for earthquake hazard and risk estimation and eruption forecasting.Nevertheless, different patterns of high b-values exist for different volcanic areas, mainly due to the different magmatic composition and style of eruption.
The forecasting of volcanic eruptions has a distinct advantage to the forecasting of earthquakes.In volcanoes, the magma moves upwards from a depth and this movement can be detected with many different methods (i.e.seismicity, geodesy, electromagnetism, chemistry of fumaroles gases, gases in soils, conductivity).
A combination of many methods might eventually lead to a better understanding of the preparatory process of volcanic activity.In a recent review paper, De Lauro et al. (2011) stated that the Independent Component Analysis (ICA) method of deconvolution of seismic waveforms and the Natural Time Analysis method of analyzing seismic activity (Varotsos et al., 2011) have success in the monitoring and forecasting of volcanic processes in Stromboli (Italy), Erebus (Antarctica), Volcán de Colima (Mexico) and Izu (Japan).
In this study, we investigated the seismicity of the Santorini volcanic island complex using the earthquake catalogue of NOA, in order to contribute to the mosaic of a volcanic and seismic hazard estimates for Santorini.The limited and diffused seismicity of the region does not allow for a more detailed temporal and spatial investigation at this point.However, it is apparent that the detectability and reporting of the NOA seismological network in the Santorini region have recently improved and the seismological data will prove to be valuable for future research studies.

Fig. 2 .
Fig. 2. Seismicity recorded in the Santorini region from 1966 to the end of 2011 (14 December), from the earthquake catalog of the National Observatory of Athens.

Fig. 3 .
Fig. 3. Cumulative seismicity curves of the Santorini region as a function of time.Clustered (blue curve) vs. Declustered (red curve).

Fig. 4 .
Fig. 4. Comparison of the rates as a function of magnitude for the two periods, which are printed at the top.The numbers are normalized by the duration of the periods.(a) Frequency-magnitude curve.(b) Non-cumulative numbers of events as a function of magnitude.(c) Magnitude signature.

Fig. 5 .
Fig. 5. Z-value map showing the rate changes after the recent period of detection increase.(note that negative z-values by definition indicate seismicity rate increases).

Fig. 6 .
Fig. 6.The maximum of the derivative of the frequency magnitude distribution of the declustered catalogue is used to determine the magnitude of completeness (green curve) and the b-value is determined by the weighted least squares method (red curve).

Fig. 7 .
Fig. 7.The temporal variation of the M c in NOA seismicity catalogue for Santorini.

Fig. 8 .
Fig. 8. b-value map showing the spatial variation in the Santorini island volcanic complex from 1966 till 2012, calculated from the declustered seismicity catalogue of NOA.

Fig. 9 .
Fig. 9. b-value cross section showing the spatial variation after the detection increases in 2010 and 2011, along a NE-SW profile across the volcanic complex.