Natural Hazards and Earth System Sciences Probabilistic seismic hazard assessment in Greece – Part 1 : Engineering ground motion parameters

Seismic hazard assessment represents a basic tool for rational planning and designing in seismic prone areas. In the present study, a probabilistic seismic hazard assessment in terms of peak ground acceleration, peak ground velocity, Arias intensity and cumulative absolute velocity computed with a 0.05 g acceleration threshold, has been carried out for Greece. The output of the hazard computation produced probabilistic hazard maps for all the above parameters estimated for a fixed return period of 475 years. From these maps the estimated values are reported for 52 Greek municipalities. Additionally, we have obtained a set of probabilistic maps of engineering significance: a probabilistic macroseismic intensity map, depicting the Modified Mercalli Intensity scale obtained from the estimated peak ground velocity and a probabilistic seismic-landslide map based on a simplified conversion of the estimated Arias intensity and peak ground acceleration into Newmark’s displacement.


Introduction
Probabilistic seismic hazard (PSHA) maps are particularly useful to present the potentially damaging phenomena associated with earthquakes.One important consideration in compiling such maps is which set of ground shaking measures shall be mapped to best serve future engineering needs.Currently, the majority of hazard maps are based upon peak ground acceleration (PGA).This is a simple and convenient way to characterize the ground shaking for many purposes, but it is often found that PGA fails to furnish information about the important characteristics of ground shaking.Characteristics such as duration, frequency energy content and seismic pulse sequences are very important for measur-Correspondence to: G-A.Tselentis (tselenti@upatras.gr)ing the earthquake damage potential and provide the earthquake engineers with crucial design criteria.
Earthquake damage potential increases with amplitude however the relation is complex because of the nonlinear inelastic response of sedimentary deposits and structures to damaging levels of motion.Structures, and in some places, sedimentary deposits responding to ground shaking in a resonant manner cause relatively large deformations and stresses to result if the shaking include several cycles of motion with frequencies close to the resonant frequencies of the structure or deposit.
Ground shaking poses a hazard not only for structures, but may also trigger tsunami, landslides, liquefactions, rockfalls, fire, and in many cases these secondary effects account for a significant proportion of the total earthquake damage.From an engineering perspective, there is a continuous trend to improve the measures of the ground shaking damage potential and these improved measures of ground shaking can provide the input for a more appropriate decision making process in earthquake risk mitigation research.
In recent years, two ground motion parameters have got extensive geotechnical and structural application: Arias intensity (I a ), Arias (1970) and cumulative absolute velocity (CAV), (EPRI, 1988).I a , was proposed as a measure of earthquake intensity based on the instrumental records and was found to be the most efficient intensity measure of the earthquake induced landslide as well as liquefaction potential.I a correlates well with the Newmark's displacement and adequately characterizes the stiff, weak slopes (Bray, 2007).I a is calculated through the integration over the entire length of the acceleration time history described by: where a(t) is the recorded ground acceleration; g is the acceleration due to gravity (both in m/s 2 ); T is the total duration of earthquake motion (in s).
G-A. Tselentis and L. Danciu: Probabilistic seismic hazard assessment in Greece The second ground motion parameter namely CAV, was found to be related with the onset of structural damage and can be used for determining the exceedance of the Operating Basis Earthquakes (OBE) described by (Reed and Kassawara, 1990) whereby the calculation implies the sum of the absolute value of the recorded acceleration.A recent study of Mitchell and Kramer (2006) has shown that a modified version of CAV, computed with a 0.05 g acceleration threshold hereby indicated as CAV 5 , appears to better reflect the longer period (low frequency) components of the motions.Since pore pressure generation is known to be closely related to strain amplitude, which is proportional to particle velocity reflecting longer period components of a ground motion than PGA, CAV 5 might have a closer relation to pore pressure generation than PGA and I a .CAV 5 can be calculated as: In the framework of landslides hazard analysis I a is considered as one of the basic ground motions measures.Abdrakhmatov et al. (2003) have computed a probabilistic seismic hazard analysis in terms of Arias intensity for the territory of Kyrgyzstan, and Peláez et al. (2005) for south-eastern Spain.In addition Del Gaudio et al. (2003) have proposed a way to incorporate the time factor in seismic landslide hazard assessment, based on a predictive model of I a , critical acceleration and Newmark displacement.
It is well known that Greece is characterized with high seismic hazard (e.g.Stavrakakis and Tselentis, 1987), and the assessment of seismic hazard in terms of I a and CAV 5 will provide additionally insight to understand the hazard associated to ground shaking as well as to the possible trigger of landslides or liquefaction.
A first attempt in this respect was performed by Tselentis et al. (2005) who developed probabilistic seismic hazard maps in terms of I a for Greece.The seismic hazard maps were determined using rock site conditions for a period of 50 and 100 years with 90% probability of non-exceedance were based on a ground motion predictive model empirically derived from a limited number of records.
The major aim of the present investigation is to incorporate the engineering ground motion parameters into a consistent seismic hazard analysis.In addition to I a and CAV 5 we also have selected PGA and peak ground velocity PGV, because of their common and widespread usage.Results referring to the acceleration response spectra (SA) and elastic input energy spectra (VE i ) are presented in Part 2 of this investigation (Tselentis et al., 2009).

PSHA methodology
In this study a probabilistic approach initially proposed by Cornell (1968) and improved by Esteva (1970), is used to estimate the seismic hazard in Greece in terms of various ground motion parameters.Detailed review of the existing PSHA methodology is given by the extensive literature on this subject (Kramer, 1996;McGuire, 2004;Reiter, 1990;Thenhaus and Campbell, 2003) and only some basic elements of the PSHA will be recalled for the scope of this investigation.
The methodology follows the four steps: I) sources identification; II) assessment of earthquake recurrence and magnitude distribution; III) selection of ground motion model and IV) the mathematical model to calculate seismic hazard.The state of the practice is to represent the temporal occurrence of earthquakes as well as the occurrence of ground motion at a site in excess of a specified level by a Poisson process.We also assume that: (i) earthquakes are spatially independent; (ii) earthquakes are temporally independent; and (iii) probability that two seismic events will take place at the same location and at the same time approach zero.
Assuming a Poissonian model, the probability P (Y > y * ;t) that at a given site a ground motion parameter, Y , will exceed a specified value, y * , during a specified time period, t, is given by the expression: where v Y (y * ) is the annual frequency of exceeding ground motion level y* and this may alternatively be expresses as: where λ(m i ) is the frequency of earthquakes on seismic sources "i" above a minimum magnitude of engineering significance (m 0 ); P Y ≥ y * |m,r,ε is the probability that, given a magnitude M i earthquake at a distance R i from the site, the ground motion exceeds a value y * ; f M (m) i represents the probability density function associated to the likelihood of magnitude of events (m 0i <M i <m ui ) occurring in source "i"; f R (r|m) i is the probability density function used to describe the randomness epicenter locations within each source "i"; and f E (ε) i is the probability density function that the ground motion deviates epsilon (ε) standard deviations from its median value.The later, represents the number of standard deviations above or below the median ground motions estimated from a predictive equation of the following form: where f M,R,θ fault,site is the functional form of the ground motion predictive model as a function of magnitude, distance, θ fault,site represents a vector of explicit model parameters (fault mechanism, soil category, stress drop, focal depth, slip distribution, etc.), σ ln(Y ) is the standard deviation of the predictive model and epsilon (ε) is expressed as: The functional form of f M (m) i is derived from the unbounded Guttenberg-Richter (1954), relation, which implies that the earthquake magnitudes are exponentially distributed and leads to infinite energy release.In practice, this relationship is truncated at some lower and upper magnitude values which are defined as the truncation parameters related to the minimum (m 0 ) and maximum (m u ) values of magnitude, obtained by different methods in the region under analysis (Cornell and Vanmarcke, 1969).The truncated exponential density function is given by: The function f R (r|m) i for source-to-site distance is computed conditionally on the earthquake magnitude and is obtained by discretization of the cumulative distance probability relationship using a suitable step size.Equation (4) gives the total annual frequency of exceedance, or its reciprocal return period, of each different ground motion level for each ground motion parameter of interest.This relationship between ground motion level and annual frequency of exceedance is called ground motion hazard curve -which is traditionally the standard PSHA output.Another output is the constant-probability, or uniform hazard spectra (UHS) which illustrates the variation of the response spectra amplitudes at a constant return period.

Seismic sources
The broad area of Greece constitutes the overriding plate of the Africa-Eurasia convergent plate system, defining one of the most active plate tectonic regimes in Western Europe and is characterized with high seismic hazard.In brief, the most prominent features of tectonic origin are from south to north, the Mediterranean Ridge, the Hellenic trench, the Hellenic arc, which consists of the outer sedimentary arc and the inner volcanic arc, and finally the back-arc Aegean area, which includes the Aegean Sea, the mainland of Greece, Albania, south FYROM, South Bulgaria and Western Turkey.
The generally accepted assumption in the PSHA methodology is that past seismicity predicts future seismicity and this relies heavily on the quality of the employed seismicity catalogues for such investigations.The seismicity of Greece has been extensively studied, and the major sources of earthquake data are the catalogues and bulletins of ITSAK and the researched catalogues of Papazachos et al. (2000) and Burton et al. (Burton et al 2004).These catalogues span the twenty century for Greece giving an accurate description of seismicity in the region and are shown to be sufficiently homogeneous in magnitude to facilitate the application to PSHA and seismic zoning for Greece.
The catalogue used in this study comprises information on a large number of Greek earthquakes with the time span covering from 550 BC through June 2008, as was published by the Geophysical Laboratory of the University of Thessaloniki.The catalogue consists of a large number of shallow events, homogeneous in term of moment magnitude, herein denoted M. For the region of investigation, the catalogue is complete for magnitudes M ≥8.0 since 550 BC, for M ≥7.The reported uncertainty for the epicenters associated with the historical (prior to 1911) and instrumental events in the catalogue is 30 km and 20 km, respectively.The uncertainty in magnitude estimation is 0.25 magnitude units for the instrumental catalogue, while for the historical catalogue the uncertainty of magnitude may rich 0.5 magnitude units (Papazachos et al., 2000).
Figure 1, depicts the high seismicity of Greece, and this information cannot be related always to specific faults precisely enough for use in seismic hazard analysis.This is due to the occurrence of a large number of events offshore.Therefore, the location of possible earthquakes is represented by seismic zones.Papaioannou and Papazachos (2000), on the basis of the available seismicity, geological and geomorphologic information, have proposed a delimitation of Greece and adjacent regions in 67 seismogenic sources of shallow earthquakes.We have selected these seismic zones as an optimum representation for adequate spatial distribution of the seismicity in the region.Moreover, we have investigated the type of faulting in each seismogenic source zone and the dominant fault mechanism was accordingly assigned.A map depicting these seismogenic source zones were adopted herein and the geographical boundaries of the zones and delimitations along the dominant fault mechanism for each one are present in Fig. 1.Epicenter locations and magnitudes of all earthquakes from the entire catalogue with magnitude of 5.0 M or greater are also shown.
Surface sources instead of point or line sources are used in the present PSHA analysis.Seismic sources are modeled with the parameters reflecting the characteristics of each source.These seismicity parameters include: the a-and b-value parameters of the Gutenberg-Richter frequency relationship in the area of each source, the lower (m 0 ) and upper bound magnitude (m u ) and the annual rate λ(m) of exceedance of earthquakes with M>5.An attempt was done to compute these seismicity parameters for the selected catalogue.The adopted parameters to characterize the seismicity of each seismic zones were obtained from the study of Papaioannou and Papazachos (2000).Because their catalogue is until December 1999 and we wanted to investigate the influence that later events might have on the final values of the seismicity parameters.
For this purpose, we have declustered only the instrumental catalogue using the algorithm proposed by Reasenberg (1985).Working on the declustered instrumental catalogue we have computed the earthquake recurrence statistics for the most active seismic zones using the unequal observation period for different magnitude ranges and the maximum likelihood method (Weichert, 1980).Surprisingly, we found that the differences between the selected and the new values of the seismicity parameters are negligible.Therefore we have selected for each source the seismicity parameters proposed by Papaioannou and Papazachos (2000).No uncertainty was taken into account in the magnitude-frequency parameters therefore these relationships only provide a median estimate of seismicity rates.

Minimum magnitude
The lower bound earthquake magnitude represents the threshold to remove the non-damaging earthquakes from the hazard analysis.In this study this threshold is set at M=5 for all source zones, because earthquakes bellow this magnitude are not potentially damaging for well-engineered structures Recently, the study by EPRI (2006) indicates that the PSHA results are sensitive to the selection of the lower bound magnitude, and they have proposed a criterion based on CAV to cut-off these small non-damaging events.This CAV parameter was found to be the most suitable parameter for use in predicting the threshold of ground motion damage potential (Reed and Kassawara, 1990).
The damage potential threshold was conservatively formulated to ensure that no damage will occur to buildings of suitable design and construction and a CAV value of 0.16 g-s was found.In the framework of PSHA, the probability of exceeding a CAV value of 0.16 g-s is used to remove earthquakes that are not potentially damaging from the hazard analysis.The CAV filtering procedure requests region specific relationships to estimate CAV from few independent parameters including PGA, duration, amplification factors, magnitude and distance.At this time, it is pointed out that there are no previously published results on such relationships for the region of Greece.

Maximum magnitude
This is an important source parameter which can dominate the ground motion assessments particularly at long return period hazard (low-probability hazard).Despite of its importance, the estimation of the maximum magnitude expected in a source zone is not a straightforward task.Several approaches have been proposed including the maximum historical earthquake procedure and the likelihood method proposed by Kijko and Sellevoll (1989).For the present investigation the maximum magnitude for each source region was retained from the study of Papaioannou and Papazachos (2000).
Following the rule-of thumb in the hazard analysis we have increased the values by half-unit in the magnitude scale to yield the maximum expected magnitude.The effect of the maximum magnitude on the final PSHA results will be investigated using a sensitivity analysis of the hazard results.
Another source parameter is the focal depth of the earthquakes.This parameter is poorly known for historical events, due to lack of data.However this investigation focuses on the shallow earthquakes and we have neglected deep events, because few data is available for them.A mean depth of approximately 10 km was assigned to all sources.

Ground motion model
The selection of ground motion predictive equations suitable for a region of interest is of great importance, and a seismic hazard analysis must consider all the potentially applicable ground motion predictive equations for that region.Cotton et al. (2006) have defined a set of criteria to justify the selection of appropriate ground motion models for a particular target area.The criteria of selecting ground motion models were selected on the basis of rejection reasons: (i) the model is derived from a clearly irrelevant tectonic regime; (ii) the model is not published in an international peer-reviewed journal; (iii) the documentation of the model and the underlying dataset is insufficient; (iv) the frequency range and the functional form of the model is not appropriate; and (v) the regression method and regression coefficients are inappropriate.We have adopted these criteria as guidance in our selection process and the selected ground motion models are described in the next paragraphs.
According with the aim of the present investigation, we have primarily selected the regression models proposed by Danciu and Tselentis (2007), as appropriate to estimate the set of ground motion parameters.They have empirically derived from a set of 335 records from 151 shallow earthquakes a set of predictive equations for various ground motion parameters, including PGA, PGV, I a , CAV, CAV 5 , SI, SA, and VE i .For a reliable PSHA, the rule of thumb is to consider at least two alternative predictive equations.It is worth noticing that for the region of interest, there is a lack of regional predictive equations for ground motion parameters, such as I a or CAV, and therefore alternative equations might be considered.For PGA and PGV we have selected the regression models proposed by Margaris et al. (2002) and Skarlatoudis et al. (2003).These candidate predictive equations were empirically derived from the data recorded in Greece in the last decades.
The selected ground motion predictive equations are compatible in terms of magnitude scale and source-to sitedistance definition and no additional conversion is required.The magnitude scale involved in the ground motion predictive equations is moment magnitude M, and the source to site distance is reported as the epicentral distance.The equations also take into account the local site conditions.The main difference between the equations appears from the type of-faulting parameter.
The model proposed by Margaris et al. (2002) was derived mainly using events with normal fault mechanism, whereas the other two models explicitly incorporate the typeof-faulting parameter.Another difference between the two models may arise from the way the horizontal components are treated in the regression process, i.e. arithmetic mean of the two components (Danciu and Tselentis, 2007), including both components as individual data points (Margaris et al., 2002;Skarlatoudis et al., 2004).This implies that the later two predictive models would generally bias as to the increase the results in the estimated peak ground acceleration or velocity compared with the former due to the increased variance.
For CAV 5 , the only alternative predictive model available is the one adopted herein and proposed by Kramer and Mitchell (2006).It has to be mentioned that these alternative predictive models for I a and CAV 5 were in agreement with the selection criteria mentioned above.The models are actually similar and, well documented and were obtained from a large dataset.The models have the same functional form, including the fault mechanism factor, and are homogeneous in terms of moment magnitude.The models have used the arithmetic mean of the two horizontal components and the same regression technique.In the following sections the next acronyms will be used to refer the selected predictive equations: Margaris et al. (2002) denoted as MA02; Skarlatoudis et al. (2004) denoted as SK04, Danciu and Tselentis (2007) denoted as DT07, Travasarou et al. (2003) denoted as TR03, Mitchell and Kramer (2006) denoted as MK06.
The main differences between the TR02 and MK05 predictive equations and DT07 arises from different source to site distance definition.The former uses source-to-site distance defined as the closest distance to the seismic fault, which is the rupture distance, whereas the later uses epicentral distance.This difference may lead to major differences in the resulted ground motion in the near field and large events.Scherbaum et al. (2004) have presented a statistical model to convert different source-to-site distances for extended sources.They have shown that the conversion of distances increases the variability of the modified ground motion predictive model and so become magnitude and distance dependent.To overcome this, they have proposed empirical laws for metrics conversion based on simulated scenarios.
Although this method is promising, we have not considered it in this investigation mainly because it implies definitions of scenarios based on converting magnitude and rupture characteristics (length, width and area) and surface displacement which is not always evident for the region of Greece dominated by offshore events.These objections are valid mainly for near field and large magnitude events, for small events, the epicentral distance is equal to or larger than the rupture distance.This will imply that the predictive equations which uses the epicentral distance will yield a higher estimated ground motion than that if the rupture distance pawww.nat-hazards-earth-syst-sci.net/10/1/2010/Nat.Hazards Earth Syst.Sci., 10, 1-15, 2010  rameter is used (Reiter 1990).In this case the equations of TR03 and MK06 would evidently underestimate the I a and CAV 5 values at near field distances.
A simplified attempt to convert the rupture distance to the epicentral distance is also made in this study.A dataset of 112 horizontal components recorded in Greece with reported epicentral and rupture distances was gathered.We then fit a linear regression on the selected dataset and we have obtained the following relations between the two distance definitions: where R E is the epicentral distance and R r is the rupture distance in km.
The predictive equations involved in the present hazard computation are described in Table 1.Since the input of hazard computations is required as a matrix of median estimated ground motion predictive equations, we have represented the predictive equations as 3-D surface plots, rather than as a log-log representation, as depicted in Fig. 2, which shows the overall behavior of the selected predictive equations and also reveals the feasibility of identifying regions of equal intensity for different combinations of magnitude and distances.form.The uncertainties associate to the ground motion can be characterized as aleatory variability and epistemic uncertainty and it was found to have great influence of the final PSHA results, particularly at low annual frequencies of exceedance (Bender, 1984;Bommer and Abrahamson, 2006).Aleatory variability is introduced by true randomness in nature while epistemic uncertainty is due to lack of knowledge.These types of uncertainties are treated differently in advanced PSHA: aleatory variability is explicitly integrated in the hazard calculation, while epistemic uncertainty is treated by multiple hypotheses, models, or parametric values.The later, was already quantified by considering various ground motion predictive models, which were incorporated in a pseudo-logic tree approach.For the purpose of this study it is called pseudo-logic tree approach because it has only one decisional node, at the level of predictive equations.The logic tree approach requests that the sum of the probabilities associated to each branch to approach unity.Therefore, for PGA and PGV models, a 0.35 probability to SK04 and DT07 and a 0.3 probability to MA02 was attributed; for I a models the DT07 was assigned to a 0.60 probability while TR03 to a 0.40; for CAV models the probability associated DT07 was 0.40 while for MK05 was 0.40.
The aleatory (random) variability associated with the ground motion predictive equations can be characterized by sigma.Sigma (σ ) is the standard deviation of the residuals, which are generally assumed to follow a log-normal distribution.It has a great influence on the final PSHA results, for a given level of seismicity, effectively controlling the shape of the seismic hazard curve.
It has become an obligatory practice to incorporate σ into the PSHA calculation, and to impose some limits on the maximum value of ε.More specifically, if σ is incorporated into the hazard integral it was found that the lognormal distributed ground motion does not saturate but grows infinitely particularly at long return periods.As a result, the lognormal distribution for ground motions is truncated at a specified number of standard deviations, ε.
The ε values are arbitrary adopted and most often values of 2 to 4 have been proposed, however there is a lack of consensus concerning the physical basis of these values.The study conducted by EPRI (2006) concluded that there is no basis of truncating the ground motion distribution at an epsilon value of less than 3 and there are observations of ε values greater than 3. Similar findings were recently highlighted by Strasser and Bommer (2008) and Abrahamson and Bommer (2008).They have investigated the physical basis for choosing the maximum number of ε, by examining the large residuals from the dataset used to derive one of the ground models of the Next Generation Attenuation (NGA) project.Their conclusions were that it is not possible to obtain direct physical constrains on the epsilon, even for well-document data.Both studies have emphasized that the probability plot of the residuals may provide insight and give guidance about the level of truncation.

Seismic hazard computation
In order to obtain the probability of exceedance of the selected ground motion parameters, it is necessary to convolve over all possible magnitudes, distances and values of ground motion parameters in a probabilistic manner.The complete treatment of these aspects requires complex specialized computer packages and the widely accepted and used are EQRISK (McGuire, 1976), SEISRISK III (Perkins, 1998), FRISK (McGuire, 1978) and Crisis2003 (Ordaz et al., 2003).In the present investigation, we have selected the software Crisis2003 version 3.1.This software allows assigning different predictive equations to different seismic source zones.
We have decided to estimate the ground motion parameters in 50 year time period, corresponding to the design lifetime for buildings and we have chosen 10% of exceedance that would lead to a return period of 475 years.We have not incorporated the local site condition on the hazard computation therefore the results are valid for an ideal "bedrock" (V s >800 m/s) local site conditions.Since the predictive equations take into account the fault mechanism, we have identified the predominant type-of-faulting for each seismogenic source, as shown in Fig. 1.The fault mechanism is the criterion used to assign the predictive equations to the corresponding seismogenic sources.Since the seismic hazard is site specific, rather than region specific, the geographical territory of Greece spanning the area 19 • W-30 • E and 34 • S to 42 • N was divided into a mesh of grid points with an interval of 0.1 • (about 10 km) both in latitude and longitude.The integration over all the grid points was carried out by running the Crisis2003 computer package and the results are presented in the next section.

Seismic hazard results
Prior to the presentation of the final PSHA results, the implication of some factors, including minimum and maximum magnitude, uncertainty in the ground motion estimation, truncation levels have to be quantified.Among all these factors, it was found that the strongest influence on the hazard results is exhibited by the uncertainty in the ground motion estimates.An average increase of almost 35% in all ground motion estimates was observed when the standard deviation of the predictive model was taken into account.However, it has become a mandatory practice to incorporate the ground motion variability into the PSHA.Increasing the maximum magnitude by half-unit we find the hazard to increase by 18%, whereas by lowering the minimum magnitude by half-unit the hazard is found to increase by 10%.A decrease of 5% of the seismic hazard was also observed when the rupture distance was converted to an epicentral distance.
The corresponding hazard maps were obtained by interpolation between the grid points with the contouring algorithm implemented in the GMT mapping package (Wessel and Smith, 1998).
From the derived hazard maps we can identify regions of high seismic hazard: (i) in the western Hellenic arc with the maximum hazard reported for various ground motion parameters in the Cephalonia Island; (ii) in the northern Aegean, at the intersection of two different tectonic regimes in the extension of the Aegean Arc and the North Anatolian Fault; (iii) in central Greece which includes the region of Corinth Gulf; and (iv) in the South-Western part of Crete Island.It has to be mentioned here that the ground motion parameters estimated in the South-Western part of Crete Island may be underestimated since intermediate depth earthquakes which dominate this region, were not taken into account since that only shallow events were considered.
Additionally to the hazard maps we have determined and present in Table 2 the estimated ground motion parameters for 52 Greek municipalities.Among these municipalities, the most prone to high seismic hazard is the municipality of Argostolion, with the highest values for the estimated ground motion parameters.The seismic hazard results presented in  Figure 4 illustrates the spatial distribution of the mean PGA values computed at each grid point based on the weighted predictive models MA02, SK04 and DT07 and produced for a return period of 475 years.The region of the highest estimated values, greater than 0.45 g is observed in the region of the Ionian Islands, which is the most seismic active region in Greece.
At this point it might be of interest to validate our probabilistic hazard analysis by a comparison with other studies that have used different approaches.In this respect, we have selected the studies of Burton et al. (2003) and Mäntyniemi et al. (2004).Burton et al. (2003) assessed seismic hazard in Greece using the Gumbel's asymptotic distributions of extreme values, and two different completeness periods for the Greek catalogue: 1900Greek catalogue: -1999Greek catalogue: and 1964Greek catalogue: -1999. .The median PGA values were estimated for a return period of 475 years (they reported as 90% non-occurrence in 50 years) for six Greek cities.The corresponding median PGA estimated on stiff soil with the aim of the predictive equations proposed by Theodulidis and Papazachos (1992) were found to be: Athens -0.24 g, Thessaloniki -0.22 g, Patras -0.23 g, Corinth -0.36 g, Rodhos -0.28, and Heraklion -0.29 g.Mäntyniemi et al. (2004) have proposed probabilistic seismic hazard curves in terms of PGA for five Greek cities including Athens, Thessaloniki, Patras, Volos and Heraklion by using the parametric-historic procedure.The corresponding median PGA values are based on the Margaris et al. (2002) predictive model and the result for each one of these cities are: 0.24 g for Athens, 0.53 g for Heraklion, 0.30 g for Patras, 0.35 for Thessaloniki, and 0.30 g for Volos.Comparing these values with those presented in Table 2 and given the differences in the various components of the analyses, the results are notably similar and corroborative.The updated version of the Greek design code (EAK, 2003) divides the country in three seismic zones, namely Zone I (PGA=0.16g) Zone II (PGA=0.24g), and Zone III (PGA=0.36g) as can be seen in Fig. 5. Comparing this map with the probabilistic estimated mean PGA map it can be observed that the latest highlights an increased seismic hazard.
The PGV maps presented in Fig. 6 show a similar pattern with the mean PGA maps and minor differences exist towards the southern part of Crete Island.High hazard however is also observed in the Ionian Islands, as is observed for the case of PGA.The probabilistic PGV maps are more sensitive to high magnitudes than the PGA maps, relevant for a frequency band of about 1 Hz and thus provide information on the relative levels of expected shaking for larger, more flexible structures with lower resonant frequencies.However, PGV was found to be a very good descriptor of earthquake intensity, and a versatile parameter to compute instrumental intensity for rapid damage assessment tools such as ShakeMaps (Wald et al., 1999).
With the ShakeMap approach, a rapid estimation of the instrumental intensity can be accomplished with the help of the estimated PGV in a fictitious network.With the help of the recently derived relationships between PGV and Modified Mercalli Intensity scale-MMI proposed by Tselentis and Danciu (2008) we can convert the probabilistic PGV values to a probabilistic MMI map.Without neglecting the large dispersion introduced in such an approach, the probabilistic intensity map was obtained by applying the following equation at each grid point: The probabilistic MMI map obtained from the estimated median PGV values is presented in Fig. 7.This type of MMI maps are useful in comparing the hazard results with historical information, because the later spans a long time period.It is of our interest to compare the present estimated MMI values for several cities with the probabilistic MMI values obtained by Papaioannou and Papazachos (2000).We observe that the values reported in the present study are slightly larger than the probabilistic MMI.For Athens and Thessaloniki, we obtained respectively 7.41 and 7.38 while Papaioanou and Papazachos (2000) reported 7.12 and 7.17 for the same cities.
The highest MMI value, 9.08 was obtained for the city of Argostolion, which is slightly larger than the value of 8.43 reported on the same study.
Another engineering significance of the probabilistic PGA and PGV maps arises from the investigation of the PGV/PGA ratio, defined as a frequency index (Tso and Zhu, 1992).
The PGV/PGA ratio is direct proportional with the width of the acceleration region of an elastic response spectrum in the tripartite logarithmic format and it is directly connected to the corner period (T c ) of the design spectrum.As a result a higher PGV/PGA ratio provides for a wider accelerationsensitive region.It is also observed that the near-field records with directivity effects tent to have high PGV/PGA ratio which can influence seriously their response characteristics.Based on the two probabilistic mean PGA and PGV maps, the range of the PGV/PGA ratio was found to be ranging from 0.05 s to 0.15 s.That is a relatively moderate value covered by the width of the acceleration region presented in EAK (2003) for rock local site condition.
The maps depicting the mean of the estimated I a are presented in Fig. 8.There is not much spatial variation of the distribution of the seismic hazard estimated in terms of I a as compared with the PGA and PGV parameters.The mean I a exceedance values for the entire region oscillate from 0.1 m/s to 7 m/s.Keefer and Wilson (1989) proposed three groups of slope instability based on the threshold values of I a : 0.11 m/s -falls, disrupted slides and avalanches; 0.32 m/s -slumps, block slides and earth flows; and 0.54 m/s -lateral spreads and flows.Considering these values, it appears that in the regions where high topographic relief is combined with the high I a estimated values there is significant slope instability potential.However, the mean I a map shown in Fig. 8 represents a rather simplified way to characterize the seismiclandslide hazard, because there is not an immediate connection between the level of shaking and its effects on slope instability.Also, factors such as the local ground conditions and topographic relief have to be taken into account.
A measure of the permanent displacement caused by shaking along a slide surface was proposed by Newmark (1965).The Newmark's displacement caused by earthquakes can  down-slope movement, and PGA is peak ground acceleration in g.We have arbitrary assumed that the ratio between the a c and PGA is fixed at 0.5, and for each grid point we have estimated the D n using Eq. ( 10) and the previously estimated I a and PGA values for a return period of 475.
The probabilistic Newmark's displacement map obtained as a function of ground shake intensity (I a ) and dynamic slope stability (ratio between critical acceleration and PGA) is presented in Fig. 9. Wilson and Keefer (1985) suggested threshold values of D n equal 10 cm, for triggering coherent slides (slumps, block slides, slow earth flow) and 2 cm for disrupted slides (rock falls, rock slides, rock avalanches).Judging from the D n map it appears that the whole region exhibits a prone landslide risk, particularly for disrupted slides, since the threshold value of 2 cm is exceeded everywhere.It can be observed that in the region of Corinth Gulf there is a high risk of coseismic landslide hazard for both coherent and disrupted slides (Tselentis and Gkika, 2005).However, the landslide hazard estimation is heavily slope-driven, because the dynamic slope stability has large spatial variability in nature even within geologic units.
Due to the arbitrary assumption of the acceleration ratio, the presented seismic-landslide hazard map is only an approximation.A genuine seismic-landslide hazard evaluation relies on (i) a detailed inventory of slope processes, (ii) the study of those processes in relation to their environmental setting, (iii) the analysis of conditioning and triggering factors, and (iv) a representation of the spatial distribution of these factors (National Research Council, 1996).This is beyond the aim of the present study, whereas we have portrayed an example of rapid conversion of probabilistic PGA and I a maps into an engineering product: a landslide susceptibility map.The spatial distribution of the ground shaking described by CAV 5 is illustrated in Fig. 10.As was expected, due to the fact that hazard is governed by the seismic zonation, there are no large difference between the pattern of CAV 5 values and the other ground motion parameters.The spatial distribution of the high seismic hazard is corroborative with the previous hazard maps, with the maximum hazard occurring in the Cephalonia Island.Very high values, as large as 100 m/s were obtained in the region of Ionian Islands, as well as in the region of Corinth Gulf.
Because CAV 5 was found to be a suitable liquefaction descriptor, the probabilistic CAV 5 map may be seen as a liquefaction-assessment map that depicts how often a level of ground shaking, sufficient to cause liquefaction, is likely to occur.The CAV 5 map can be combined with a liquefactionsusceptibility map that shows the distribution of sedimentary units, and therefore to obtain a probabilistic map of the liquefaction-potential for the given region.Papathanassiou et al. (2003) focused their study on the coseismic effects such as ground failures like rock falls, soil liquefaction, ground cracks and slope failures over the island of Lefkada in the Ionian Sea.They have observed that these phenomena appeared all over the island, rock falls are more concentrated on the northwestern edge and liquefaction on the quaternary deposits at the northern and eastern island coasts, which are considered to be of high potentially liquefiable type.Also, they observed evidences of liquefaction occurrence as sand boils and ground fissures with ejection of mud and water mixture at the waterfront areas in the town of Lefkada and in the villages of Nydri and Vassiliki.For Lefkada city the present study estimated a mean values for I a and CAV 5 equal to 4 m/s and 60 m/s, respectively.

Conclusions
A PSHA was conducted for the region of Greece.The PSHA relies on the seismic regionalization methodology, assuming that the future location of major seismic activity will be limited to the boundary of the seismotectonic zones proposed by Papaioannou and Papazachos (2000); and the ground shaking was described by various ground motion parameters: PGA, PGV, Arias, and CAV 5 .The output of the regional hazard analysis was consisted of a set of probabilistic hazard maps for PGA, PGV, I a and CAV 5 .The probabilistic approach was based on the well known Cornell methodology (1968).
The epistemic and the aleatory uncertainty was incorporated in the analysis including the standard deviation into the direct computation and considering more than one predictive model for the ground motion parameters.The PGA and PGV models included the regression model proposed by Margaris (2002), Skarlatoudis et al. (2004), and Danciu and Tselentis (2007).Along with the Danciu and Tselentis ( 2007) predictive equations we have selected the equations proposed by Travasarou et al. (2003) for I a and Mitchell and Kramer (2006) for CAV 5 , respectively.The predictive equations were assigned to each zone according to the predominant fault mechanism, and the ground motion parameters were estimated on an ideal rock category.
The output of the hazard computation consisted of probabilistic hazard maps estimated for a fixed return period of 475 years.From these maps the estimated parameters values were extracted and reported for 52 Greek municipalities.Additionally, we have obtained a set of probabilistic maps with engineering significance: a probabilistic macroseismic intensity map, depicting the Modified Mercalli Intensity scale (MMI) obtained from the estimated PGV and a probabilistic seismic-landslide map based on a simplified conversion of estimated I a and PGA into Newmark's displacement.
The present study relies on the seismogenic source zones and compares well with the free zoning approach used by Burton et al. (2003), and with the parametric historic procedure used by Mäntyniemi et al. (2004), despite the difference in the hazard assessment methodologies.The major limitation of the present methodology is that the PSHA relies on the seismogenic source zones and the associated statistics.A rigorous estimation of the maximum potential magnitude and the consideration of the local site conditions could lead to more realistic results.

Figure 3
Figure 3 depicts the probability plot of the residuals of the DT07 predictive equations.It can be observed that the residuals follow a lognormal distribution very closely, right up to values of ±3σ .However, for the present investigation we have imposed the truncation level at ε=3 for all the selected ground motion predictive equations.

Table 2 .
Probabilistic engineering ground motion parameters values for major Greek municipalities.
Table 2 allow rapid comparison with previous hazard assessment proposed by various studies.All the results are calculated having 10% chance of being exceeded over 50 year exposure period corresponding to 475 year return period.www.nat-hazards-earth-syst-sci.net/10/1/2010/Nat.Hazards Earth Syst.Sci., 10, 1-15, 2010