A detailed seismic zonation model for shallow earthquakes in the broader Aegean area

In the present work we propose a new seismic zonation model of area type sources for the broader Aegean area, which can be readily used for seismic hazard assessment. The definition of this model is based not only on seismicity information but incorporates all available seismotectonic and neotectonic information for the study area, in an attempt to define zones which show not only a rather homogeneous seismicity release but also exhibit similar active faulting characteristics. For this reason, all available seismological information such as fault plane solutions and the corresponding kinematic axes have been incorporated in the analysis, as well as information about active tectonics, such as seismic and active faults. Moreover, various morphotectonic features (e.g. relief, coastline) were also considered. Finally, a revised seismic catalogue is employed and earthquake epicentres since historical times (550 BC–2008) are employed, in order to define areas of common seismotectonic characteristics, that could constitute a discrete seismic zone. A new revised model of 113 earthquake seismic zones of shallow earthquakes for the broader Aegean area is finally proposed. Using the proposed zonation model, a detailed study is performed for the catalogue completeness for the recent instrumental period. Using the defined completeness information, seismicity parameters (such as G–R values) for the 113 new seismic zones have been calculated, and their spatial distribution was also examined. The spatial variation of the obtained b values shows an excellent correlation with the geotectonic setting in the area, in good agreement with previous studies. Moreover, a quantitative estimation of seismicity is performed in terms of the mean return period, Tm, of large (M ≥ 6.0) earthquakes, as well as the most frequent maximum magnitude, Mt, for a typical time period (T = 50 yr), revealing significant spatial variations of seismicity levels within the study area. The new proposed seismic zonation model and its parameters can be readily employed for seismic hazard assessment for the broader Aegean area.


Introduction
Even before the first steps of modern Seismology, the reliable assessment of seismic hazard was traditionally considered as one of the most important tasks of Seismology.Especially since the pioneering work of Cornell (1968) it was obvious that this assessment depends on several factors, among which perhaps the most important, and often more poorly understood, is the seismic source model.The compilation of a reliable seismic source model requires the use of high quality (and density) data from various disciplines such as Geology (for the location and configuration of potentially active faults), Geophysics and Geodesy (for the definition of seismic sources not readily observed at the Earth's surface), Active Tectonics (for the rate and type of deformation) and Historical Seismology (for the elaboration of historical records of strong past earthquakes).
Modern studies on probabilistic seismicity and seismic hazard assessment at any scale are usually based on several procedures and algorithms, which typically require the study area to be divided into seismogenic source zones.In several cases, this corresponds to the definition of a set of sources (e.g.sets of polygons for areal seismic sources) over the seismically active region which is examined.In the case of areal sources, the polygons that are often em-Published by Copernicus Publications on behalf of the European Geosciences Union.while the southern Aegean Benioff zone isodepths (white solid lines), the southern Aegean volcanic arc and major basins within the Aegean Sea are also presented (modified from Papazachos et al., 1998;Karagianni et al., 2005;Mountrakis et al., 2012).
ployed may have a complex geometry, reflecting the complexity of major tectonic trends and typically encompass areas within which seismicity may be considered rather homogeneously distributed in space and stationary in time.In most cases, the geometrical definition of the seismic sources is followed by the estimation/adoption of a statistical model for the magnitude distribution of the earthquakes within each source, which is parametrized usually by a Gutenberg-Richter type frequency-magnitude relation (Gutenberg and Richter, 1944).
Published studies on the compilation and adoption of proper seismic sources model include various types of sources, such as point sources (e.g.catalogues of earthquakes), linear sources (e.g.faults) or areal type of sources (polygon-type active regions, with known tectonic and seismic activity parameters).In several cases seismic zonation models have been used for seismicity or earthquake forecasting studies.A first global zonation attempt was performed by Flinn and Engdahl (1965), who proposed a seismic regionalization scheme that was refined later (Flinn et al., 1974) and revised in 1995 (Young et al., 1996).
In the present work we focus our attention in the broader Aegean Sea region (Fig. 1).This area includes the whole Greek area, Albania, FYROM, as well as significant segments of Bulgaria and western Turkey.From a geotectonic point of view the Aegean belongs to the Eurasian conti-nental lithospheric system, at the convergence boundary between Africa and Eurasia.Due to this setting, the Aegean hosts more than 60 % of the total European seismicity.This seismicity is a result of the deformation induced from three main geodynamic sources: (a) the northward motion of the African lithosphere, resulting in the subduction of the Eastern Mediterranean lithosphere under the Aegean (Papazachos and Comninakis, 1970;Le Pichôn and Angelier, 1979), (b) the westward motion of the Anatolia/Turkey microplate along the North Anatolia faults and its continuation in the Northern Aegean (McKenzie, 1970(McKenzie, , 1972)), and (c) the anticlockwise rotation of the Apulia/Adriatic microplate and its collision with Eurasia along the Adriatic, Albania and Epirus (NW Greece) coast (Anderson and Jackson, 1987).These geodynamic phenomena are responsible for the complex seismotectonic setting of this region, which results in the generation of a large number of very strong shallow and intermediate-depth events, with magnitudes up to M = 8.3 and M = 8.0, respectively (Papazachos, 1990).
Several researchers have worked on the seismic zonation problem for the broader Aegean area during the last four decades (Papazachos, 1980(Papazachos, , 1990;;Hatzidimitriou et al., 1985;Papazachos and Papaioannou, 1993;Papaioannou and Papazachos, 2000;Slejko et al., 2010).In a first attempt, Papazachos (1980) used the distribution of earthquake epicentres, the strike and type of every significant known seismic fault, azimuth of P and T stress axes provided by the available focal mechanisms and b parameter values of G-R relation.In the following years, additional information was taken into account by several researchers for seismic source definitions, such as geological data, geomorphological settings and the attenuation pattern observed in macroseismic isoseismals.
More recently, Papaioannou and Papazachos (2000) proposed a model for the south Balkan area containing 67 seismic zones related with shallow earthquakes.At the same time, Papazachos et al. (2001) studied all recent geological and geophysical data in combination with seismological information and determined 159 active faults (rupture zones), related with the known strong (M ≥ 6.0) shallow earthquakes from 550 BC up to 2000 for the broader Aegean area.Based on this seismic source model, Papaioannou (2001) proposed a hybrid model consisting of faults and shallow seismic zones at the same time and used in the compilation of the 2001 revision of the Greek seismic code.Recently, SEAHELLARC Working Group (Slejko et al., 2010) developed a new seismogenic zonation in order to evaluate and compute seismic and tsunamis hazard and risk for SW Peloponnesus (SW Greece).
The seismic zonation models have been incorporated in the seismicity analysis, either in a qualitative or a quantitative manner.The initial attempts were rather qualitative, mostly based on the geographical distribution of earthquake epicentres.For the Aegean area several such studies have been presented for shallow (Galanopoulos, 1963;Makropoulos, 1978;Comninakis and Papazachos, 1978;Papazachos andPapazachou, 1997, 2003), and intermediate-depth earthquakes (Galanopoulos, 1953;Papazachos and Comninakis, 1971;Comninakis and Papazachos, 1980).On the other hand, the quantitative study of seismicity in the Aegean area during the last 50 years was based on the geographical distribution of parameters such as the a and b values of the G-R frequency-magnitude relation, moment rate release, maximum earthquake magnitude (M max ), number of events with a magnitude larger than a specific magnitude, mean return period (T m ) of strong earthquakes, most probable magnitude in a certain time period (M t ), or the probability of occurrence of earthquakes with magnitudes over a certain value (Delibasis and Galanopoulos, 1965;Comninakis, 1975;Makropoulos and Burton, 1985;Stavrakakis and Tselentis, 1987;Papazachos, 1990Papazachos, , 1999;;Papadopoulos and Kijko, 1991;Hatzidimitriou et al., 1994;Tsapanos et al., 2003, among others).
It is well known that both the aleatory and epistemic uncertainties influence the results of any hazard analysis.In order to reduce the epistemic uncertainties related to the geometry of the seismogenic sources and their seismicity parameters, an attempt is made in the present paper to propose a new model of seismic sources for the broader Aegean, merging the knowledge accumulated during the last 30 years on active tectonics, seismotectonics, seismicity, and so forth, using a high-quality data set.The starting point for this effort is based on a model of active tectonics in the broader area (Papazachos et al., 1998), which defines five generalized zones on the basis of the dominant faulting type, which are: (1) the reverse faulting/compressional zone across the Adriatic Sea and the North Ionian Sea (e.g.Anderson and Jackson, 1987), (2) the continuation of this reverse faulting zone along the southern Aegean Hellenic Arc (e.g.Papazachos and Delibasis, 1969), (3) the extension faulting zone across the Hellenides mountain range (e.g.Papazachos et al., 1984), (4) the extension faulting zone in the back-arc region with general N-S direction, located in the inner part of the arc, that geographically covers the area from S. Bulgaria and FYROM, northern and central Greece and follows the curvature of the volcanic arc in the South Aegean to southwestern Turkey (e.g.McKenzie, 1970McKenzie, , 1972)), and (5) the strike-slip faulting zone starting from the western part of the Anatolian fault and its continuation in the North Aegean Trough (e.g.Taymaz et al., 1991) and the corresponding strike-slip zone in the central Ionian islands (Scordilis et al., 1985;Papazachos et al., 1994).
The geographical distribution and density of the seismological stations in the broader Aegean until quite recently (last decade) did not allow for detection of small enough hypocentral uncertainties, in order to correlate earthquakes with specific known faults and proceed in the determination of fault-type sources and the corresponding seismicity parameters.For this reason.we have based the new zonation model on areal-type polygonal sources; however, these zones are strongly associated with the seismotectonic features of the study area.More specifically, the proposed zonation model includes 113 new shallow seismic sources which have been defined on the basis of focal mechanisms, P and T kinematic axes, seismicity, as well as information on seismic and active geological faults.Furthermore, a revised homogeneous earthquake catalogue for the period 550 BC-2008 is employed for quantitative estimations.Using the revised earthquake catalogue, a detailed study is made for magnitude completeness since 1981, applying the ZMAP software (Wiemer, 2001), while for older periods the results from Papazachos and Papazachou (2003) are adopted.Using the complete data set, we compute the typical Gutenberg-Richter magnitude-frequency relation parameters for the 113 new shallow sources and present the spatial distribution of several seismicity parameters.The results illustrate areas of high seismicity (central Ionian islands), as well as lowseismicity regions (part of FYROM, central part of the southern Aegean Sea, Thrace and Bulgaria).The obtained new model of seismic sources and its parameters can be readily employed for reliable seismic hazard assessment studies at various scales.

Generation of a geological-geophysical database for
the new zonation model generation

Fault plane solution information
Fault plane solutions constitute a very useful tool for the study of the source characteristics of an earthquake.Focal mechanisms illustrate the connection between earthquakes and specific known faults, the general and local stress field, as well as their relation with other possible faults not already known, whose existence we can often only hypothesize.The contribution of fault plane solutions in the progress of seismotectonic understanding in the broader Aegean/southern Balkan area is very important, as they historically allowed the determination of the reverse faulting zones of the convex part of the Hellenic Arc and Albanian coast (Papazachos and Delibasis, 1969), of the extension field in continental Greece and surrounding southern Balkan area (McKenzie, 1970(McKenzie, , 1972) ) and its continuation to the north of the North Aegean (Papazachos et al., 1979;Taymaz et al., 1991).Dextral strike-slip faults in the Cephalonia area (Scordilis et al., 1985;Louvari et al., 1999) and the North Aegean Trough (Galanopoulos, 1967;McKenzie, 1972) are some other examples of the important contribution of fault plane solutions.
For the new proposed seismic zonation model, fault plane solutions concerning 177 shallow earthquakes with magnitudes M ≥ 5.5 (Papazachos and Papazachou, 2003) have been initially collected for the period 1953-1999.In order to improve this data set, we decided to enhance it with more recent and smaller events (M < 5.5).For this reason, fault plane solutions were collected from several Balkan and international seismological centres.The Global Centroid-Moment-Tensor (Global CMT) project, European-Mediterranean Seismological Centre (EMSC-CSEM) and International Seismological Centre (ISC) are the main sources of this data type, as they collect and publish information from numerous international (e.g.ETH, INGV) and Balkan (AUTH, NOA, UPSL, KOERI) seismological centres.Using these sources, we collected 364 fault plane solutions determined from waveform modelling, mainly concerning earthquakes with magnitude M ≥ 4.5 for the period 1978-2009.
Additional information was incorporated from Benetatos ( 2007), who has compiled a similar catalogue from alternative sources and the Global CMT project and collected information about fault plane solutions of 33 strong earthquakes (M 5:0) for the years 1980-2005.Moreover, a significant number of additional fault plane solutions (184 events) were calculated by Mountrakis et al. (2006), who studied the active tectonics and stress field in northern Greece, using waveform data from the AUTH Geophysical Laboratory seismological network (P wave first motions for events with magnitudes 3.0 ≤ M ≤ 5.8), with epicentres located in northern Greece and southern Albania, FYROM and Bulgaria for the years 1989-1999.In addition, nine focal mechanisms derived by the analysis of reliable information about seismic faults with major surface ruptures, which were associated with strong earthquakes (M ≥ 6.7) in the years 1861-1953 (Papazachos and Papazachou, 2003), were also added to the database.Following this procedure, a new database with fault plane solution information of 767 shallow earthquakes was created, containing events mostly from the period 1953-2009, strongly weighted towards the last 10 years.
In Fig. 2, the collected fault plane solutions are presented with different colours according to the rupture type related with every event -green for normal faults (a), brown for reverse faults (b) and blue for strike-slip faults (c).Information provided by those fault plane solutions is of special importance, as we can easily recognize regions with similar earthquake faulting characteristics, suggesting a specific seismotectonic setting and a dominant stress regime.More specifically, the fault plane solution distribution of Fig. 2 confirms the general pattern of reverse ruptures along the coastline of Montenegro and Albania, in the Ionian Sea and the outer part of subduction zone of the Hellenic Arc.Strike-slip focal mechanisms are mainly recognized in the North Aegean, the Marmara Sea and the northwestern part of Turkey, following the North Anatolia fault branches and its continuation in the North Aegean Trough, as well as in the central Ionian Sea and a significant part of NW Peloponnesus.Finally, normal faults associated with corresponding focal mechanisms dominate the largest part of continental Greece and SW Turkey.This pattern agrees with current knowledge on the seismotectonic model and the geodynamics in the broader area of the Aegean Sea (e.g.Papazachos et al., 1998).We note that in a few cases where the two nodal planes show a different strike-slip and dip-slip (thrust or normal) character, the focal mechanism may appear as a strike-slip type.However, it is clear that these mechanisms (e.g. three thrust/strike-slip events in the coast of Albania and one normal/strike-slip in inner-northern Albania) are in excellent agreement with the corresponding spatial distribution of pure thrust or normal faults.
Exploring the fault plane solution database, several useful information can be derived from the local stress field, by analysing the corresponding P and T principal kinematic axes.According to the rupture type related with each fault plane solution (normal, reverse, strike-slip), principal kinematic axes exhibit a significant horizontal or vertical component.The azimuthal distribution of these axes can help to divide the earthquakes into groups of homogeneous rupture type, which provides an important criterion for the division of the study area in smaller zones of common seismotectonic characteristics.The spatial distribution of principal compression P axes and extension T axes is illustrated in Fig. 3a and  b, with red and green arrows, respectively.In both figures only sub-horizontal axes, e.g. with a dip angle less than 45 • , are presented.
The geographical distribution of P axes (Fig. 3a) suggests that the compression stress field is located exclusively in the convergence zone between Mediterranean and Eurasian tectonic plate along the outer Hellenic Arc (see also Fig. 1) and its continuation north of the Cephalonia fault in NW Greece and coastal Albania, with a general orientation normal to the arc (yellow-shadowed polygon).A second area includes the branches of the strike-slip North Anatolia fault-North Aegean Trough (North Aegean-NW Turkey), with a general NW-SE direction (red-shadowed polygon).In the remaining part of the study area, sub-horizontal compressive axes are almost completely absent.
Concerning the T axes distribution (Fig. 3b), it is clear that the dominant stress field is extensional for the largest of the study area.Geographically, this extension appears to have a general N-S to NNE-SSW direction in the eastern part of the study area (West Turkey, North Aegean), which is locally modified to a NNW-SSE to NW-SE direction, as we move to the west towards the continental part of Greece, in agreement with earlier neotectonic (e.g.Mercier et al., 1989) and seismological (Papazachos and Kiratzi, 1996) observations.The direction of extension changes completely in the area adjacent to the Hellenic Arc subduction zone and the Adriatic-Eurasia southern collision system, since the T axes exhibit a more or less E-W direction across the Albanides and Hellenides mountain chain, up to the central part on southern Peloponnesus and the area of Cythera-Crete-Carpathos-Rhodes (South Greece).

Information on active faults and seismicity
Major active faults and rupture zones were also considered for the separation of the study area in smaller, seismotectonically coherent areas, in order to create a new seismic zonation model.Papazachos et al. (2001) published a list of 159 major faults related with the strong known earthquakes in Greece since the 6th century BC.These faults constituted the base for the new seismic zonation model and are depicted in black in Fig. 4. Additional information from Mountrakis et al. (2010) concerning active neotectonic faults in the broader area of Greece, derived from a large number of neotectonic studies, were also added in the fault database.These faults are also presented in Fig. 4, with a different colour according to the proposed fault type.
As expected, most of the faults are identified by both research works, hence most of the linear faults observed in Fig. 4 are practically identical.However, local differences are observed between the two databases, partly reflecting the different procedure that was adopted for the fault determination.In the first case, Papazachos et al. (2001) collected information about active "seismic" faults, determined on the basis of mainly seismological and to a lesser extent geologicalmorphotectonic data, associated mostly with strong historical earthquakes.Fault dimensions were mainly determined on the basis of the corresponding maximum earthquake magnitude observed.On the other hand, the Mountrakis et al. (2010) database includes information about all active, already recognized and mapped neotectonic faults, which exhibited a significant surface fault trace.
Epicentres from all recorded earthquakes from 550 BC up to 2008 from the catalogue of Papazachos et al. ( 2010) (http: //geophysics.geo.auth.gr/ss/CATALOGS/seiscat.dat) consist of the most important ingredient, in order to determine areas exhibiting of a common seismotectonic/seismicity status.The magnitudes in this catalogue have been revised by Vamvakaris (2010), who employed a large number of magnitude scaling relations to convert local and other magnitudes from a large number of agencies to equivalent moment magnitudes estimates.In general, the results of Vamvakaris (2010) showed that new equivalent moment magnitudes were slightly smaller compared to the original Papazachos et al. ( 2010) estimates for M < 4.5, while showing a Historical and recent earthquakes have a critical contribution to the new seismic zonation model, as they express the final "result" of the active tectonic deformation in the study region.Of course, more recent, instrumentally recorded data are of higher accuracy, hence they participate with a different weighting in the zonation procedure.In Fig. 5 epicentres of all shallow earthquake of the catalogue used with magnitude M ≥ 4.0 are presented.Solid circles of different colours depict epicentres for the instrumental period 1911-2008, while open circles represent epicentres of historical earthquakes (before 1911, when the first seismograph was installed in Athens).From the seismicity distribution shown in this figure, several areas with high and low epicentre density can be recognized, illustrating high and low seismicity level regions, respectively.Additional information for the seismicity level features of the broader Aegean/southern Balkan area is given in the following section, where the new seismic zonation model is presented in detail.

Development of the new seismic zonation model
As mentioned previously, the separation of the study area into smaller, seismotectonically homogeneous zones is based on criteria mainly related with the epicentre distribution and the types of faults identified in this area.Moreover, additional information from the stress field distribution and other geological, geotectonic and morphotectonic characteristics has been also taken into account.In the present work, we decided to employ simple geometric shapes for the definition of the seismic source model, in compliance with the criteria described previously.In general, the combination of all available information can help to determine areas of common seismotectonic characteristics.The type and strike of major faults and the stress regime were primarily considered for the definition of the general shape and grouping of the proposed zones, while the epicentre distribution was considered for the detailed assessment of the transition borders between neighbouring seismic zones.
In Fig. 6 we present an example of zonation for the northwestern part of the study area, for the region containing S. Montenegro, Albania, Greece (Epirus and W. Macedonia) and West FYROM.In this map, seismic faults (Papazachos et al., 2001) are depicted with black and different patterns, according to their rupture type, while neotectonic faults (Mountrakis et al., 2010) are shown with different colours (white for normal, black for reverse and red for strike-slip faults).Moreover, all available information regarding principal P (red arrows) and T (black arrows) kinematic axes derived from fault plane solutions is also shown in the same figure.Finally, the epicentres of all earthquakes before and after 1911 are presented with open and solid circles, respectively (4.0 ≤ M < 6.0, white circles; 6.0 ≤ M < 7.0, pink circles; and M > 7.0, red circles).
In the northwestern part of that area we can identify several large reverse faults following the coastline of Montenegro, with a general NW-SE direction.These faults continue in West Albania with branches of NNW-SSE to N-Strending faults, and return again to NW-SE directions as they follow the coastline of SW Albania and NW Greece (Epirus).The dipping direction of these faults is towards the east and there is an obvious connection with the compressive stress field detected in the Adriatic Sea (e.g.Anderson and Jackson, 1987), which reaches the central Ionian Sea and is interrupted by the dextral strike-slip fault of Cephalonia (Scordilis et al., 1985;Papazachos et al., 1994).All these common characteristics suggest that a group of seismic zones of reverse ruptures should be considered, that could be separated into four more detailed seismic zones, based mainly on their geographic location, epicentre distribution, local seismotectonic characteristics (e.g.secondary ruptures), as well as geomorphological features, such as relief variations (e.g.Quaternary basins).These four seismic zones (T-A1, T-A2, T-A3 and T-A4) are clearly associated with the compressive stress field across the West Albania/NW Greece coast.To the west of the two southern zones there is a sparser concentration of earthquake epicentres.For this reason, an additional seismic zone (T-A5), belonging to the same group of reverse faulting zones (T-group), was defined, due to its distance from the main faulting branches, and its lower seismicity level.The southern limit of zone T-A4 was defined in the area where a major change of the stress field was observed.This change is a result of the different seismotectonic settings described by the active faulting, the fault plane solutions, the major stress axes and the earthquake epicentres in the Amvrakikos Bay and the northern part of Lefkada Island.
East of this group of zones which are controlled by a compressive stress regime, across the Albanides and Hellenides mountains, we can recognize a narrow region with significant normal faults of a general N-S direction.These faults were created by an extensional E-W stress field (Papazachos et al., 1984;Kiratzi et al., 1987;Armijo et al., 1992).Major local neotectonic features, like the Ohrida rupture zone (associated with strong earthquakes in the past) and the small but significant variation in fault orientations, as well as the distribution of local earthquake epicentres, have led us to the definition of four smaller seismic zones (N-B1, N-B2, N-B3 and N-B4) along this region of E-W extension.These zones follow, from north to south, the major active faults and the local seismicity, constituting a transition zone, separating the strong ENE-WSW compressive stresses of the Adriatic-Ionian coast with the N-S extensional stress field which dominates in the largest part of North Greece, FY-ROM and Bulgaria, where we find significant normal faults of an E-W to ENE-WSW general direction.
In the map of Fig. 6 a large NE-SW normal fault dominates the area of Skopje, associated with strong historical as well as recent seismicity.This area was defined as a separate zone (N-E1) and is clearly separated in the south by a different zone, where we can find the N-S-trending fault of Bitola and the local concentration of earthquake epicentres around Prespes lakes, in the border region between Greece, Albania and FYROM.This zone (N-E3) is interrupted to the east, where a significant drop of seismicity levels is observed for a large area.The lack of earthquake epicentres and the corresponding absence of major seismotectonic features (e.g.active faults) have led us to the definition of this area as a separate zone, characterized by very low seismicity levels (N-E4).
To the south, in the West Macedonia region, active rupture zones are again modified.In this area we can find multiple rupture zones of NE-SW direction associated with strong historical earthquakes, mostly in the broader area of the city of Kastoria.The boundaries of zone N-E11 have been determined by the general spatial pattern of seismicity, which is mainly located in the zone's southern section, due to the recent (1995) Kozani-Grevena earthquake sequence.The broader area associated with this specific rupture forms zone N-E14, which is clearly separated from the normal faults of North Thessaly.This North Thessaly faulting system exhibits a different fault strike direction (generally E-W), which justifies the separate definition of the corresponding seismic zones, taking also into account the significant seismicity in the area of the cities of Trikala and Karditsa during historical times.
Following a similar approach, the broader area of the Aegean has been separated into 113 seismic zones with common seismotectonic features.The new seismic zonation map is presented in Fig. 7, where all available data employed for the analysis are also superimposed.Nevertheless, the 113 new seismic zones exhibit several common characteristics, hence it is possible to provide a general grouping of the in-

Detailed description of the new proposed seismic zonation model
The T-A group of zones geographically corresponds to the area of Adriatic-Albania/W.Greece convergence, where reverse ruptures with a general NW-SE strike dominate, due by the convergence compressive stress field.This stress field is expressed with sub-horizontal P axes that, in general, follow a direction normal to the coastline of Montenegro, Albania and NW Greece (Epirus) along the Adriatic and Ionian Sea.This group includes five seismic zones, mainly divided on the basis of their geographical location, fault direction and seismicity level, as previously explained in detail.Zone T-A5 is slightly different than the other zones (Figs. 6 and 9) and represents the outer part of zones T-A3 and T-A4, with a lower seismicity level, without important seismotectonic characteristics.
A narrow area of normal faults with a NNW-SSE to N-S direction is located in the eastern part of Albanides and Pindus mountains.This faulting pattern modification is in agreement with the stress field change in this area, recognized by fault plane solutions and the corresponding T axes direction (almost E-W).This transition between the NW-SE thrust faults in the west and the normal ∼ E-W-trending faults commonly found in FYROM, Bulgaria and Macedonia, defines the N-B group which has been separated in four seismic zones, as previously described in detail.
The S-C group includes the northern central Ionian islands (Cephalonia and Lefkada) and exhibits the highest seismic- ity levels in the whole study area.Large dextral NE-SW to NNE-SSW strike-slip faults are found in this area, interrupting the reverse faulting pattern transition from the southern Adriatic and NW Greek coast to the South Ionian Sea and the western section of the Hellenic Arc.The division of this group in five seismic zones was mainly based on the seismotectonic setting details (faults, focal mechanisms, earthquake epicentres, etc.) as we move from the western to the eastern part of these islands.
The T-D group consists of zones characterized by significant thrust faults and strong earthquakes related with corresponding seismotectonic activity.The roughly NW-SE strike of reverse faults remains rather stable and does not follow the Hellenic Arc trend.Along the Hellenic Arc we have defined 13 seismic zones with common seismotectonic characteristics.An exception of this rule is the area of Rhodes-Karpathos (zones T-D10 and T-D12), where a large sinistral strike-slip fault is also observed, associated with several historical and recent strong earthquakes (e.g.1957 M = 7.2 event and its M = 6.8 preshock).In the outer part of the Hel-lenic Arc another group of lower seismicity level zones, following the inner ones, has been identified.
In a wide area from FYROM to the Black Sea and from central Bulgaria to the North Aegean Sea, normal faults with a general E-W to NE-SW direction are found.Different local seismotectonic settings have been used to distinguish 26 individual seismic zones in this N-E group of zones.Among these zones, there are areas where seismicity spatially clusters around specific fault zones (e.g.zones N-E5, N-E6, N-E9, N-E14, N-E16 and N-E17), as well as areas where a relative seismicity gap is observed in both instrumental and historical times (e.g.zones NE4, N-E8, N-E18 and E24).
The N-F group of zones includes all normal faults of Central Greece, the Corinth Gulf and Attica.Seismicity is mainly located along the S. Thessaly fault zone, as well as Maliakos and Corinth gulfs.A rather uniform extension stress field of ∼ N-S direction is recognized for the whole area.The spatial distribution of seismicity data and other local seismotectonic characteristics, such as morphology and active faults, have been employed for the division of this group into 13 individual seismic zones.
Two seismic zones in the northwestern part of Peloponnesus form the N-H group.The main reason for the separation of these two zones from the remaining North Peloponnesus region is the presence of characteristic dextral strikeslip faults of NE-SW direction in Ilia and Kyparissiakos Gulf, similar to the neighbouring S-C group.This faulting pattern (e.g. the recent Andravida-Patras 2008 event) overlaps and modifies the general trend of E-W normal faulting observed to the north and the N-S normal faulting observed further to the south (Kalamata area).
The N-L group of zones covers the remaining part of Peloponnesus, Crete and Karpathos Island.In the nine individual seismic zones of this group medium to large N-S to NNW-SSE normal active faults are observed, resulting in the generation of several strong historical and recent shallow earthquakes.Seismicity levels, focal mechanisms and fault parameters were the main criteria for the zone separation.Zones N-L1, N-L3, N-L8 and N-L9 exhibit high seismicity levels, while zones N-L4, N-L5, N-L6 and N-L7 were included in the same group due to their location and the similar kind of fault ruptures observed.
In the S-I group, significant dextral strike-slip faulting is observed, related to the horizontal displacement of the Anatolia microplate parallel to the Anatolian fault, towards the west, in the northern Aegean area (e.g.Taymaz et al., 1991).15 seismic zones have been defined for this region, with 10 of them located in the North Aegean Trough and the islands of the NE Aegean Sea.The remaining five seismic zones include the area around the Marmara Sea, from Gulf of Saros to Istanbul and from East Thrace to Bursa.Dextral ruptures depicted from the corresponding focal mechanisms and stress field, as well as high seismicity levels with strong earthquakes, are the main features of this area, with the NE-SWtrending North Aegean Trough fault branches acting along the North Aegean area up to North Sporades islands.Several ruptures with similar strike are detected to the south of the trough.On the other hand, the dominant fault orientations in the Marmara Sea have a general E-W trend, though earthquakes that have occurred during the last century are generally few and of relatively small magnitude.
In the Central Aegean Sea the significant lack of strong earthquakes and large active faults resulted in the definition of four individual seismic zones, forming the N-J group.Due to its significant local seismic activity, partly seismovolcanic, the Milos island is considered as a separate, individual seismic zone (N-J4).
The last N-K group contains 17 seismic zones mainly located in W. Turkey and the eastern Greek Aegean islands (e.g.Chios, Samos, Kos, Amorgos).In general, normal faults are observed in this group, with some local exceptions where strike-slip faults are also detected (e.g.dextral faults in N-K2 zone, sinistral faults in N-K10).The general fault strike is roughly E-W for zones along the Aegean coastline (N-K1, N-K2, N-K7, N-K8 and N-K14 zones), although ENE-WSW (N-K3, N-K11 and N-K12 zones) and WNW-ESEtrending branches (N-K5, N-K6 and N-K9 zones) are also found.A special seismic zone in this group is the Santorini Island (NK11), where the local seismicity is mainly related with the volcanic activity that differs from the other zones of N-J group, though it exhibits several common geotectonic characteristics with the neighbouring zones of Amorgos and Kos (N-K12 and N-K13, respectively).
The definition of the geometrical features of the seismic zonation model was followed by an effort to quantitatively assess the seismicity level for each seismic zone.This assessment requires the study of the corresponding seismicity catalogue completeness (with respect to magnitudes), as well as the determination of several quantitative seismicity parameters, as presented in detail in the following.

Data completeness
The quantitative assessment of seismicity requires the knowledge of the completeness of the employed seismicity catalogue.The study area shows extremely variable magnitude completeness, both spatially and temporally (e.g.D' Alessandro et al., 2011), reflecting the corresponding changes of the historical information sources, as well as of the network coverage and density.Of course the network configuration simultaneously affects not only the data completeness, but also the accuracy of the corresponding seismological information (e.g.hypocentral errors, etc.).Moreover, the data completeness critically affects several seismological scientific targets, such as earthquake prediction, seismic hazard, etc., as it controls the completeness and accuracy of the seismicity parameter information determined from catalogue data (e.g.Gutenberg-Richter, G-R, parameters, large earthquakes associated with specific active faults, etc.).In the present work, we adopted the catalogue completeness for M ≥ 4.5 and data before 1981, which has been extensively studied by Papazachos et al. (2000) and was published as part of the corresponding earthquake catalogue.This catalogue, which is included in the catalogue of the present study, contains complete data for shallow earthquakes for the following periods and magnitudes; M ≥ 4.5 since 1970, M ≥ 5.0 since 1950, M ≥ 5.2 since 1911, M ≥ 6.5 since 1845 and M ≥ 7.2 since 1500 (Papazachos and Papazachou, 2003).A detailed study was performed in the present work for smaller earthquakes for the period 1981-2008, when the monitoring network in the broader Aegean region was gradually densified and a much larger number of smaller magnitude events was recorded throughout the entire study area.
Since a large number of small-magnitude events was available for the analysis, we performed the completeness study separately for each one of the 113 new seismic zones proposed in the present work.Several methods have been proposed for the determination of catalogue completeness (Wyss et al., 1999;Wiemer and Wyss, 2000;Woessner and Wiemer, 2005), essentially relying on the estimation of the completeness cut-off from the G-R relation as the magnitude, M C , at which the cumulative frequency-magnitude distribution stops following the linear equation (Zuniga and Wyss, 1995).In the present work, completeness analysis was performed by applying the ZMAP software (Wiemer, 2001), based on the techniques of Maximum Curvature (Wyss et al., 1999), and the goodness-of-fit test (level of fit, 90 and 95 %).For each case ZMAP selects the best combination of these techniques resulting in more robust values, on the presumption that there are sufficient data for the study area.In the cases of fewer data, where ZMAP could not produce reliable results, we performed an estimation from the simultaneous visual inspection of the number, n, of earthquakes with magnitude M, combined with the logarithm of the cumulative number, log N, of earthquakes with magnitude ≥ M for different time periods.This approach essentially identifies the break of the log N −M curve, similar to the Maximum Curvature method, while observing the exponential variation of the n(M) curve and the good linear fit of the logN − M variation (similar to the approach of Wiemer and Wyss, 2000).Following this approach, the completeness was defined for each zone since 1981 and in a few cases additional completeness level were defined for later intervals.In most cases the completeness magnitude remained rather constant for each zone, gradually improving for recent time periods.
Magnitude completeness information for each zone since 1981 is given in Table 1.Wherever additional completeness results (lower completeness magnitude) were available, such information is also included in Table 1.A typical result is presented in Fig. 9, where the corresponding spatial variation of the completeness magnitudes for years 1981 and 1990 is presented.Both plots are created by assigning the corresponding completeness magnitude at the centre of each zone and using a simple triangulation interpolation, without any spatial smoothing.From Fig. 9 it is clear that despite the fact that the completeness magnitude, M C , was defined using a combination of analytical methods through ZMAP and hand-picking for zones with fewer data separately, the obtained results show a strong spatial coherence, with smooth variations, indicative of the robustness of the M C determination procedure.
In general, useful information can be determined from the geographical distribution of completeness magnitude.The M C variation for the time period  shows that there are areas with relatively large completeness magnitudes in the range of 4.0-4.5 (e.g.south of Crete and Rhodes), whereas in other areas where we can find complete data for much smaller magnitudes, in the range of 2.3-2.7 (e.g. for the broader Thessaloniki, Athens and Marmara Sea region).
In general, we have a complete catalogue after 1981 for magnitudes M ≥ 3.2 for the largest part of the continental section of the study area.This observation is clearly justified by the geographical distribution of seismological networks operating in Greece and surrounding areas.Moreover, it is expected to have improved detectability around the major metropolitan centres such as Thessaloniki, Athens and Istanbul, where civil and financial reasons led to the installation of dense monitoring networks of seismographs.On the other hand, the poor coverage and azimuthal distribution of the network in some areas (e.g. the outer Hellenic Arc, Libyan Sea) cannot provide data with completeness in magnitudes less than 4.0, due to the increased uncertainties of the data and the reduced detectability of the networks.Although the completeness level has improved for recent years (see completeness after 1990), especially for some areas (e.g.Ionian islands, North Aegean Trough, outer Hellenic Arc), the main features of the spatial distribution of the completeness magnitude remain the same.

Seismicity parameter assessment for the new zonation model
The final determination of the G-R seismicity parameters is not a straightforward task.Although these parameters are often associated with the geotectonic characteristics of the studied area (e.g.stress regime), the determined values, and especially the G-R slope b, are very sensitive to the quality and size variability of the data set used for the calculations (Papazachos, 1974).For the estimation of the a and b Gutenberg-Richter parameters we employed the revised earthquake catalogue for the period 550 BC-2008 and the completeness estimations for each one of the 113 shallow seismic zones, as described before.The magnitudefrequency distribution of each seismic zone has been studied using different completeness intervals, and we employed the mean value method (Milne and Davenport, 1969) to reduce data of different completeness periods to a common time window. .Moreover, we also applied the maximum likelihood method (MLM) (Weichert, 1980) in order to calculate independently the b values for each one of the 113 shallow seismic zones.In order to examine the origin of the b value variability, in both cases, we studied the effect of the employed magnitude range of the complete data set for each zone (dM = M max − M min ) on the G-R estimation.It is well known that data with small dM values often lead to unreliable estimates of the b value (e.g.Papazachos, 1974).This effect is clearly seen in Fig. 10, where obtained b values are plotted against dM for all 113 seismic zones.It is clear that for small dM, unrealistic (mainly high) b values are obtained, as a result of the small magnitude range available for these seismic zones.Comparing the results from both least-squares and MLM (fig.10) we note that b values derived from MLM show a much wider distribution not only in cases with small dM (< 1.9).Moreover, even cases of large dM which should result in more stable b values exhibit a broad value distribution.This comparison suggests that the employed least-squares method provides more stable results for enough large dM range, with b values concentrating in a narrow area around the value −1.0, contrary to the MLM results which show a much wider spread of values is shown.For this reason we decided to use the results from the least-squares method and more specifically only cases with a larger magnitude range (dM ≥ 1.9), which excluded less than 12 % of the b values considered for the final calcula- tions.In order to cover the resulting spatial gaps, the selected b values were spatially interpolated using an adapted linear variogram and a kriging gridding method.From this spatial interpolation, final values for the seismic parameter b were re-calculated for each one of the 113 shallow seismic zones.
The obtained b values from the previous procedure are presented in the contour map of Fig. 11.The observed spatial variation is in very good agreement with older studies (Hatzidimitriou et al., 1985(Hatzidimitriou et al., , 1994;;Papazachos, 1999), despite the different methods and data sets used.Concerning the northern part of the study area, b values seem to be strongly related with the seismotectonic settings (Papazachos, 1999) and the distribution of geological zones (Hatzidimitriou et al., 1985(Hatzidimitriou et al., , 1994)).These values are gradu- Using the distribution of b values presented in Fig. 11, we have calculated the corresponding a values from the complete seismicity data.Moreover, we reduced the estimated a values which refer to the reduced time period, t, of complete data and the seismic zone surface, S, to the corresponding annual values, a 1 , for a typical area of 10 000 km 2 (a 1 = alog(tS)+4).The corresponding values are given in Table 1.In order to examine the spatial variation of seismicity in the study region, we employed several typical seismicity measures, such as the mean return period, T m , of earthquakes with magnitude greater or equal to a specific magnitude, M: or by the most probable maximum magnitude, M t , observed for a specific time period, T : (2) All these calculations have been performed with the a 1 parameter (reduced to a time period of 1 year and an area of 10 000 km 2 ) in order to reflect the spatial variation of the corresponding measures, without being influenced by the zone dimensions or completeness interval.In Figs. 12 and 13 we present the spatial distribution of the mean return period for M ≥ 6.0 (a typical damaging earthquake, occurring almost everywhere throughout the study area) and most probable maximum magnitude for a time period of 50 years (typical infrastructure life/design period).Both plots provide interesting information on the seismicity level variations within the study region.In Fig. 12 we can recognize areas of high seismicity, identified by the low mean return period values for shallow earthquakes with M ≥ 6.0.More specifically, in the Central Ionian Sea (islands of Lefkada, Cephalonia and Zakynthos) earthquakes of magnitude greater than 6.0 for an area of 100 km × 100 km occur approximately every 5 years or even less.Values of the order of 10-20 years are observed for the remaining parts of the Ionian Sea, the northwestern Peloponnesus, the Corinth and Patras gulfs and the North Sporades islands in the northwestern Aegean.Return periods of the order of 20-50 years are recognized for the outer Hellenic Arc, the Adriatic-Albania/NW Greece collision zone, the eastern Aegean Sea and Asia Minor coast, including the Amorgos-Santorini zone, as well as southern Thessaly and the Serbo-Macedonian massif in northern Greece.On the other hand, very long return period values (> 200 years) are found for the broader Cyclades Plateau, the western Cretan Basin, the area north of the North Aegean Trough (broader Thasos Island area) and large parts of Thrace and FYROM.These results are in good agreement with the study of Papazachos (1999), despite the fact that he used fewer data, with different characteristics (completeness, accuracy, etc.) and a different method that involved no regional zonation.Similar conclusions can be derived from the most probable maximum magnitudes for a return period of 50 yr.High seismicity areas are also shown in this map, such as the Central Ionian Sea, where shallow earthquakes of maximum magnitude M > 6.6 (locally reaching 7.0) are expected to occur every roughly 50 years.On the contrary, for low seismicity areas (large parts of the Cyclades Island plateau, Cretan Basin, broader Thasos area, western Thrace, eastern Macedonia, southern Bulgaria and FYROM), the statistically most probable larger shallow earthquake for a 50-year period (and an area of 10 000 km 2 ) is expected to be smaller than M = 5.0.

Conclusions
We present a new seismic zonation model for the broader Aegean area, using all available seismotectonic information such as focal mechanisms, P and T kinematic axes, seismic and neotectonic faults, as well as a revised earthquake catalogue for the period 550 BC-2008.The proposed seismic zonation model consists of 113 shallow seismic zones for the broader Aegean/southern Balkan area, specified on the basis of mainly seismotectonic and seismicity criteria.11 groups of zones and three clusters of groups were identified, on the basis of their common seismotectonic features and the location of each zone.These three clusters of zone groups are associated with predominantly normal, thrust and strike-slip ruptures occurring within the 113 shallow seismic zones.
Using the revised earthquake catalogue, a detailed study was performed for the completeness magnitude after 1981, while older completeness results (Papazachos and Papazachou, 2003) were used for the period before 1981.Working on complete data, a and b seismicity parameters of the Gutenberg-Richter magnitude-frequency relation were estimated for each one of the 113 shallow zones.The spatial distribution of b values shows a strong relation with the seismotectonic settings and the geological units of the study area, in good agreement with previous studies (Hatzidimitriou et al., 1985(Hatzidimitriou et al., , 1994;;Papazachos, 1999), with a spatial variation from −1.2 to −0.85, as we move from the southern Hellenic Arc to the northwestern part of Greece and southern Bulgaria.
The calculated a and b seismicity parameters of each zone were used for the detailed quantitative estimation of seismicity in the study area.The mean return period, T m , of earthquakes with magnitude greater than or equal to a specific magnitude and the most probable maximum magnitude, M t , for a specific time period, were calculated for M ≥ 6.0 and T = 50 years, respectively.The spatial distribution of both T M>6.0 and M 50 yr , illustrate large contrasts between areas of high and low seismicity.For example in the central Ionian islands (Cephalonia, Zakynthos, Lefkada) strong earthquakes of magnitude greater than 6.0 occur more often than every 5 years for a typical area of 100 km × 100 km.NW Peloponnesus, the Corinth-Patras Gulf and the North Aegean Trough also exhibit very high levels of shallow seismicity.On the other hand, the largest part of FYROM, sections of northern, central and southwestern Aegean Sea (e.g.north of the North Aegean Trough, western Cretan Basin), as well as Thrace and southern Bulgaria, show much lower shallow seismicity levels, with return periods exceeding 200 years for a typical strong earthquake with M6.0 occurring in the same area (100 km × 100 km).These observed strong spatial seismicity variations, as well as the associated seismotectonic information inherently contained in the definition of each seismic zone, suggest that the new proposed seismic zonation model and its quantitative parameters are an important tool that can be efficiently used for seismic hazard assessment of the broader Aegean region.

Figure 2 .
Figure 2. Fault plane solutions for major (M ≥ 4.0) shallow earthquakes (a: normal, b: thrust and c: strike-slip) collected and used for the new seismic zonation model determination.

Figure 4 .
Figure 4. Map of the 159 major active faults (in black) for the study area(Papazachos et al., 2001) associated with large historical and instrumental earthquakes.Major neotectonic faults(Mountrakis et al., 2010) are also presented with different colours, according their rupture type (yellow for reverse faults, green for normal, and red for strike-slip faults).

Figure 5 .
Figure 5. Map of the epicentres of all known shallow (h < 50 km) earthquakes (M ≥ 4.0) which occurred in the broader Aegean/southern Balkan area during the time period 550 BC-2008.The colour and type of circles denote the magnitude scale and the time period, respectively.

Figure 6 .
Figure 6.A typical example of seismic zone separation for the region of Albania, Northwest Greece, South Montenegro and West FYROM.All available data of seismic and geological faults, stress axes, and earthquake epicentres are combined into a single map.Yellow frames depict the new seismic zone boundaries derived from this area.

Figure 7 .
Figure 7. Map of the study area, with the 113 new seismic zones determined by the combined use of all available data.

Figure 8 .
Figure 8. Map of the three general seismotectonic clusters and the 11 groups in which the 113 new shallow seismic zones were grouped.The code names of the 113 new seismic zones defined in this work are also presented.

Figure 10 .
Figure 10.Variation of the obtained Gutenberg-Richter b values for all 113 shallow seismic zones using the least-squares values (left) and MLM (right) against the magnitude range (dM = M max − M min ) of the complete data set of each zone.Notice the large scatter of values for dM < 1.9 using the least-squares method and the larger b values scatter (even for large dM) of the MLM approach.

Figure 11 .
Figure 11.Map of the spatial distribution of b value based on the 113 shallow seismic zones.Results are presented using data only for zones with M max − M min ≥ 1.9.

Figure 12 .
Figure 12.Spatial distribution of mean return period, T m , for shallow earthquakes of magnitude M ≥ 6.0 within the broader Aegean area.

Figure 13 .
Figure 13.Spatial distribution of the most probable maximum magnitude, M t , of shallow earthquakes for a return period of 50 years within the broader Aegean area.

Table 1 .
Details for the 113 shallow zones of the proposed zonation model.Information for the completeness magnitude since 1981 (additional completeness is also listed, when available) and the estimated b and a 1 (per 10 000 km 2 ) G-R parameters are also presented.