An alternative approach to probabilistic seismic hazard analysis in the Aegean region using Monte Carlo simulation

Tectonophysics 492 (2010) 253–278

Contents lists available at ScienceDirect

Tectonophysics j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / t e c t o

An alternative approach to probabilistic seismic hazard analysis in the Aegean region using Monte Carlo simulation Graeme Weatherill a,b,⁎, Paul W. Burton a a b

Seismic Risk Group, School of Environmental Sciences, University of East Anglia, Norwich, NR4 7TJ, United Kingdom European Centre for Training & Research in Earthquake Engineering (EUCENTRE), Via Ferrata 1, 27100, Pavia, Italy

a r t i c l e

i n f o

Article history: Received 23 October 2009 Received in revised form 22 June 2010 Accepted 30 June 2010 Available online 8 July 2010 Keywords: Seismic hazard assessment Monte Carlo simulation Aegean Epistemic uncertainty

a b s t r a c t The Aegean is the most seismically active and tectonically complex region in Europe. Damaging earthquakes have occurred here throughout recorded history, often resulting in considerable loss of life. The Monte Carlo method of probabilistic seismic hazard analysis (PSHA) is used to determine the level of ground motion likely to be exceeded in a given time period. Multiple random simulations of seismicity are generated to calculate, directly, the ground motion for a given site. Within the seismic hazard analysis we explore the impact of different seismic source models, incorporating both uniform zones and distributed seismicity. A new, simplified, seismic source model, derived from seismotectonic interpretation, is presented for the Aegean region. This is combined into the epistemic uncertainty analysis alongside existing source models for the region, and models derived by a K-means cluster analysis approach. Seismic source models derived using the K-means approach offer a degree of objectivity and reproducibility into the otherwise subjective approach of delineating seismic sources using expert judgment. Similar review and analysis is undertaken for the selection of peak ground acceleration (PGA) attenuation models, incorporating into the epistemic analysis Greek-specific models, European models and a Next Generation Attenuation model. Hazard maps for PGA on a “rock” site with a 10% probability of being exceeded in 50 years are produced and different source and attenuation models are compared. These indicate that Greek-specific attenuation models, with their smaller aleatory variability terms, produce lower PGA hazard, whilst recent European models and Next Generation Attenuation (NGA) model produce similar results. The Monte Carlo method is extended further to assimilate epistemic uncertainty into the hazard calculation, thus integrating across several appropriate source and PGA attenuation models. Site condition and fault-type are also integrated into the hazard mapping calculations. These hazard maps are in general agreement with previous maps for the Aegean, recognising the highest hazard in the Ionian Islands, Gulf of Corinth and Hellenic Arc. Peak Ground Accelerations for some sites in these regions reach as high as 500–600 cm s−2 using European/NGA attenuation models, and 400–500 cm s−2 using Greek attenuation models. © 2010 Elsevier B.V. All rights reserved.

1. Introduction Seismic hazard assessment has been an important objective of earthquake research in the Aegean region for many decades. The seismotectonic complexity of the region combined with the economic and urban growth of Greece, Turkey and many of the Balkan states, affords geologists and seismologists little complacency when considering the nature of earthquakes in the Aegean and their potential socio-economic impact. Seismic hazard maps, in particular, have been crucial in the seismic zonation of Greece and Turkey, which underpins the implementation of seismic resistant building design codes. Analyses of seismic hazard in the Aegean region, here defined as the area enclosed by the 18° E and 30° E meridians and 34° N and 43° ⁎ Corresponding author. Seismic Risk Group, School of Environmental Sciences, University of East Anglia, Norwich, NR4 7TJ, United Kingdom. E-mail address: [email protected] (G. Weatherill). 0040-1951/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.tecto.2010.06.022

N parallels, have developed substantially over the last 30 years largely as a result of the increasing quality and quantity of data pertaining to the seismicity of the Aegean region (Papazachos and Papazachou, 1997; Burton et al., 2004b). The seismic hazard assessment methodology has evolved less markedly than the seismicity data, largely adopting the zonation PSHA methodology (Hatzidimitriou et al., 1985; Papazachos, 1990; Grünthal et al., 1999; Papaioannou and Papazachos, 2000; Tselentis and Danciu, 2010a,b; Tselentis et al., 2010) of Cornell (1968) and McGuire (1976), or an extreme value approach (Makropoulos and Burton, 1985a,b; Burton et al., 2003, 2004a) to analyse seismic hazard. The juxtaposition of alternative methodologies is a useful check on seismic hazard results and limitations. This present paper will evolve the PSHA methodology to a further level for the Aegean region. The primary focus of this work is in the application of an alternative method of probabilistic seismic hazard analysis: the Monte Carlo approach. Early applications of the Monte Carlo method (Shapira, 1983;

254

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

Johnson and Koyanagi, 1988) demonstrated how random simulations of seismicity could be applied for the purposes of seismic hazard analysis. They were, however, limited in scope due to the overall limitations of computational efficiency. More comprehensive applications have since been developed, and have generally followed one of two categories: simulation of synthetic earthquake catalogues via random sampling (with replacement) of observed seismicity (Ebel and Kafka, 1999); or simulation of synthetic earthquake catalogues via random sampling within uniform seismic source zones and from the cumulative distribution function of an appropriate magnitude-frequency recurrence model (Musson, 1999a; 2000; Musson and Sargeant, 2007; Wiemer et al., 2009). In either case seismic hazard is estimated by calculating the expected ground motion at a site from synthetic catalogues of seismicity, with typical durations on the order of 106 or 107 years, using empirical attenuation models. The ground motion with a P % probability of being exceeded in T years is then derived from the model of Poisson distributed ground motions (Y) exceeding a threshold level (Y0). Ebel and Kafka (1999) make the calculation directly from the synthetic catalogues of durations (T0) using the relation:

is the 21 July A.D. 365 Gortyna (South Aegean) event (Mw ≈ 8.3). It contains no events smaller than Mw 6.0. Focal depth information is not available for most of the earthquakes in the catalogue. Where earthquake depths are recognised as being intermediate or deep they have been assigned arbitrary depths, usually 75 km (intermediate) or 150 km (deep). Earthquake magnitudes in the catalogue are homogenised into both the moment-magnitude (Mw) and surface-wave magnitude (MS) scales, of which the former is used in the seismic hazard analysis. Where it is necessary to convert to MS within the synthetic catalogues, an adaptation of the equations of Papazachos and Papazachou (1997) is used:


 −T ∑H ðY≥Y0 Þ P ðY≥Y0 Þ = 1− exp T0

Differences in the properties of the component catalogues and the historical quality of instrumental recording networks and macroseismic observations mean that the completeness of the earthquake catalogues is temporally variable. The homogenisation of several different catalogues results in abrupt changes in completeness magnitude (MC) at certain time periods. Analysis of completeness has been approached using different methods. Initially the Stepp (1971) technique was applied, using magnitude bins of 0.5 Mw width. For each interval the lower bound completeness magnitudes are given in Table 1. Lower bound magnitude completeness for the early instrumental period (1900–1964) is estimated as being Mw 5.0. From 1964, the catalogue corresponds to that of the International Seismological Centre (ISC), so the completeness drops to Mw 4.5. From 1998 records are taken from the National Observatory of Athens (NOA); hence the further drop in completeness to Mw 4.0. These values of MC are estimates given the magnitude bin width, except in the case of the NOA era where actual completeness is lower but the catalogue is truncated above this magnitude. To refine estimates of MC further, other approaches were applied (Wiemer and Wyss, 2000; Woessner and Wiemer, 2005). These results consistently indicate that from 1964 to 2005, MC = Mw 4.8 and for 1900 to 2005 MC ≈ MW 5.2 ± 0.1, and it is these values that are adopted in the following work. For time-independent seismic hazard analysis, non-Poissonian events are removed from the catalogue. The declustering algorithm of Reasenberg (1985) is used for this purpose. A confidence probability of 0.95 is used here, assuming errors on epicentres and depths are 1.5 km and 2 km respectively. This algorithm is implemented using the ZMAP software for statistical analysis of seismicity (Wiemer, 2001). A total of 1285 events were identified as non-Poissonian, and these events were subsequently removed from the earthquake catalogue. This resulting declustered catalogue can be seen in Fig. 1.

ð1Þ

Where H is the Heaviside step function, which assumes a value of 1 if Y ≥ Y0, 0 otherwise. This adapts the Poisson model for application to a discrete distribution of ground motion values, defining the rate of the Poisson model in terms of the number of times ground motion Y0 is exceeded within the synthetic catalogue. Musson (1999a) instead considers the maximum Y in each synthetic year and ranks them in descending order, extracting annual probabilities from the rank of each ground motion level in the synthetic catalogue. In this approach we calculate hazard directly from the Poisson model, in the manner of Ebel and Kafka (1999). The Monte Carlo method has many benefits, some of which have provided the motivation for this work. It is a more transparent process than that of Cornell (1968) — McGuire (1976), allowing for ease of communication to non-technical audiences (e.g. planners and politicians). Its main advantage is the manner in which it can handle uncertainty in the key parameters of the hazard models (Musson, 1999a,b, 2000). It is also flexible in the manner in which different models of seismicity or assumptions can be integrated (Musson, 1999a). This incorporation of uncertainty and model flexibility is something that we intend to explore in this research. 2. The earthquake catalogue 2.1. Composition of the catalogue The instrumental earthquake catalogue used in this analysis is that of Burton et al. (2004a), which covers the area enclosed by the 18°E and 30°E meridians and the 34°N and 43°N parallels. This incorporates all of Greece, Albania and the Former Yugoslav Republic of Macedonia (FYROM), as well as southern Bulgaria and western Turkey. The instrumental catalogue covers the time period 1900 A.D. to 1999 A.D. We have extended the catalogue to 2005 A.D. using data from the National Observatory of Athens Online instrumental catalogue (NOA, 2007 — www.gein.noa.gr/services/info-en.html). Magnitudes for the 2000 to 2005 period are homogenised with the rest of the catalogue using the procedure described in Burton et al. (2004a) and Weatherill (2009). The Burton et al. (2004a) catalogue is truncated at Mw 4.0. This is to remove microearthquakes that are inconsistently reported, whilst retaining events of possible engineering significance. In addition to the instrumental catalogue, a catalogue of historical Aegean earthquakes is included in the analysis (Papazachos and Papazachou, 1997). This catalogue contains historical earthquakes for the period 550 B.C. to 1900 A.D. The largest earthquake in this catalogue

MW = MS

ð2aÞ

for 6:0≤Ms

MW = 0:56MS + 2:66

for 4:2≤Ms ≤6:0

ð2bÞ

2.2. Completeness and declustering

2.3. Maximum magnitude (MMAX) The model of magnitude recurrence preferred here is the doubletruncated Gutenberg and Richter (1944) relation, expressed as a

Table 1 Estimates of completeness magnitude for the Aegean earthquake catalogue using the Stepp (1971) method. Time period (all dates A.D.)

MC

1999–2005 1965–2005 1900–2005 1800–2005 1700–2005

4.0 4.5 5.0 6.0 6.5

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

255

Fig. 1. Aegean earthquake catalogue declustered using the Reasenberg (1985) algorithm.

cumulative distribution function (Kramer, 1996; Utsu, 1999; Tselentis and Danciu, 2010a,b; Tselentis et al., 2010): FM = P ½MbmjMmin ≤m≤MMAX  =

1− exp½−βðm−Mmin Þ 1− exp½−βðMmax −Mmin Þ

ð3Þ

For each seismic source zone, and for the Aegean as a whole, it is necessary to define an upper bound magnitude, and a lower bound as we shall consider later. This may be constrained from physical fault dimensions, or else by consideration of the statistical properties of the earthquake catalogue. Some of the source models used in this analysis are taken from published literature (Papazachos, 1990; Papaioannou and Papazachos, 2000), whilst others are developed by original means. In either case, for each source zone MMAX must be calculated via analysis of the earthquake catalogue. When implementing an algorithmic approach to MMAX estimation, there are two logical constraints that limit the range of upper bound magnitudes. The first constraint is that MMAX cannot be less than the maximum observed earthquake in the catalogue (allowing for some uncertainty in historical earthquake magnitudes). The second constraint is that upper bound magnitude estimates must be consistent with both the fundamental limits of the breaking strength of rock (i.e. not unrealistically large Mw N 9.5), and the geodynamics of the zone (e.g. Mw N 9 in an extensional environment). For rapid estimation of MMAX, the cumulative moment method is used. This is an adaptation of the cumulative strain energy release method (Makropoulos and Burton, 1983). This method has several advantages over alternative catalogue-based estimators of MMAX. The main advantage is that it is non-parametric, and therefore does not require the fit of a probability density function to earthquake magnitude recurrence. This avoids propagation of errors from the fit of the probability density function to the estimation of maximum magnitude. Cumulative moment is also influenced much more strongly by large events than smaller events, and is therefore insensitive to small variation in completeness magnitude. This means that the recurrence of larger earthquakes in the historical record will influence the estimation of MMAX, rather than the relatively short duration instrumental record. An alternative estimator of MMAX is the non-parametric Gaussian method of Kijko (2004). This method

offers many of the same advantages as cumulative moment, albeit with more computational complexity. Furthermore the value of the smoothing factor in the Kernel function has a minor, but non-trivial, influence on MMAX (Kijko et al., 2001). This, and the computational speed of the cumulative moment method, is the primary reason why cumulative moment is preferred here. Uncertainty on the estimation of MMAX is difficult to determine when non-parametric methods are invoked. The uncertainty derives from two sources (Kijko, 2004): the extent to which the cycle of strain accumulation and release is captured by the historical catalogue, and erroneous determination of magnitude value. In the case of the Aegean, there is insufficient information to quantify the magnitude uncertainty terms for non-parametric estimation of MMAX; hence only the best estimates are considered. A regional maximum magnitude has been estimated using cumulative moment and can be seen in Fig. 2 (cumulative moment plot). Using the whole Aegean earthquake catalogue a maximum magnitude of Mw 8.6 is determined, with a corresponding maximum repose time (TMAX) of approximately 1000 years. Seismic moment is calculated using the formula of Hanks and Kanamori (1979). 3. Modelling the seismic source 3.1. Issues in the development of source models The development of a seismic source model can be one of the most subjective and controversial elements of seismic hazard analysis. There exists little consensus on how best to approach the task, and what weighting particular types of data should be given. This problem is complicated further when addressing the issue of data quality, particularly with regard to earthquake catalogues. The uncertainty inherent in seismic source zonation is far from trivial when applied to probabilistic seismic hazard analysis (PSHA). When producing a seismic source model there is often no formal methodology to guide the delineation of seismic sources. It is common for source models to be constructed by a panel of experts, who take into consideration geology, seismotectonics and seismicity, each weighted according to their own judgement. The results are often subjective,

256

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

Fig. 2. Regional maximum earthquake magnitude determined via the cumulative moment method of Makropoulos and Burton (1983).

representing a compromise of expert opinion rather than each expert's best judgement. Alongside the issue of subjectivity in the development of the seismic source model, is the lack of objectivity in assessing the quality of the model. This is frequently done by expert judgement. Quantitative approaches to assessment of a source model's quality have been suggested (Musson, 2004), where judgement is usually based on the ability of a model to replicate observed seismicity. Even these quantitative approaches may contain bias, depending on the relative role that seismicity has played in the model's delineation. If observed seismicity primarily guides the delineation of the seismic source, be it point, line or zone, then this source model will naturally provide a better fit to observed seismicity than a model derived using alternative geological and seismotectonic information. The application of the Monte Carlo method to PSHA allows for more flexibility in the manner in which the seismic source is characterised. Previous implementations have simulated the synthetic seismicity via random sampling (with replacement) of observed epicentres (Ebel and Kafka, 1999), random sampling of points within uniform zones (Musson, 1999a; Wiemer et al., 2009) or even random simulation of fault planes within a zone (Sinadinovski et al., 2005). When considering seismic hazard across a region on the scale of a country or group of countries it is possible, even likely, that a uniform source model for the region will span areas of high-seismicity where active sources lie in close proximity (such as the Ionian islands or Gulf of Corinth) and are well constrained by observed geological, geodetic and seismotectonic information. They may also span areas of low to moderate intraplate seismicity, where large earthquakes have longer return periods and the locations of active sources may be harder to constrain due to rapid sedimentation or erosion of faulting features. In such areas the assumption of spatial uniformity is not well-founded. It may therefore be necessary to use observed seismicity as a basis for simulating synthetic hypocentres, in which case the random resampling approach of Ebel and Kafka (1999) may be appropriate. The simultaneous application of both uniform zones and distributed seismicity allows for yet more flexibility in the process of source delineation. Uniform sources can be implemented where seismotectonic and geological information may justify such delineation. Elsewhere, rather than ascribe a uniform or a “background” zone, the Monte Carlo simulations can re-sample observed seismicity, taking the uncertainty in source location into consideration by the use of Gaussian distributed scatter. This approach represents a compromise between guiding hazard according to historical probabilism, where hazard is defined on the basis of previous earthquakes, and seismotectonic

probabilism, which incorporates geological and seismotectonic information (Muir-Wood, 1993). Whilst the latter approach may be considered by some the ultimate objective of seismic hazard analysis, for many regions the geological and seismotectonic information is not yet sufficient to derive reliable hazard estimates. The most widely used source model in the Aegean region is that of Papaioannou and Papazachos (2000), hereafter referred to as PP2000. This model separates the region into 67 shallow (Fig. 3a — PP2000) and 7 intermediate-depth zones. This model is derived by expert judgement, based on the interpretation of seismological, geological and geodetic data. The decision making process is not well-described in the published literature for this model (Papazachos and Papazachou, 1997; Papaioannou and Papazachos, 2000). There are some consistencies with previous models, especially around the Adriatic–Ionian–Hellenic margin. In other regions such as northern Greece or the southern Balkans the justification for the size and shape of many of the zones is not clear. Relevant parameters for seismic hazard analysis (a and b value, maximum magnitude, major axis of isoseismals etc.) are also given in PP2000. However, these parameters have been determined using a kernel smoothing method (Papazachos, 1999b), which is designed to reduce the variation in a- and b-value arising from transient features of the localised portion of the earthquake catalogue within each zone. Prior to the PP2000 model, the source model of Papazachos (1990) was widely used, hereafter referred to as PZ1990. This model contains 36 uniform shallow seismic source zones (Fig. 3b — PZ1990) and 5 intermediate-depth zones. The obvious difference between the two models is in the areas ascribed as “background”. PZ1990 leaves much of the southern Aegean Sea and northern Greece as a “background” zone, with no obvious control on epicentres and with a b-value of 0.8. The delineation of low seismicity uniform zones in PP2000 was a response to the shortcomings of the PZ1990 in leaving seismically active areas as “background” seismicity. As has been suggested here, the Monte Carlo method of PSHA could use the PZ1990 model in a different manner. Here, rather than consider “background” seismicity to be uniform or trivial, seismicity outside the uniform source zones can be incorporated into the simulations by random sampling, with replacement, of the observed seismicity within the “background” region. This constitutes a hybrid approach whereby both uniform zones and distributed seismicity are considered concurrently. For this reason, the PZ1990 model shall be compared in this analysis with the other seismic source models of the Aegean region. 3.2. Seismic source delineation via K-means cluster analysis Given the inherent subjectivity in the delineation of uniform seismic source zones, Weatherill and Burton (2009) presented a method of delineation using K-means cluster analysis. This automated patternrecognition approach partitions the spatial distribution of, separately, hypocentres as point sources, and fault ruptures as line sources, into a given number of clusters (K). The optimum number of clusters is defined according to partition quality indices and fit to observed seismicity. The spatial distribution of the clusters then directs the delineation of source zones, with seismicity parameters defined according to the subset of earthquakes within each cluster. Attempts to quantitatively validate the partitions, by way of numerical partition quality indices, and the resulting seismic source models did not identify a unique “ideal” partition or zone model. Instead those partitions that resulted in “better” scores were appraised in the context of Aegean seismotectonics. The K-means approach clearly represents a different perspective on the issue of seismic source delineation. There are still many issues that need to be resolved, not least of all the transition from partitioned clusters to seismic source zones. The tessellation approach implemented by Weatherill and Burton (2009) was designed to illustrate how this can be done in an algorithmic fashion. This allows for identification of the optimum number of clusters (a fundamental point of uncertainty in cluster analysis) in the context of the application to seismic hazard.

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

257

Fig. 3. Papazachos and Papaioannou (2000) source model (a) and Papazachos (1990) source model (b).

Appraisal of the partitions in this context is a useful tool, given a welldistributed set of data such as an earthquake catalogue. It allows for source delineation to be treated as an inverse problem, one in which the improvement in partition fit that accompanies higher numbers of clusters is balanced against the increased uncertainty that comes with delineating zones with few earthquakes. An alternative to tessellation where zones are delineated according to the fit of a spatial Gaussian function has also been explored (Weatherill and Burton, 2009). However, the increased complexity and uncertainty that resulted from fitting the spatial function did not drastically alter the delineation of zones when compared with the simpler tessellation function. For that reason the tessellation method is preferred here, but alternative approaches could have many uses when applying this method elsewhere. Two of the K-means partitions will be used in this work: the 27 zone rupture partition and the 29 zone rupture partition, shown in Fig. 4a and b respectively. These particular delineations provided the best fit to observed seismicity across the Aegean. Appraisal of these partitions in the context of Aegean seismotectonics also showed that they partitioned seismicity appropriately given the existing knowledge of Aegean tectonics. The fit of these two particular models is very similar; hence they are afforded equal weighting in the analyses of epistemic uncertainty. 3.3. A hybrid seismic source model for the Aegean A further novel source model is introduced here that is designed as a hybrid of uniform zones and distributed seismicity. Uniform source

zones are delineated where we believe there is sufficient information from seismotectonics, geology and seismicity to do so. Elsewhere the observed seismicity is used to characterise the seismic source. This model represents a subjective interpretation of Aegean seismotectonics resulting from consultation of the following data: GPS data and geodynamic models (McClusky et al., 2000; Nyst and Thatcher, 2004; Reilinger et al., 2006), observed seismicity (from existing catalogue), focal mechanisms catalogues (Kiratzi and Louvari, 2003; Vannucci et al., 2004) and published literature (various sources cited where appropriate). The source model, hereafter referred to as WT2006, is shown in Fig. 5. A region-by-region explanation of the decisions made in delineation can be found in the Electronic Supplementary Material to this paper (Appendix B).

4. Selection of strong ground motion attenuation models for use in the Aegean 4.1. Selection of PGA attenuation models The translation of earthquake recurrence interval into strong motion return period via predictive strong motion attenuation relations is another complex issue in PSHA. Regardless of whether it is the Cornell (1968)–McGuire (1976) method, extreme value analysis or Monte Carlo methods, empirical attenuation relations are usually needed to describe the decay of ground motion with distance from an earthquake of a given magnitude. The selection of candidate

Fig. 4. K-means derived zones for the Aegean region. (a) 27-zone model and (b) 29-zone model (Weatherill and Burton, 2009).

258

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

Fig. 5. WT2006 Hybrid Source model for the Aegean region (see Electronic Supplementary Material for a detailed description).

attenuation relations is not a trivial procedure, and can often be oversimplified (or at least poorly justified). Strong motion attenuation in the Aegean region has been an important focus of seismic hazard research. Many PGA attenuation models have been published in the last two decades that are derived from exclusively Greek strong motion records (Theodulidis and Papazachos, 1992; Margaris et al., 2001; Skarlatoudis et al., 2003; Danciu and Tselentis, 2007). However, Greek and Turkish strong motion records also form a significant portion of the strong motion data set used in the derivation of European and Middle Eastern PGA and spectral acceleration attenuation (e.g. Ambraseys, 1995; Ambraseys et al., 1996, 2005; Bommer et al., 2007). The prevalence of Aegean records in the European Strong Motion Database (Ambraseys et al., 2004) indicates that European PGA attenuation models may reasonably be applied to the Aegean in addition to the Greek models. The recent publication of Next Generation Attenuation (NGA) models has highlighted many interesting issues in the selection of empirical attenuation models. NGA models are derived from a standardised database (Chiou et al., 2008) of strong motion records collected from across the globe. Although there is a geographical bias towards Californian earthquakes, and a substantial number of records from the 1999 Chi-Chi, Taiwan, earthquake, Greek and Turkish strong motion records are included in the data set. Although the NGA dataset contains earthquakes from around the globe, it cannot be immediately assumed that NGA relations are equally applicable across the globe. Discussion of the application of NGA relations in Europe can be found in Campbell and Bozorgnia (2006) and Stafford et al. (2008). Some new issues emerge when considering the use of these relations for seismic hazard mapping rather than site-specific studies. Many of the NGA relations introduce new parameters into the attenuation model (such as footwall/hanging wall terms or depth to the 2.5 km s−1 shearwave velocity horizon). These parameters may not be well-constrained across an area, at least not on the resolution required for hazard mapping, often leading to PGA hazard maps based on fixed or assumed values for these terms. Site condition is one such example of a parameter that is not easily constrained in the Aegean. For many of the strong motion records in the European and Greek databases, site investigations

were not nearly as extensive as those for other regions (e.g. Japan, California). As such, site condition may often be characterised by NEHRP site class (as is found in Greek or European attenuation models), or Vs30 for some sites. Of the five NGA models presented, the Boore and Atkinson (2007) model will be considered in this application. This relation depends on the fewest number of parameters (Mw, JoynerBoore Distance, Vs30, Fault Type, aleatory uncertainty), which allows for the best comparison with other attenuation models. It should not be assumed that alternative NGA models could not be applied to sitespecific hazard analyses in the Aegean region. The selection of ground motion models for use in the Aegean region is done using the procedure suggested by Cotton et al. (2006), also favoured by Tselentis et al. (2010). Candidate models are initially rejected on the basis of five criteria: i) the model is clearly from an irrelevant tectonic regime, ii) the model is not published in an international peer-reviewed journal, iii) documentation of the model and the underlying dataset is insufficient, iv) frequency range and the functional form are inappropriate and v) the regression method and regression coefficients are inappropriate. Additional candidates such as Theodulidis and Papazachos (1992) and Makropoulos and Burton (1985b) were immediately rejected on the basis of erroneous characterisation of site (in the former) and inadequate modelling of aleatory uncertainty (in the case of the latter). Unlike Tselentis et al. (2010) and Mantyniemi et al. (2004), we also reject Margaris et al. (2001) as it does not appear in peer reviewed literature, instead Skarlatoudis et al. (2003) is preferred in its place. Six attenuation models passed the search criteria for applicability to the Aegean region. These are Ambraseys et al. (1996) (hereafter cited as Am96), Skarlatoudis et al. (2003) [Sk03], Ambraseys et al. (2005) [Am05], Danciu and Tselentis (2007) [DT07], Bommer et al. (2007) [Bm07] and Boore and Atkinson (2007) [BA07]. Ambraseys et al. (1996) is preferred in place of Ambraseys (1995). Key properties of these attenuation models are given in Table 2 and the relations given in Appendix A. A comparison of the attenuation models for Mw 5, 6, 7 and 8 is shown in Fig. 6. There are many differences in the attenuation models that complicate direct comparison. The most notable is the difference in the components of horizontal motion used. Both Am96 and Am05 use

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

the larger component of horizontal motion as their regression variable. DT07 and Bm07 use the mean (arithmetic and geometric, respectively) of the horizontal components whilst Sk03 uses each horizontal component separately. BA07 uses a slightly different parameter in the form of orientation-independent geometric mean (Boore et al., 2006), although the scaling between this and the more conventional geometric mean is close to unity (Beyer and Bommer, 2006). All the models except Am96 use Mw as the magnitude parameter. Eq. (2) is used to convert from Mw to MS in the synthetic earthquake catalogues. The greatest distinction between the “Greek” relations (Sk03 and DT07) and the others is in their use of epicentral distance (REPI) rather than Joyner–Boore distance (RJB) (see Appendix A). To convert to Joyner–Boore distance, empirical relations presented in Montaldo et al. (2005) are compared: RJB = −3:5525 + 0:8845REPI

Ambraseys and Bommer (1991)

ð4Þ

RJB = −5:0497 + 0:9433REPI

Ambraseys et al. (2004)

ð5Þ

RJB =

h i  2 0:566−0:114MS 2 1 =2 REPI + 25 10 −5:8

Montaldo et al. (2005) ð6Þ

Eq. (10) requires conversion of magnitude to MS, which adds an additional uncertainty into the calculation and as such will not be used in the final calculations. Characterisation of site and fault condition is also inconsistent across these relations. All except BA07 use NEHRP site class as the site variable in the regression. In all cases only three classes are considered: Rock (NEHRP B), stiff soil (NEHRP C) and soft soil (NEHRP D). BA07 use Vs30 as a direct parameter. For comparison with the other relations an “average” Vs30 value is estimated for each soil type: Vs30 = 760 m s−1 (Rock), 570 m s−1 (stiff soil), 300 m s−1 (soft soil). There is less consistency in the definition of fault type. The “Greek” relations consider only two fault classes: normal and thrust/strike-slip. No fault term is included in Am96, but four conditions are included in Am05: normal, strike-slip, thrust and oblique. Bm07 make a distinction normal, strike-slip and thrust faulting but not oblique faulting. BA07 considers three fault conditions (normal, thrust and strike-slip) and an “unknown” mechanism. Differences in the classification of each fault type within the attenuation models also make direct comparison more difficult (Bommer et al., 2003). 4.2. Quantifying the fit of attenuation models to strong motion data The maximum likelihood procedure of Scherbaum et al. (2004) is used to assess the quality of fit of the relations to observed strong

259

motion data from the Aegean. The observed strong motion data is a subset of the European Strong Motion Database (Ambraseys et al., 2004), corresponding exclusively to records from Greece, Albania, FYROM and western Turkey. This provides 146 records from 37 different earthquakes. Some of these records come from earthquakes whose magnitudes exceed the validity range for some equations, such as those from the 1999 Izmit (Mw 7.6) and Duzce (Mw 7.2) earthquakes. Although extrapolation of attenuation models to larger earthquakes is not recommended, the magnitudes of the Izmit and Duzce events are less than or equal to the maximum magnitudes for some parts of the Aegean. Since large events will be simulated in the synthetic catalogue it is important to see how well “Greek” attenuation relations compare with observed attenuation for large events. Where conversion of the distance metric is required, the impact of the empirical Eqs. 4–6 will be considered separately. The classification scheme suggested by Scherbaum et al. (2004) is also applied here, as shown in Table 3. For horizontal PGA, the BA07 model provides the best fit when the Ambraseys and Bommer (1991) conversion (Eq. (4)) from epicentral to Joyner–Boore distance relation is used. Hence this distance conversion relation is adopted within the hazard analysis. DT07 and Sk03 both fall within category B, although the absolute mean and median of the normalised residuals is smaller for DT07. Both the Bm07 and Am05 models are classified as category C, despite having good likelihood values. This low ranking arises due to the large standard deviation of the normalised residuals. It is not clear exactly why this is the case, although it should be noted that of the relations only Am05 and Bm07 contain a heteroscedastic aleatory variability term. The influence of this term, especially with respect to strong ground motion records from smaller earthquakes, may account for the greater dispersion of residuals. However, the Am96 model doesn't have a heteroscedastic aleatory variability term and provides a similarly poor fit. The emergence of BA07 as a viable attenuation relation for the PGA hazard in the Aegean is an interesting development. This result is in contrast to that of Stafford et al. (2008), who find that although BA07 generally provides a good fit to European strong motion data, the Am05 relation should be preferred. It should be noted that in their analysis, Stafford et al. (2008) implement a likelihood fitting approach that considers intra- and inter-event aleatory uncertainty separately. The absence of separate intra- and inter-event aleatory uncertainty terms in the Am96 model prevents us from making the same comparison here. From both the quantitative and qualitative analyses it is possible to make a judgement as to which attenuation relations should be implemented in a seismic hazard analysis. All six relations considered here will be compared by the creation of separate hazard maps and curves. For the purposes of analysing epistemic uncertainty (be it via logic tree or Monte Carlo methods) both the Am96 and Sk03 relations

Table 2 Key properties of the six PGA attenuation models used in the hazard analysis. Horizontal component

M

MLOWER MUPPER R

RLOWER RUPPER Site (km) (km) class.

Fault

σ

2-stage

Larger

MS

4

225

2-stage

Separate

595

135

Random effects Larger

335

151

Europe & 997 Middle East Worldwide 3552

Attenuation model

Region

No. of No. of Regression records earthquakes

Ambraseys et al. (1996) Skarlatoudis et al. (2003) Ambraseys et al. (2005) Danciu and Tselentis (2007) Bommer et al. (2007) Boore and Atkinson 2007)

Europe & 1260 Middle East Greece 1000

619

Europe & Middle East Greece

7.3

RJB

1

200

B, C, D

None

Single

Mw 4.5

7

REPI 5

160

B, C, D

2 (N, SS/T)

Single

Mw 5

7.6

RJB

0

100

B, C, D

Random effects Arithmetic mean Mw 4.5

6.9

REPI 0

136

B, C, D

4 (SS, N, R, O) Two – Mw dependent 2 (N, SS/R) Two

289

Random effects Geometric mean Mw 3

7.6

RJB

0

100

B, C, D

175

Random effects

8

RJB

0

200

Vs30

Orientationindependent geometric mean

Mw 5

Two – Mw dependent 4 (U, N, SS, R) Two 3 (N, SS, R)

260

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

Fig. 6. Comparison of the six PGA attenuation models for a “rock” site with strike-slip faulting assumed. a) Mw 5, b) Mw 6 and c) Mw 7. Skarlatoudis et al. (2003) and Danciu and Tselentis (2007) relations extrapolated to Mw 7.

are excluded from these analyses. When considering horizontal PGA, both the Am05 and Bm07 relations are preferred over Am96, with the Bm07 relation given a stronger weighting than Am05 where applicable. Similarly, DT07 is preferred as a representative Greek relation to Sk03 due to the improved fit to strong motion data. Weighting schemes for epistemic uncertainty analysis will be discussed in Section 6.

Table 3a Fit of attenuation models to strong motion data from the Aegean region using the classification scheme (Scherbaum et al., 2004) in Table 3b. E(Z), Z50 and Zσ are the mean, median and standard deviation of the observed ground motions normalised via the attenuation model. LH50 is the median likelihood value. Attenuation Relation

Magnitude

Distance

E(Z)

Z50

Zσ

LH50

Class

Am96 Am96 Am96 Sk03 Am05 Am05 Am05 DT07 Bm07 Bm07 Bm07 BA07 BA07 BA07

MS MS MS Mw Mw Mw Mw Mw Mw Mw Mw Mw Mw Mw

RjbAB RjbAmb RjbMon Repi RjbAB RjbAmb RjbMon Repi RjbAB RjbAmb RjbMon RjbAB RjbAmb RjbMon

−0.48 −0.47 −0.33 0.14 −0.62 −0.62 −0.47 0.05 −0.32 −0.32 −0.19 0.05 0.05 0.21

−0.42 −0.45 −0.26 0.17 −0.46 −0.46 −0.25 0.05 −0.23 −0.23 −0.05 −0.02 0.01 0.17

1.30 1.30 1.35 1.20 1.44 1.44 1.49 1.22 1.42 1.42 1.47 1.12 1.13 1.16

0.35 0.37 0.35 0.39 0.41 0.42 0.42 0.38 0.43 0.43 0.40 0.53 0.53 0.53

C C C B C C C B C C C A B B

5. Seismic hazard maps for the Aegean region 5.1. Producing the hazard maps The previous discussion of the source, recurrence and strong ground motion attenuation models begins to illustrate the scale of epistemic uncertainty in seismic hazard analysis in the Aegean. Judgements can be made as to which combination of models are preferred, and will be done so in the next section, but here the focus is simply on the influence of each model. The impact of particular model combinations can be best understood by considering each of the models on the seismic hazard maps and hazard curves for selected sites. We opt to consider seven possible source models and, initially, six attenuation models. The seven possible source models are: PP2000 (with recurrence parameters determined from present catalogue), PP2000 (with recurrence parameters given in Papaioannou and Papazachos, 2000), PZ1990 (recurrence from present catalogue), PZ1990 (recurrence parameters given in Papazachos (1990)), WT2006, K27 and K29.

Table 3b Scherbaum et al. (2004) classification scheme (font style indicates class). Class

LH50

|E(Z)|

Zσ

Validity

A B C D

≥0.4 ≥0.3 ≥0.2 ≤0.2

≤0.25 ≤0.5 ≤0.75 ≥0.75

≤1.125 ≤1.25 ≤1.5 ≥1.5

Excellent fit–applicable to region Good fit–some error Moderate–poor fit Unacceptable–Should not be used

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

Although the source zones may be the same when comparing our recurrence parameters with those of the initial authors in the PP2000 and PZ1990 models, we consider the two scenarios separately. Tselentis and Danciu (2010a) indicate that the recurrence parameters given in PP2000 remain consistent when re-analysed with more recent earthquakes. The absence of uncertainty terms on a- and b- values makes direct comparison between the original and re-analysed models difficult. Since uncertainty is not characterised here, we assume a conservative value of σb = 0.01 and σa = 0.2. In the other models the uncertainty terms on a and b are determined using the Burton et al. (2004a) catalogue. The spatial distribution of seismicity is simulated via a hybrid approach. Epicentres are uniformly distributed within zones where they are defined. Outside of these regions synthetic hypocentres are created by random re-sampling (with replacement) of the observed hypocentres. A two-step sampling process is used. Initially the distribution of the seismicity is sampled as described. The final synthetic epicentre is then sampled from a 2D Gaussian function using the initial synthetic epicentre as the mean position. The standard deviation σx = σy corresponds to the degree of smoothing, which is effectively equivalent to the smoothing term in the spatially smoothed seismicity approach of Frankel (1995). This term is used to take into consideration uncertainty on the observed hypocentre location and to soften the boundaries of zones to reduce the ridge effect of the zone models (Musson and Henni, 2001). In these analyses σx = σy = 0.1∘ which broadly corresponds to the uncertainty in epicentre location for the early 20th century earthquake catalogue. A minimum magnitude of Mw 5.2 is assumed. There are several reasons this value is chosen. Firstly it corresponds to the completeness magnitude for the whole 1900–2005, allowing for a direct correspondence check between the synthetic seismicity rate and observed seismicity over the longest period. The second reason is that the lower bound of applicability of many of the attenuation models considered here is equal to or below Mw 5.0. It has been suggested that not only should attenuation models not be extrapolated below their applicability range, there is also some question as to whether they provide a good fit to strong motion from earthquakes at the limits of their applicability range (Bommer et al., 2007). Given these conditions Mw 5.2 would seem to be a reasonable choice. Finally, assuming existing macroseismic attenuation models for the Aegean (Papaioannou and Papazachos, 2000), Mw ≈ 5 is the approximate threshold of economic damage (EMS VI) for earthquakes in the extreme near-field. It also corresponds to the approximate damage threshold suggested for Nuclear Power Plants of Cumulative Absolute Velocity (CAV) 160 cm s−1 (O'Hara & Jacobson, 1991), assuming the median CAV relation of either Danciu and Tselentis (2007) or Kramer and Mitchell (2006). The maximum magnitude for each zone is determined using the cumulative moment method for observed seismicity in each zone, except where the models use parameters given by the original authors. This nonparametric approach will ensure that MMAX ≥MOBS MAX, regardless of the location. For many zones, maximum magnitude may be equal to, or only slightly larger than MOBS MAX. To allow for the inclusion of earthquakes with REGIONAL magnitudes in the range MOBS , magnitudes from a subMAX b Mw ≤MMAX sample of earthquakes are randomly sampled from the CDF of the regional Gutenberg and Richter (1944) relation (Eq. (3)). 5.2. Calculating hazard from the synthetic catalogues The six preferred attenuation models are Am96, Sk03, Am95, DT07, Bm07 and BA07. Hazard maps presented initially correspond to the definition of horizontal PGA according to each model; hence no conversion between geometric mean horizontal PGA and larger component horizontal PGA is made at this point, but will be done so later. Aleatory variability is modelled by log normally distributed scatter around the median PGA. To allow for comparison between Am96, Sk03 and other models, the scatter corresponds to the total aleatory variability, and distinction between inter- and intra-event

261

terms is not made. The aleatory variability term is truncated at 3σ, on the guidance of previous studies into the impact of this term and the poor fit of the model to residuals whose scatter exceeds 3σ (Abrahamson, 2000; Bommer et al., 2004; Strasser et al., 2008; Tselentis and Danciu, 2010a). Faulting, where defined, is presumed to be strike-slip or normal/strike slip, according to the attenuation model. Again, the impact of this term will be explored later in the paper. The hazard analysis also assumes a “rock” site (NEHRP Class B with Vs30 = 760 m s−1), which will also be explored in the next section. The Monte Carlo procedure allows for ground motion from any earthquake to be calculated at any point in the study area. In practise this is a computationally expensive procedure that will result in excessive calculation of trivial ground motions. Furthermore, none of the attenuation models are intended for use at source-site distances in excess of 250 km. This value is used as the “hazard radius”, i.e. the radius beyond which earthquakes are not considered to contribute to hazard at a site. This distance also corresponds to the maximum radius of MMI VI damage originating from the assumed regional magnitude earthquake (Mw 8.6), assuming the intensity attenuation model of Papaioannou and Papazachos (2000). Large earthquakes have caused damage over greater source-site distances in other areas of the world (e.g. Mexico City, 1985). This is usually due to non-linear site amplification on poorly consolidated sediments. There is little evidence to suggest that historical earthquakes in the Aegean have produced damage at this sort of distance, hence 250 km is a reasonable compromise. Given this condition, we restrict the spatial extent of the hazard maps by a similar amount. This is to avoid producing erroneously low hazard estimates around the periphery of the study area due to the lack of consideration given to seismic sources outside the source zones. This is particularly true in the Marmara region of Turkey, where large earthquakes further east along the North Anatolian Fault also contribute to hazard in this region, but are not considered here. Hazard maps in Fig. 7 illustrate the impact of varying the source model, and Fig. 8 the impact of varying the attenuation model. PGA (cm s−2) with a 10% probability of being exceeded in 50 years is used as the preferred hazard parameter. 5.3. Influence of the source model on hazard maps It is immediately clear that the choice of source model heavily influences the appearance of the hazard maps (Fig. 7). All of the maps identify the Ionian Islands and Gulf of Corinth as areas of high seismic hazard, and the south Aegean Sea and Xanthi region of northeastern Greece as areas of lower seismic hazard. Most of the maps agree on the lower seismic hazard of FYROM, although the K27 and K29 models produce higher hazard in the west of the country due to the active source zones along the Adriatic coast. Disparities in the hazard owing to different source models are more apparent in the Hellenic arc. All the source models recognise the high hazard along the Hellenic Arc, yet the hazard presented in the maps also reflects the delineation of source zones. The PP2000 and WT2006 models both produce higher hazard along the Cretan and Dodecanese sections of the arc than in the southern Greece section. When PZ1990 is used, the contrast between the western section of the arc and the Cretan/Dodecanese sections is not as great. The distinction is also less obvious in the K-means models. For western Turkey, both the PP2000 and PZ1990 models delineate seismic zones in the form of large rectangles. Whilst hazard is relatively high in this region this masks areas of localised high seismic hazard along active fault systems. The WT2006 model best captures this variability as it generates synthetic epicentres by re-sampling the observed seismicity. The parameters of earthquake frequency magnitude do not have as great an influence on the seismic hazard maps as source and attenuation models. For the PP2000 and PZ1990 source model we consider two sets of magnitude frequency parameters: “automatic”

262

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

263

Fig. 8. Influence of the attenuation model on PGA (cm s−2) with a 10% probability of being exceeded in 50 years, assuming a “rock” site (Vs30 760 m s−1) and normal faulting using the PP2000 source model. Influence of the attenuation model: a) Am96, b) Sk03, c) Am05, d) DT07, e) Bm07, f) BA07 attenuation models. PGA represents the horizontal motion as defined according to the original publications (Table 2). No conversion is made.

parameters are calculated for each zone using the input catalogue described in Section 2, “manual” parameters for each model are those given in the literature. For the WT2006, K27 and K29 models, earthquake magnitude frequency parameters are determined for each zone using the current catalogue. Although small differences are apparent, the overall impact on seismic hazard is low. This may illustrate the comparative dominance of source and attenuation model uncertainty in the seismic hazard analysis. However, as

recurrence parameters are sampled from a probability distribution with a given uncertainty, then differences between assigned and automatically calculated recurrence parameters may be smaller. 5.4. Influence of the attenuation model Fig. 7 shows how the choice of source model affects the spatial variation in hazard across a region. The actual level of PGA hazard is,

Fig. 7. Influence of the source model on PGA (cm s−2) with a 10% probability of being exceeded in 50 years, assuming a “rock” site (Vs30 760 m s−1) and normal faulting using the BA07 attenuation model and different source models: a) PP2000 (automatic), b) PP2000 (manual), c) PZ1990 (automatic), d) PZ1990 (manual), e) WT2006, f) K27 and g) K29.

264

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

Fig. 9. Formulation of the logic tree framework implemented in the Monte Carlo Analysis of Epistemic Uncertainty (MCMAEU).

however, dependent on the choice of attenuation model as illustrated in Fig. 8. Two trends are obvious. The “Greek” attenuation models (Sk03 and DT07) produce hazard maps with lower PGA across the entire region. Conversely, Am05 produces seismic hazard maps with considerably higher PGA than other relations. The remaining relations Am96, Bm07 and BA07 show some differences, but there is greater agreement between the hazard maps when these models are used. The comparison of the PGA attenuation models in Fig. 6 does provide some insight into why the models influence hazard in the manner that they appear to. There is some convergence between the six attenuation models at Joyner–Boore distances greater than 10 km, when considering magnitudes in the range 6 b Mw b 7. The greatest disparities between the models are found in the near-field (RJB b 10 km). For smaller earthquakes both DT07 and Sk03 predict lower PGA than for the other models, whilst Am05 predicts much higher PGA to a greater source-site distance. This suggests that near-field seismicity has a great influence on seismic hazard across the region. As near-field strong motion records are not as abundant in Europe, and show great variability due to directivity, ground motion in this region is not as well-constrained by existing attenuation models as it is in the intermediate- and far-field region. Another key issue to consider is the impact of the aleatory uncertainty term of the attenuation model. Both Bm07 and Am05 contain a heteroscedastic uncertainty term, which shows decreasing sigma with increasing magnitude. Noting briefly that the minimum magnitude used here is within the applicable magnitude range for all the attenuation models, the influence of this correlation should permeate into the hazard analysis. This may also account, at least partially, for the high PGA hazard when the Am05 relation is used. Musson (2009) has demonstrated the influence of the Am05 model on idealised hazard scenarios, which illustrates the insensitivity of the hazard curve to changes in maximum magnitude. However, the similarity in hazard between the Bm07 and BA07 relations would suggest that the size of the aleatory uncertainty term is not as influential. This is especially true when taking into account that at the minimum magnitude (Mw 5.2) the total aleatory uncertainty in the Bm07 relation is ≈ 0.34, which is much smaller than the fixed sigma in BA07 of 0.5. 6. Quantitative analysis of epistemic uncertainty 6.1. Characterising epistemic uncertainty The hazard maps shown in Figs. 7 and 8 demonstrate the influence that different combinations of source and attenuation model have on the seismic hazard across the Aegean. No judgement has yet been made as to which combination of combinations should be preferred for the “output” of the seismic hazard analysis. There are many

qualitative criteria upon which a judgement could be based (e.g. the most/least conservative map, the map that is most consistent with previous estimates, the map derived from exclusively Greek models etc). Many of these criteria are subjective or, if using previous hazard estimates as a basis for validation, are based on circular reasoning. The dilemma in selecting particular seismic hazard models comes from the shortcomings in seismic hazard validation. Quantitative means are available to assess the fit of the source model to observed seismicity (Musson, 2004), and the fit of the attenuation model to observed strong motion (Scherbaum et al., 2004). Validation of the whole hazard result is a more complicated process, and is one that is hugely limited by available information. Attempts to constrain hazard estimates using strong motion records (Ordaz and Reyes, 1999; Beauval et al., 2008) are certainly promising, but validation can only really be limited to a small number of sites, and for a relatively short time period, due to the limitations of the strong motion records. For seismic hazard over a “geological timescale” precariously balanced rocks (PBRs) have been shown to provide upper bound constraints to the seismic hazard (Brune, 1996; Anooshehpoor et al., 2004; Stirling and Anooshehpoor, 2006). This method too is limited to a small number of sites, and to our knowledge no PBR studies have been published for the Aegean region. Macroseismic intensity records may offer the best means of hazard validation in terms of spatial coverage and historical duration. However, correlation between macroseismic intensity and engineering based strong motion parameters show considerable scatter (Tselentis and Danciu, 2008). Uncertainty in validating data on such a scale as this would likely result in a substantial degree of non-uniqueness in the model fit, making it hard to identify even the better models. Quantitative validation of the whole seismic hazard model should be an important focus for future research. In lieu of this, separate judgements have to be made about which models are preferred. Where several models are applicable it is necessary to combine this into an analysis of epistemic uncertainty. The most common approach to analyse epistemic uncertainty is the use of a logic tree (Coppersmith and Youngs, 1986; Reiter;, 1990; Kramer, 1996). This is a method of formalising the different models and indicating the spread of their respective hazard values. Each node of the tree represents a part of the hazard analysis where different models are considered, whilst the branches of the tree at each node represent the models themselves. In common seismic hazard practise, each branch of the logic tree at a given node should satisfy the criteria of mutual exclusivity and collective exhaustiveness. The mutual exclusivity condition requires that each model [branch] be recognised as independent of other models, such that no model is considered redundant. Collective exhaustiveness assumes that a “true model” (a model that fully

Table 4 Weighting schemes for the (a) attenuation model and (b) the source model. a) Attenuation model

PGA class

EXP

EXP+

LIN

SOC

BA07 DT07 Bm07 Am05

A B C C

0.5 0.25 0.15 0.1

0.7 0.2 0.08 0.02

0.4 0.3 0.2 0.1

0.25 0.25 0.25 0.25

b) Source model

χ2

Weight (χ2)

Weight (SOC)

PP2000 (automatic) PZ1990 (automatic) WT2006 K27 K29 PP2000 (manual) PZ1990 (manual)

1.368 1.264 1.157 1.294 1.282 1.605 1.360

0.06 0.12 0.3 0.2 0.2 0.04 0.08

0.1 0.1 0.2 0.2 0.2 0.1 0.1

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

describes the observations, within the accepted imprecision) can be found within the set of models considered. In practise, there may often be some degree of balance between selecting a sufficient number of models to satisfy the collective exhaustiveness criterion without violating mutual exclusivity (Bommer and Scherbaum, 2008). The issue of model independence within this hazard analysis shall be addressed shortly. Weights are assigned to each branch to represent the relative degree of confidence in each model. Each permutation of models is then given an overall confidence score (the term probability is avoided as the weights are usually assigned by judgement). The formulation of the logic tree used in this analysis is given in Fig. 9. The question of how to extract hazard values from the logic tree still remains an area of debate. The assimilation of competing, weighted, branches into a single distribution is often used, raising the issue of whether to take the mean or median value. Advocates of the median value suggest that the mean value is unfairly weighted by the least likely models when considering very low annual probabilities (Abrahamson and Bommer, 2005). Conversely, it is also argued that fractile levels (e.g. median or 84th percentile) are not truly probabilistic as the probability of the particular branch corresponding

265

to that fractile level is itself low (Musson, 2005). A further question arises when considering the number of branches that should be placed into the logic tree. If a particular branch is apportioned very small weighting (P b 0.05), is there any real merit in including the model at all? In framing a logic tree in the manner undertaken here, decisions have to be made to determine the weighting given to each model within the layer of the logic tree. Six PGA attenuation models have been considered in the previous section. We opt to prune the logic tree by removing two of these models from further consideration. The Am96 and Sk03 models are removed from further analysis. In both cases other models are included that display similar behaviour but have been derived from an updated data set, these are Bm07 and DT07 respectively. As BA07 gave the best fit to the observed PGA data (Section 5) this model will also be included. Am05 represents something of a borderline case for inclusion. The likelihood fit to observed strong motion data is similar to Bm07, yet the hazard maps for which it is used clearly show higher accelerations. Hazard using this model appears to be controlled by smaller near-field earthquakes, even at the relatively high cut-off magnitude of Mw 5.2. We opt to include the Am05 model, albeit with a

Fig. 10. PGA (cm s−2) with a 10% probability of being exceeded in 50 years on a “rock” site, using the MCMAEU method for analysis of epistemic uncertainty. The BA07 attenuation model is used, and the χ2 (a) and SOC (b) zone weighting scheme compared.

266

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

lower weighting than the others. We recognise that the similarity in geographical scope and model functional form between the Am05 and Bm07 models may suggest that the two models cannot be considered independent, thus violating the criterion of mutual exclusivity. However, given the extension of the Bm07 model to lower magnitudes and the resulting change in the input data set, along with the contrasting ground accelerations that each model produces, there remains a case for inclusion of the Am05 model. Four weighting schemes for the four PGA attenuation models are tested. The rank of the models (from best fit to poorest fit) is based on their classification according to the likelihood fit in Section 5 (Scherbaum et al., 2004). Weights are apportioned in accordance with the rank of each model, albeit with a different approach. The four schemes are defined as follows, with the weights indicated in Table 4a: i) Exponential (EXP): best fitting models are strongly weighted ii) Exponential-plus (EXP+): best fitting models are very strongly weighted, iii) Linear (LIN): weighting from best to poorest fitting models decreases in linear increments. iv) Socratic (SOC): all models are assigned equal weighting Clearly expert judgement could be used to add on more weighting schemes (Scherbaum et al., 2004). We shall consider only these schemes as representative of how weighting can be apportioned so that a qualitative judgement can be made as to which models are “better” fitting than others. When comparing models in an epistemic uncertainty analysis the differences in the respective forms of the attenuation models must be taken into account. Firstly, two models define horizontal PGA according

to the geometric and arithmetic mean of the horizontal components (Bm07 and DT07, respectively). BA07 uses the orientation-independent geometric mean (Boore et al., 2006), although scaling between this and geometric mean is shown to be close to unity (Beyer and Bommer, 2006). Am05 uses the larger horizontal component of PGA, and as such is scaled down for comparison with geometric mean using the scaling factors presented in Beyer and Bommer (2006). Assigning weights to the source models is a more complex process. The χ2 process described in Weatherill and Burton (2009) is used to assess how well each of the seismic source models replicates observed seismicity in the Aegean. It compares the distribution of the largest ground motions predicted across a region using the observed earthquake catalogue and the DT07 Arias Intensity attenuation model, with those ground motions determined in the same fashion using multiple synthetic earthquake catalogues. A comparison of the fit of each of the seven models (five source models, but with automatically calculated and author suggested recurrence parameters treated separately) is given in Table 4b. This indicates that WT2006 and PZ1990 hybrid model provide the best fit to the observed seismicity of the Aegean. It should be noted, however, that since the hybrid models sample from observed seismicity in parts of the Aegean, the statistic is biased towards these models. For this reason we are cautious about applying strong weightings that will result in these models dominating the ensuing hazard analysis. As judgements regarding the ranking of each model are harder to make, or at least the case for higher ranking less certain, we implement only two weighting schemes: χ2 (where weighting is apportioned on the basis of the χ2 fit) and SOC (where the five zone models are weighted equally and weightings for the PP2000 and PZ1990 models are evenly split between the two recurrence parameter schemes). These schemes are shown in Table 4b.

Fig. 11. As Fig. 10, using the PP2000 source model, with attenuation model weighted via the a) EXP, b) EXP+, c) LIN and d) SOC weighting schemes.

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

Once again, the issue of mutual exclusivity is raised in opting to include both the K27 and K29 model within the epistemic uncertainty analysis. We recognise that these models may also violate the criterion of mutual exclusivity, although true independence of any pair of source models is difficult to demonstrate. There are no zones that are identical in terms of spatial extent or earthquake recurrence in both the K-means models. Whilst the data set input into the cluster analysis is the same, the delineation of the zones for each cluster analysis is not influenced by the delineation resulting from other values of K. It is inevitable that similarities in the models will exist, but independence between two partitions can be neither proven nor disproven. We opt to consider multiple K-means partitions as there was little basis, from both a numerical quantification and a seismotectonic interpretation, for selecting one model in clear preference to the other. In this respect it may be necessary to reiterate that we apportion model weightings on a relatively subjective basis that reflect the degree of confidence in each model notionally quantified using the methods described previously. It is therefore more appropriate to consider the weightings to be subjective indications of relative merit (Musson, 2005), in which case complete independence is neither a relevant nor realistic expectation. 6.2. Monte Carlo Method for the Analysis of Epistemic Uncertainty (MCMAEU) The use of Monte Carlo techniques to explore epistemic uncertainty within the hazard analysis is well established (Cramer et al., 1996; Smith, 2003). Here we develop this further by randomly sampling pathways through the logic tree, within the seismic hazard analysis itself. This is done by dividing the Td-year synthetic catalogue into Ns

267

sub-catalogues of equal duration. For each of the sub-catalogues the source and attenuation model are then randomly selected from the logic tree using a roulette wheel selection method (Coley, 2005). This means that the likelihood of selecting a given model is proportional to the weight of model within the logic tree framework. The properties of the synthetic earthquakes in each sub-catalogue correspond to those of the source model selected. Likewise, the ground motions arising from the earthquakes in each sub-catalogue are modelled using the attenuation relation selected for the sub-catalogue. The resulting full synthetic catalogue, from which the hazard is determined, contains all NS sub-catalogues. As regional seismicity is used to anchor overall seismicity rate, the inter-event time characteristics of the synthetic catalogues are deemed irrespective of source model. However, the relative proportion of earthquakes simulated within each zone is controlled by the seismicity rate for the zone, as given within the model selected for the sub-catalogue. The Monte Carlo method of analysis of epistemic uncertainty (MCMAEU) is subject to many of the same considerations that are required for the conventional application of a logic tree. Where it differs is in the use of Monte Carlo sampling to integrate the epistemic uncertainty within the hazard calculation, rather than producing separate hazard analyses for each logic tree path. The output of the MCMAEU is similar to that of the mean hazard for the logic tree, albeit that no formal selection of mean or median need be made. This particular approach is unique to the Monte Carlo method of seismic hazard analysis; hence direct comparison cannot necessarily be made between this and logic trees applied to PSHA using the standard Cornell (1968) and McGuire (1976) approach. It may, however, offer an advantage for developing the hazard input required in seismic design codes. It allows for a single hazard map to be produced for

Fig. 12. As Fig. 11, using the χ2 zone model weighting scheme.

268

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

which epistemic uncertainty is already integrated. Design earthquakes can still be identified in the manner described by Musson (1999b), with the additional advantage that it is possible to identify the source model and attenuation model associated with each earthquake producing the expected ground motion at a site. This may present a clearer selection of scenario events than can be determined from disaggregation of the mean hazard, whilst also allowing the end user to identify the relative contribution of each source and attenuation model to the hazard analysis. To allow for a more open scrutiny of the epistemic uncertainty analysis, we consider the impact that different weighting schemes may have on each layer of the logic tree. This is achieved by implementing the MCMAEU method separately for each layer. Initially the attenuation model is fixed and the source model is sampled from the layer within the logic tree (Fig. 10), then vice-versa (Fig. 11). We then look at the impact on the hazard map of apportioning the different source model and attenuation model weighting schemes described in Section 6.1. The PGA hazard maps are generated for 10% probability of being exceeded in 50 years, assuming a “rock” site and normal or normal/strike-slip faulting. These maps are illustrated for discussion as follows: 1) Comparison of χ2 (weighting based on χ2 fit to observed hazard) and SOC source model schemes using the BA07 attenuation model (Fig. 10). 2) Comparison of the EXP, EXP+, LIN and SOC attenuation model weighting schemes using the PP2000 source model (Fig. 11). 3) Comparison of the EXP, EXP+, LIN and SOC attenuation model weighting schemes using a variable source model with χ2 weighting (Fig. 12).

4) Comparison of the EXP, EXP+, LIN and SOC attenuation model weighting schemes using a variable source model with SOC weighting (Fig. 13). The faulting regime is assumed to be normal or normal/strike-slip, depending on how it is defined in the respective attenuation relation. The site is considered to be a “rock” site (NEHRP Class B, Vs30 760 m s−1). Horizontal ground acceleration is assumed to be geometric mean of the two components; hence other models have been adjusted accordingly using the scaling factors of Beyer and Bommer (2006). The most obvious feature when these maps are compared is that the hazard maps display few differences when different weighting schemes are used. In the case of the zone model the differences between the weighting schemes were comparatively small to begin with. There is greater disparity between the maps when considering different attenuation model weighting schemes. Hazard is higher when the SOC scheme is used, which is due to the increased weighting of the Am05 model. Otherwise the differences in hazard are small. These maps illustrate that MCMAEU can give hazard analysis results that are comparable to logic tree methods of epistemic uncertainty analysis. The main question to emerge is whether the MCMAEU method is compatible with the demands of seismic hazard analysis as set out in building codes? The application of MCMAEU to PGA hazard shown here provides results that are compatible with the single model Monte Carlo seismic hazard analysis summarised in Figs. 7 and 8. In that respect design earthquakes can also be identified using the disaggregation method of Musson (1999b). These can be used to produce design spectra if required. It should be noted, that within the disaggregation, particular magnitude–distance pairs are decoupled from specific source models. For compatibility with

Fig. 13. As Fig. 11 using the SOC zone model weighting scheme.

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

building codes it will be necessary for the source zone and attenuation models to be recorded for each earthquake contributing to the specific level of hazard at a site. This may help identify or eliminate models that are not well suited to hazard analysis at a site. 7. Sensitivity analyses and further extensions to the method 7.1. Minimum magnitude The reasons for selecting MMIN = Mw 5.2 have already been given. The limitations of the attenuation models would, with the exception of Bm07, prevent the inclusion of earthquakes with magnitudes as low as Mw 4.0. To illustrate the impact that MMIN may have on the hazard analysis we consider a reduction of MMIN from Mw 5.2 to Mw 4.8. This value is within the magnitude range of DT07 and Bm07, and requires only a small degree of extrapolation beyond the ranges of Am05 and BA07. A comparison of two ensuing hazard maps is shown in Fig. 14. There appears to be very little difference in hazard across the region when MMIN is reduced to Mw 4.8. Across most of the Aegean this has the

269

impact of raising the hazard by a small amount. Several studies that have considered the impact of MMIN on the seismic hazard analysis (Bender and Campbell, 1989; Grünthal and Wahlström, 2001; Beauval and Scotti, 2004), particularly in regions of low to moderate intra-plate seismicity (north Eastern U.S., Germany and France respectively). These studies consistently showed a small increase in the level of hazard when MMIN is reduced, which is more pronounced for PGA and shorter period motion than for longer period motion. The effect of reducing MMIN can be seen in the PGA hazard curves for six selected sites (Fig. 15). These curves agree well with the trends observed in previous studies, although the overall impact on the hazard analysis is small at the 475 year return period considered in this analysis. 7.2. Intermediate-depth earthquakes The previous hazard analyses have considered hazard only resulting from shallow seismicity. For consistency with previous Aegean source models (Papazachos, 1990; Papaioannou and Papazachos, 2000), the shallow source zones have an assumed depth of 60 km. However, the Hellenic Arc is a subduction zone where

Fig. 14. PGA (cm s−2) with a 10% probability of being exceeded in 50 years for a “rock” site, using the PP2000 source model and BA07 attenuation model, with (a)MMIN = Mw 4.8, and (b) the change in hazard due to reduction of MMIN from Mw 5.2 to Mw 4.8.

270

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

intermediate-depth earthquakes (60–100 km) have produced damaging ground shaking in the past. This presents a particular hazard to engineered structures with long resonant periods. There are two significant challenges in assessing the hazard from intermediate depth earthquakes. The first is the uncertainty arising from the absence of or large error in depth estimates for the historical and early instrumental earthquake catalogue. Errors on depths may be on the order of several tens of kilometres, which means there is a significant likelihood that some intermediate-depth earthquakes have been erroneously classified as shallow, and vice-versa. It is unclear whether there is any persistent bias, so no judgement is made as to whether intermediate depth earthquakes are over- or underestimated in the earthquake catalogue. The PZ1990, PP2000 and WT2006 all contain uniform source zones to model seismicity arising from the subducting African lithosphere. The K27 and K29 models don't have source zones, but as with shallow seismicity we opt to implement a

hybrid approach where intermediate-depth earthquakes are generated in synthetic catalogues via random re-sampling from the observed distribution. The second, and most important, challenge arises from the lack of suitable PGA attenuation models for the intermediate-depth earthquakes in the Aegean and a paucity of strong motion records from such events. Konstantinou et al. (2006) analyse the fit of the Youngs et al. (1997) subduction attenuation model to strong motion records from the 2006 Kythera earthquake. They find that the Youngs et al. (1997) relation over-estimates the observed PGA. We opt to import the Atkinson and Boore (2003) subduction zone attenuation model, but it should be noted that no Hellenic strong motion records contribute to the data set of either subduction attenuation model. Whilst a few strong motion records are available from which a quantitative likelihood fit could be attempted, the magnitude and distance range is too small to identify systematic bias in the records.

Fig. 15. Comparison of hazard curves for six selected cities in the Aegean, using the BA07 (blue) and Bm07 (red) attenuation relations, when MMIN is reduced from Mw 5.2 (continuous line) to Mw 4.8 (dashed line).

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

271

Fig. 16. PGA (cm s−2) with a 10% probability of being exceeded in 50 years for a “rock” site, using the MCMAEU method with χ2 zone weighting and EXP attenuation model weighting, including intermediate-depth earthquakes.

To incorporate intermediate depth earthquakes into the hazard analysis, some assumptions have to be made. The PZ1990, PP2000 and WT2006 models all make the same assumption that earthquakes originating at depths of 60–100 km are thrust events from the subducting slab interface. Earthquakes in the deeper zone display a steeper dip, suggesting a decoupling of the subducting slab from the Aegean lithosphere. Seismicity in this part of the Benioff zone is lower, which may be attributed to deformation of the slab. Given this model, we designate earthquakes in the 60–100 km depth region to be interface events, following the definition of Atkinson and Boore (2003), and deeper earthquakes to be intra-slab events. The Atkinson and Boore (2003) relation uses rupture distance as the distance metric rather than Joyner–Boore distance. Given the depths of the shallow source zones, the shortest possible rupture distance for intermediate seismicity is 60 km. Although we recognise that rupture distance is not interchangeable with hypocentral distance, for distances over 60 km we consider this assumption to be fair. Hypocentral distance is therefore used as a proxy for rupture distance. A hazard map incorporating intermediate-depth events is shown in Fig. 16. Despite some localised areas of increased hazard (perhaps

most notably south of the Cyclades islands) the overall impact on the hazard maps is small, and most differences are less than 10 cm s−2. A greater increase in hazard may be expected at longer periods of acceleration, but is not visible in PGA. Give the inherent difficulty in modelling hazard from intermediate-depth earthquakes, a need for improved modelling of attenuation from intermediate-depth earthquakes in the Aegean region is demonstrated here. 7.3. Incorporating site condition into the hazard maps The hazard maps presented to this point have assumed a reference rock site (NEHRP Class B with an assumed Vs30 of 760 m s−1) and with normal or normal/strike-slip faulting assumed, depending on how it is defined in the attenuation model. Whilst it may be common in mapping applications of the seismic hazard to produce separate maps for different site conditions and fault types, it shall be illustrated here how these effects can be incorporated into the hazard analysis. The integration of a high resolution database of site condition, here via the topographically derived Vs30 proxy (Wald and Allen, 2007), allows for seismic hazard to be analysed using many of the same site condition

Fig. 17. Variation in site condition across the Aegean region using the global Vs30 database (Wald & Allen, 2009). (a)Vs30 and (b) corresponding NEHRP site class.

272

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

assumptions that are presently in use in the USGS ShakeMap software. This characterisation of site condition is preferred over characterisations implied from geological maps (e.g. Romeo et al., 2000; Slejko and Rebez, 2002) so as to achieve a degree of compatibility with the ShakeMap procedure. Information on site condition across the Aegean comes from the global Vs30 database (Wald and Allen, 2007). This provides estimates of Vs30 at a resolution of 30-arc seconds across the planet, derived from digital elevation models (DEMs), making the distinction between active tectonic and stable continental regions (of which the former is used here). The Vs30 estimates are derived from empirical correlations of Vs30 with slope, which is determined from SRTM 30 arc-second data. The distribution of Vs30 is shown in Fig. 17, alongside a map of the corresponding NEHRP site class. As computational limitations prevent the calculation of seismic hazard at this resolution, hazard is interpolated from the 0.2∘ × 0.2∘grid onto a 30 arc-second grid. For hazard analyses using a single attenuation model the seismic hazard maps may be readily scaled to take into account the site condition. For analyses of epistemic uncertainty the process is more complicated. Three of the four attenuation models used here (Am05, DT07 and Bm07) only define site condition on the basis of NEHRP site classification, with three site classes considered (rock, stiff soil and soft soil). BA07 uses Vs30 directly, so for comparison with the other relations a “reference” Vs30 type for each NEHRP site class is assigned (Stafford et al., 2008). Although topographic slope shows a good correlation with Vs30, there is still considerable error on the estimates. The resolution of the map and the strength of the correlations are such that simple categorisation into three site classes may be the most realistic depiction of how well site condition can be resolved over such a large spatial scale. No amplification is assumed offshore. A map of PGA with a 10% probability of being exceeded in 50 years using the MCMAEU method with site condition included is shown in Fig. 18. The influence on seismic hazard of some of the larger geomorphological features has become visible. In particular, hazard is raised along some of the major river basins in Western Turkey and in the low lying region of Thessaly. The 30-arc second resolution will clearly not resolve the variation in hazard due to local scale site amplification, such as in an urban area. Even though site condition may be incorporated into a hazard map, this is still no substitute for a detailed microzonation or site investigation for urban risk assess-

ments or engineering design. Finer resolution could be achieved using higher resolution DEMs, but in many areas, particularly urban regions, vegetation and urban sprawl introduce error into satellite topography calculations. Allen and Wald (2009) also note that for resolutions finer than 9 arc-seconds the increased variability in slope calculations reduces the confidence in the correlation of slope with Vs30. 7.4. Incorporating fault type into the hazard maps It has become well-established that thrust faulting earthquakes are generally associated with higher accelerations than normal or strikeslip faults (Boore et al., 1997; Campbell, 1997). In the subsequent ten years since publication of these models, the influence of fault type on ground motion became a feature of many strong motion attenuation models. The influence of fault type should also be reflected in the seismic hazard analysis. For site specific analyses where the source may be characterised by a finite number of well-mapped faults, estimation of fault parameters may be done with confidence. When considering seismogenic source zones across a region for hazard mapping applications, more crude assumptions have to be made on the basis of the analyst's interpretation of the seismotectonics. For the purpose of comparison between the four attenuation models, we consider four types of fault class: normal, strike-slip, thrust and mixed. Fault types are assigned to each zone in the PP2000, PZ1990 and WT2006 models. Our assignment for PP2000 differs slightly from that of Tselentis and Danciu (2010a), who consider only two classes: normal and strike-slip/thrust. The inclusion of style of faulting into the hazard analysis when considering epistemic uncertainty introduces additional complication. This arises due to the differences in fault characterisation for each model. DT07 consider only two categories, Am05 considers four types (normal, strike-slip, thrust and oblique), Bm07 considers three types (normal, strike-slip and reverse) and BA07 three types plus an “unknown” category. Furthermore, the classification of fault type on the basis of focal mechanism, and rake in particular, is not consistent between relations (Bommer et al., 2003). Although some adjustment schemes have been suggested (Bommer et al., 2003) there is no particular scheme that can homogenise the style-of-faulting term across these four models. In the hazard maps shown here we opt to take a different approach than simple classification of a zone as displaying only one fault type. We define a “characteristic” set of fault plane parameters (strike, dip and

Fig. 18. PGA (cm s−2) with a 10% probability of being exceeded in 50 years incorporating site condition. MCMAEU method is used with χ2 zone model weighting and EXP attenuation model weighting. Faulting is assumed to be strike-slip.

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

rake) for each zone, and some associated uncertainty. For existing models such as PP2000 characteristic strike, dip and rake are taken from the published literature. The K27 and K29 models provide an advantage to such work here as they are developed from partitions of a rupture catalogue (Papazachos, 1999a,b; Papazachos et al., 1999) that has been extended (Weatherill and Burton, 2009). Each zone contains a set of ruptures from which ruptures in the synthetic catalogue can be sampled. In both cases uncertainty in the rupture parameters strike, dip and rake is modelled by random sampling from a Gaussian distribution. Characterisation of fault type is then done on the basis of the synthetic rake using the Frohlich and Apperson (1992) classification scheme. A hazard map for PGA with a 10% probability of being exceeded in 50 years, incorporating fault type and an assumed “rock” site, is given in Fig. 19. As expected, there is an increase in hazard along the Hellenic Arc and Adriatic belt, and also along the trace of the North Anatolian fault. Hazard is also reduced across much of continental Greece and western Turkey where normal faulting is prevalent. The change in hazard is typically on the order of 20 to 60 cm s−2.

273

8. Conclusions The Monte Carlo methodology and the hazard maps shown here illustrate the impact that alternative models and approaches have on seismic hazard assessment for a region. There are several interesting results that warrant some consideration in future seismic hazard analyses of the Aegean region. The acceptance of PP2000 as the definitive source model for the Aegean may be overlooking a fundamental area of epistemic uncertainty in hazard analyses. Whilst we recognise the difficulties in objective validation of source models, this analysis has shown that there is merit in considering other models or methods of source characterisation. Analyses of epistemic uncertainty should consider some of these alternatives, even if the analyses are not as exhaustive as those produced here. Furthermore, the uncertainties on a-value, b-value and MMAX are real. Taking the median value for each of these parameters without consideration of the uncertainties may ultimately result in an underestimation of the seismic hazard. The manner in which the Monte Carlo approach can incorporate the

Fig. 19. PGA (cm s−2) with a 10% probability of being exceeded in 50 years, for a “rock” site with variable faulting. (a) MCMAEU method is used with χ2 zone model weighting and EXP attenuation model weighting. (b) Change in hazard (compared to Fig. 12a) owing to the inclusion of variable faulting.

274

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

Fig. 20. PGA (cm s−2) with a 10% probability of being exceeded in 50 years, incorporating site condition and variable faulting. MCMAEU method is used with χ2 zone model weighting and EXP attenuation model weighting.

uncertainties would suggest that this approach should be considered alongside the standard Cornell (1968)–McGuire (1976) methodology. Analyses of the use of attenuation models for the Aegean should extend beyond those derived exclusively from Greek earthquakes. The analysis of the fit of PGA attenuation models here, whilst not exhaustive, suggests that European relations and the Boore and Atkinson (2007) NGA relation could also be applied, with a significant impact on the hazard assessment. There is justification for considering these models as Aegean earthquakes contribute to the data set used for their construction, strongly in the case of Bm07 and Am05. Whilst we have included Am05 in the epistemic uncertainty analysis, the high PGA for near field earthquakes means that smaller near-field earthquakes (sampled at 0.5σ to 2σ above the median) come to dominate the hazard event at the relatively short return period. For PGA application and short period accelerations the Bm07 model is preferred. Several extensions to the hazard assessment using the Monte Carlo method herein all reveal some interesting results. The dissemination of the global Vs30 database (Wald and Allen, 2007) has been a powerful development and allows site condition to be incorporated into seismic hazard maps. It also improves compatibility with real-time ground motion estimates made following large earthquakes using the USGS Shakemap program. Similarly, simple characterisation of fault type lays the ground work for improving the representation of the seismic source in seismic hazard maps. A substantial shortcoming of this analysis is the poor constraint of ground motion from intermediate depth earthquakes in the south Aegean. Attenuation modelling of these earthquakes will improve as more strong motion records become available. A final time-independent hazard map is constructed that combines all the elements considered previously (except intermediate depth seismicity) and is shown in Fig. 20. This demonstrates the PGA with a 10% probability of being exceeded in 50 years, using the MCMAEU method with the χ2 based zone weighting scheme and EXP attenuation weighting scheme (MMIN = Mw 5.2). Fault type and site condition are both incorporated into the map. The hazard map(s) that arise from this work can be used to form the basis of seismic zonation for the implementation of Eurocode 8 in the Aegean. They may also be used for vulnerability studies in urban areas across the region. Site-specific data should not be extracted from these maps. However, the Monte Carlo methodology adopted here has been shown to be capable of undertaking such assessment, albeit a more localised approach would be needed for a robust site-specific

analysis. This work is a progressive step toward improving the estimation of seismic hazard in the Aegean region, and identification of its uncertainties. Acknowledgements We are very grateful to Dario Slejko and one anonymous reviewer, whose detailed and insightful comments greatly improved this manuscript. This work has been funded by a University of East Anglia Faculty of Science Studentship. Appendix A The following relations describe the horizontal PGA attenuation models (or zero-period spectral attenuation models where defined as such) summarised in Table 2. A schematic illustrating the different source-site distance metrics is given in Fig. A.1. A.1. Ambraseys et al. (1996) [Am96] log ðPGAÞ = −1:48 + 0:266MS −0:922 log

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  R2JB + 3:52

ðA:1Þ

+ 0:117SA + 0:124SS + 0:25P

SA and SS are dummy variables assuming a value of 1 for stiff soil or soft soil respectively, and 0 otherwise. A.2. Skarlatoudis et al. (2003) [Sk03] log10 ðPGAÞ = 1:07 + 0:45MW −1:35 log10

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi REPI + 6

ðA:2Þ

+ 0:09F + 0:06S + 0:286P

F is a dummy variable describing the faulting, assuming a value of 1 for strike-slip or thrust faults, 0 for normal faults. S is a dummy variable describing the site condition, assuming a value of 0 for NEHRP Site Class B, 1 for Site Class C and 2 for Site Class D.

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

275

Fig. A.1. Schematic of four common source-site distance metrics for a) dipping fault, and b) a vertical fault. REPI = epicentral, RHYP = hypocentral, RJB = Joyner–Boore and RRUP = Rupture distance.

A.3. Ambraseys et al. (2005) [Am05]

A.6. Boore and Atkinson (2007) [BA07]

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R2JB + 7:62 log ðPGAÞ = 2:522−0:142MW + ð−3:184 + 0:314MW Þ log

    log ðY Þ = FM ðMW Þ + FD RJB ; MW + FS Vs30 ; RJB ; MW + εσT

+ ::::::0:134SS + 0:050SA −0:084FN + 0:062FT −0:044F0 + εP

ðA:3Þ

ðA:8aÞ

FD defines the distance scaling: !

SA and SS are defined as for Ambraseys et al. (1996). FN, FT and F0 assume a value of 1 for normal, thrust and oblique faulting respectively, 0 otherwise. ε is the total variance of the model defined by:

  h  i FD RJB ; MW = c1 + c2 Mw −Mref ln

  2 2 1 =2 ε = ð0:665−0:065MW Þ + ð0:222−0:022MW Þ

Where Rref, Mref assume fixed values of 1.0 km and Mw 4.5 respectively. c1, c2 and c3 are period-dependent coefficients. R is defined as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R = R2JB + h20 , where h0 = 1.35. FM defines the magnitude scaling:

ðA:4Þ



A.4. Danciu and Tselentis (2007) [DT07] log10 ðY Þ = 0:883 + 0:458MW −1:278 log10

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R2EPI + 11:5152

+ 0:038S + 0:116F + εP

ðA:5Þ

S and F are dummy variables defined as for Skarlatoudis et al. (2003). ε is the total variability separated into inter- (τ) and intra-event (σ) components, assume values of 0.109 and 0.270 respectively, ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiwhich ε = τ2 + σ 2 . A.5. Bommer et al. (2007) [Bm07] 2

log 10 ðY Þ = 0:0031 + 1:0848Mw −0:0835Mw qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi + ð−2:4423 + 0:2081MW Þ log10 R2JB + 8:02822

ðA:6Þ

SS, SA, FN and FR are defined as for Ambraseys et al. (2005). Aleatory variability is separated into inter- (τ) and intra-event (σ) components defined as:

2

e1 U + e2 SS + e3 NS + e4 RS + e5 ðMW −Mh Þ + e6 ðMW −Mh Þ e1 U + e2 SS + e3 NS + e4 RS + e7 ðMW −Mh Þ

if

M≤Mh M N Mh

ðA:8cÞ Where U, SS, NS and RS are dummy variables assuming values of 1 for unspecified, strike-slip, normal and reverse faulting respectively, 0 otherwise. Mh is a period-dependent “hinge coefficient”, which assumes a value of 6.75 for PGA. e1 to e5 are period-dependent coefficients. The site amplification function is determined from:     V FS VS30 ; RJB ; MW = bLIN ln S30 + FNL 760

ðA:8dÞ

Where Vs30 is the average 30-m shear wave velocity, bLIN is a period dependent coefficient (−0.360 for PGA). FNL is defined as:

ðA:8eÞ pga4nl refers to the median PGA on reference rock (Vs30 = 760 m s−1), c and d are defined by:

ðA:7aÞ c = ð3Δy−bNL ΔxÞ = Δx

σ = 0:599−0:058MW

  + c3 R−Rref ðA:8bÞ

8 bNL lnð0:6Þ > pga4nl≤0:03 > > > > h  . i2 h  . i3 > pga4nl pga4nl < bNL lnð0:006Þ + c ln + d ln for 0:03≤pga4nl≤0:09 FNL = 0:03 0:03 >  .  > > pga4nl > > bNL ln 0:09bpga4nl > : 0:1

+ ::::::0:0781SS + 0:0208SA −0:0292FN + 0:0963FR

τ = 0:323−0:031MW

FM ðMW Þ =

R Rref

ðA:7bÞ

2

d = −ð2Δy−bNL ΔxÞ = Δx

ðA:8fÞ 3

ðA:8gÞ

276

where Δx = ln according to:

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278



0:09

=0:03 Þ and Δy = bNL ln(0.9). bNL is determined

8 b1 > > > > < ðb1 −b2 Þ lnðVS30 = 300Þ = lnð180 = 300Þ + b2 bNL = for > > b lnðVS30 = 760Þ = lnð300 = 760Þ > > : 2 0:0

VS30 ≤180 180≤VS30 ≤300 300bVS30 b760

ðA:8hÞ

760≤VS30

where b1 and b2 are period-dependent coefficients (−0.64 and −0.14 respectively for PGA). The PGA inter-event variability is defined separately for earthquakes with a specified fault mechanism (τM = 0.260), and for those with an unspecified fault mechanism (τU = 0.265). Intra-event variability for PGA is fixed at σ = 0.502. Appendix B. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.tecto.2010.06.022. In addition to the references cited in this paper, the following are cited in the electronic supplementary material: Abercrombie et al. (1995), Anderson and Jackson (1987), Armijo et al. (1992, 2002), Avallone et al. (2004), Baker et al. (1997), Benetatos et al. (2004), Bernard et al. (2006), Caputo (1996), Chatzipetros et al. (2005), Clarke et al. (1997), Erdik et al. (1999), Karacostas et al. (1999), Louvari et al. (1999), Moretti et al. (2003), Papadimitriou et al. (2006) and Tranos et al. (2003). References Abercrombie, R.E., Main, I.G., Douglas, A., Burton, P.W., 1995. The nucleation and rupture process of the 1981 Gulf of Corinth earthquakes from deconvolved broadband data. Geophysical Journal International 120, 393–405. Abrahamson, N., 2000. State of the practice of seismic hazard evaluation. GEOEng 2000. Australian Geomechanics Society, Melbourne, Australia, pp. 1–27. Abrahamson, N.A., Bommer, J.J., 2005. Probability and uncertainty in seismic hazard analysis. Earthquake Spectra 21, 603–607. Allen, T.I., Wald, D.J., 2009. On the use of high resolution topographic data as a proxy for seismic site conditions (VS30). Bulletin. Seismological Society of America 99, 935–943. Ambraseys, N.N., 1995. The prediction of peak ground acceleration in Europe. Earthquake Engineering and Structural Dynamics 24, 467–490. Ambraseys, N.N., Bommer, J.J., 1991. Database of European strong ground-motion records. European Earthquake Engineering 2, 18–37. Ambraseys, N.N., Simpson, K.A., Bommer, J.J., 1996. Prediction of horizontal response spectra in Europe. Earthquake Engineering and Structural Dynamics 25, 371–400. Ambraseys, N.N., Smit, P., Douglas, J., Margarsi, B., Sigbjornsson, R., Olafsson, S., Suhadolc, P., Costa, G., 2004. Internet-site for European strong motion data. Bollettino di Geofisica Teorica ed Applicata 45, 113–129. Ambraseys, N.N., Douglas, J., Sarma, S.K., Smit, P.M., 2005. Equations for the estimation of strong ground motions from shallow crustal earthquakes using data from Europe and the Middle East: horizontal peak ground acceleration and spectral acceleration. Bulletin of Earthquake Engineering 3, 1–53. Anderson, H., Jackson, J., 1987. Active tectonics of the Adriatic Region. Geophysical Journal of the Royal Astronomical Society 91, 937–983. Anooshehpoor, A., Brune, J.N., Zheng, Y., 2004. Methodology for obtaining constraints on ground motion from precariously balanced rocks. Bulletin Seismological Society of America 94, 285–303. Armijo, R., Lyon-Caen, H., Papanastassiou, D., 1992. East-west extension and Holocene normal fault scarps in the Hellenic arc. Geology 20, 491–494. Armijo, R., Meyer, B., Navarro, S., King, G., Barka, A., 2002. Asymmetric slip partitioning in the Sea of Marmara pull-apart: a clue to propagation processes of the North Anatolian Fault. Terra Nova 14, 80–86. Atkinson, G.M., Boore, D.M., 2003. Empirical ground-motion relations for subduction-zone earthquakes and their application to Cascadia and other regions. Bulletin Seismological Society of America 93, 1703–1729. Avallone, A., Briole, P., Agatza-Balodimou, A.M., Billiris, H., Charade, O., Mitsakaki, C., Nercessian, A., Papazissi, K., Paradissis, D., Veis, G., 2004. Analysis of eleven years of deformation measure by GPS in the Corinth Rift Laboratory area. Comptes Rendus Geosciences 336, 301–311. Baker, C., Hatzfield, D., Lyon-Caen, H., Papadimitriou, E., Rigo, A., 1997. Earthquake mechanisms of the Adriatic Sea and Western Greece: implications for the oceanic subduction-continental collision transition. Geophysical Journal International 131, 559–594. Beauval, C., Scotti, O., 2004. Quantifying sensitivities of PSHA for France to earthquake catalog uncertainties, truncation of ground-motion variability and magnitude limits. Bulletin Seismological Society of America 94, 1579–1594.

Beauval, C., Bard, P.-Y., Hainzl, S., Gueguen, P., 2008. Can strong-motion observations be used to constrain probabilistic seismic-hazard estimates? Bulletin Seismological Society of America 98, 509–520. Bender, B., Campbell, K.W., 1989. A note on the selection of minimum magnitude for use in seismic hazard analysis. Bulletin Seismological Society of America 79, 199–204. Benetatos, C., Kiratzi, A., Papazachos, C., Karakaisis, G., 2004. Focal mechanisms of shallow and intermediate depth earthquakes along the Hellenic Arc. Journal of Geodynamics 37, 253–296. Bernard, P., Lyon-Caen, H., Briole, P., Deschamps, A., Boudin, F., Makropoulos, K., Papadimitriou, P., Lemeille, F., Patau, G., Billiris, H., Paradissis, D., Papazissi, K., Castarede, H., Charade, O., Nercessian, A., Avallone, A., Pacchiani, F., Zahradnik, J., Sacks, S., Linde, A., 2006. Seismicity, deformation and seismic hazard in the western rift of Corinth: new insights from the Corinth Rift Laboratory (CRL). Tectonophysics 426, 7–30. Beyer, K., Bommer, J.J., 2006. Relationships between median values and between aleatory variabilities for different definitions of the horizontal component of motion. Bulletin Seismological Society of America 96, 1512–1522. Bommer, J.J., Scherbaum, F., 2008. The use and misuse of logic trees in probabilistic seismic hazard analysis. Earthquake Spectra 24, 997–1009. Bommer, J.J., Douglas, J., Strasser, F.O., 2003. Style-of-faulting in ground-motion prediction equations. Bulletin of Earthquake Engineering 1, 171–203. Bommer, J.J., Abrahamson, N.A., Strasser, F.O., Pecker, A., Bard, P.Y., Bungum, H., Cotton, F., Fäh, D., Sabetta, F., Scherbaum, F., Studer, J., 2004. The challenge of defining upper bounds on earthquake ground motions. Seismological Research Letters 75 (1), 82–94. Bommer, J.J., Stafford, P.J., Alarcon, J.E., Akkar, S., 2007. The influence of magnitude range on empirical ground motion. Bulletin Seismological Society of America 97, 2152–2170. Boore, D.M., Atkinson, G.A., 2007. Boore and Atkinson NGA ground-motion relations for the geometric mean horizontal component of peak and spectral ground motion parameters. PEER Report. University of California, Berkeley. Boore, D.M., Joyner, W.B., Fumal, T.E., 1997. Equations for estimating horizontal response spectra and peak acceleration from Western North American earthquakes: a summary of recent work. Seismological Research Letters 68, 128–153. Boore, D.M., Watson-Lamprey, J., Abrahamson, N.A., 2006. Orientation-independent measures of ground motion. Bulletin Seismological Society of America 96, 1502–1511. Brune, J.N., 1996. Precariously balanced rocks and ground-motion maps for southern California. Bulletin Seismological Society of America 86, 43–54. Burton, P.W., Xu, Y., Tselentis, G.A., Sokos, E., Aspinall, W., 2003. Strong ground acceleration seismic hazard in Greece and neighbouring regions. Soil Dynamics and Earthquake Engineering 23, 159–181. Burton, P.W., Xu, Y., Qin, C., Tselentis, G.A., Sokos, E., 2004a. A catalogue of seismicity in Greece and the adjacent areas for the twentieth century. Tectonophysics 390, 117–127. Burton, P.W., Xu, Y., Qin, C., Tselentis, G.A., Sokos, E., 2004b. Extreme earthquake and earthquake perceptibility study in Greece and its surrounding area. Natural Hazards 32, 277–312. Campbell, K.W., 1997. Empirical near-source attenuation relationships for horizontal and vertical components of peak ground acceleration, peak ground velocity, and pseudo-absolute acceleration response spectra. Seismological Research Letters 68, 154–179. Campbell, K.W., Bozorgnia, Y., 2006. Next Generation Attenuation (NGA) empirical ground motion models: can they be used in Europe? Proceedings of the First European Conference on Earthquake Engineering and Seismology, Paper No. 458, Geneva, Switzerland. Caputo, R., 1996. The active Nea Anchialos Fault System (Central Greece): comparison of geological, morphotectonic, archealogical and seismological data. Annali di Geofisica 34, 557–574. Chatzipetros, A., Kokkalas, S., Pavlides, S., Koukouvelas, I., 2005. Paleoseismic data and their implication for active deformation in Greece. Journal of Geodynamics 40, 170–188. Chiou, B., Darragh, R., Gregor, N., Silva, W., 2008. NGA Project Strong-Motion Database. Earthquake Spectra 24, 23–44. Clarke, P.J., Davies, R.R., England, P.C., Parsons, B.E., Billiris, H., Paradissis, D., Veis, G., Denys, P.H., Cross, P.A., Ashkenazi, V., Bingley, R., 1997. Geodetic estimate of seismic hazard in the Gulf of Korinthos. Geophysical Research Letters 24, 1303–1306. Coley, D.A., 2005. An introduction to genetic algorithms for scientists and engineers. World Scientific, Singapore. edn, Vol., pp. Pages,. Coppersmith, K.J., Youngs, R.R., 1986. Capturing uncertainty in probabilistic seismic hazard assessments within intraplate tectonic environments. Proceedings of the Third US National Conference on Earthquake Engineering, Charleston, South Carolina, pp. 24–28. August. Cornell, C.A., 1968. Engineering Seismic Risk Analysis. Bulletin Seismological Society of America 58, 1583–1606. Cotton, F., Scherbaum, F., Bommer, J.J., Bungam, H., 2006. Criteria for selecting and adjusting ground -motion models for specific target regions: application to Central Europe and rock sites. Journal of Seismology 10, 137–156. Cramer, C.H., Petersen, M.D., Reichle, M.S., 1996. A Monte Carlo approach in estimating uncertainty for a seismic hazard assessment of Los Angeles, Ventura, and Orange Counties, California. Bulletin Seismological Society of America 86, 1681–1691. Danciu, L., Tselentis, G.-A., 2007. Engineering ground-motion parameters attenuation relationships for Greece. Bulletin Seismological Society of America 97, 162–183. Ebel, J.E., Kafka, A.L., 1999. A Monte Carlo approach to seismic hazard analysis. Bulletin Seismological Society of America 89, 854–866. Erdik, M., Biro, Y.A., Onur, T., Sesetyan, K., Birgoren, G., 1999. Assessment of earthquake hazard in Turkey and neighbouring regions. Annali di Geofisica 42, 1125–1138.

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278 Frankel, A., 1995. Mapping Seismic Hazard in the Central and Eastern United States. Seismological Research Letters 66, 8–21. Frohlich, C., Apperson, K.D., 1992. Earthquake focal mechanisms, moment tensors and consistency of seismic activity near plate boundaries. Tectonics 11, 279–296. Grünthal, G., Wahlström, R., 2001. Sensitivity of parameters for probabilistic seismic hazard analysis using a logic tree approach. Journal of Earthquake Engineering 5, 309–328. Grünthal, G., Bosse, C., Sellami, S., Mayer-Rosa, D., Giardini, D., 1999. Compilation of the GSHAP regional seismic hazard map for Europe, Africa and the Middle East. Annali di Geofisica 42, 1215–1223. Gutenberg, B., Richter, C.F., 1944. Frequency of Earthquakes in California. Bulletin Seismological Society of America 34, 1985–1988. Hanks, T.C., Kanamori, H., 1979. A moment-magnitude scale. Journal of Geophysical Research 84, 2348–2350. Hatzidimitriou, P.M., Papadimitriou, E.E., Mountrakis, D.M., Papazachos, B.C., 1985. The seismic parameter b of the frequency-magnitude relation and its association with the geological zones in the area of Greece. Tectonophysics 120, 141–151. Johnson, C.E., Koyanagi, R.Y., 1988. A Monte-Carlo approach applied to the estimation of seismic hazard for the state of Hawaii. Seismological Research Letters 59, 18. Karacostas, V.G., Papadimitriou, E.E., Papazachos, C.B., 2004. Properties of the 2003 Lefkada, Ionian Islands, Greece, Earthquake Seismic Sequence and Seismicity Triggering. Bulletin Seismological Society of America 94, 1976–1981. Kijko, A., 2004. Estimation of the maximum earthquake magnitude, Mmax. Pure and Applied Geophysics 161, 1655–1681. Kijko, A., Lasocki, S., Graham, G., 2001. Non-parametric Seismic Hazard in Mines. Pure and Applied Geophysics 158, 1655–1675. Kiratzi, A., Louvari, E., 2003. Focal mechanisms of shallow earthquakes in the Aegean Sea and the surrounding lands determined by waveform modelling: a new database. Journal of Geodynamics 36, 251–274. Konstantinou, K.I., Kalogeras, I.S., Melis, N.S., Kourouzidis, M.C., Stavrakakis, G.N., 2006. The 8 January 2006 earthquake (Mw 6.7) offshore Kythira Island, Southern Greece: seismological, strong-motion, and macroseismic observations of an intermediatedepth event. Seismological Research Letters 7, 544–553. Kramer, S.L., 1996. Geotechnical Earthquake Engineering, 1st edition. Prentice Hall, New Jersey. Kramer, S.L., Mitchell, R.A., 2006. Ground motion intensity measures for liquefaction hazard evaluation. Earthquake Spectra 22, 413–438. Louvari, E., Kiratzi, A.A., Papazachos, B.C., 1999. The Cephalonia Transform Fault and its extension to western Lefkada Island (Greece). Tectonophysics 308, 223–236. Makropoulos, K.C., Burton, P.W., 1983. Seismic risk of circum-pacific earthquakes i: strain-energy release. Pure and Applied Geophysics 121, 247–267. Makropoulos, K.C., Burton, P.W., 1985a. Seismic hazard in Greece. I: magnitude recurrence. Tectonophysics 117, 205–257. Makropoulos, K.C., Burton, P.W., 1985b. Seismic Hazard in Greece II: ground Acceleration. Tectonophysics 117, 259–294. Mantyniemi, P., Tsapanos, T.M., Kijko, A., 2004. An estimate of probabilistic seismic hazard for five cities in Greece by using the parametric-historic procedure. Engineering Geology 72, 217–231. Margaris, B., Papazachos, C., Papaioannou, Ch., Theodulidis, N., Kalogeras, I., Skarlatoudis, A., 2001. Ground motion attenuation relations for shallow earthquakes in Greece. 12th European Conference on Earthquake Engineering. McClusky, S., Balassanian, S., Barka, A., Demir, C., Ergintav, S., Georgiev, I., Gurkan, O., Hamburger, M., Hurst, K., Kahle, H., Kastens, K., Kekelidze, G., King, R., Kotzev, V., Lenk, O., Mahmoud, S., Mishin, A., Nadariya, M., Ouzounis, A., Paradissis, D., Peter, Y., Prilepin, M., Reilinger, R., Sanli, I., Seeger, H., Tealeb, A., Toksöz, M.N., Veis, G., 2000. Global Positioning System constraints on plate kinematics and dynamics in the eastern Mediterranean and Caucasus. Journal of Geophysical Research 105, 5695–5720. McGuire, R.K., 1976. FORTRAN computer program for seismic risk analysis, in Open File Report 76-67, ed USGS. USGS, Menlo Park. Montaldo, V., Faccioli, E., Zonno, G., Akinci, A., Malagnini, L., 2005. Treatment of groundmotion predictive relationships for the reference seismic hazard map of Italy. Journal of Seismology 9, 295–316. Moretti, I., Sakellariou, D., Lykousis, V., Micarelli, L., 2003. The Gulf of Corinth: an active half graben? Journal of Geodynamics 36, 323–340. Muir-Wood, R., 1993. From global seismotectonics to global seismic hazard. Annali di Geofisica 36 (3–4), 153–168. Musson, R.M.W., 1999a. Probabilistic seismic hazard maps for the North Balkan region. Annali di Geofisica 42, 1109–1124. Musson, R.M.W., 1999b. Determination of design earthquakes in seismic hazard analysis through Monte Carlo simulation. Journal of Earthquake Engineering 3, 463–474. Musson, R.M.W., 2000. The use of Monte Carlo simulations for seismic hazard assessment in the U.K. Annali di Geofisica 43, 1–9. Musson, R.M.W., 2004. Objective validation of source models. 13th World Conference on Earthquake Engineering, Paper No. 2492, Vancouver, Canada. Musson, R.M.W., 2005. Against Fractiles. Earthquake Spectra 21, 887–891. Musson, R.M.W., 2009. Ground motion and probabilistic hazard. Bulletin of Earthquake Engineering 7, 575–589. Musson, R.M.W., Henni, P.H.O., 2001. Methodological considerations of probabilistic seismic hazard mapping. Soil Dynamics and Earthquake Engineering 21, 385–403. Musson, R.M.W., Sargeant, S.L., 2007. Eurocode 8 seismic hazard zoning maps for the U.K. Technical Report, CR/07/125. British Geological Survey, pp. 1–62. Nyst, M., Thatcher, W., 2004. New constraints on the active tectonic deformation of the Aegean. Journal of Geophysical Research 109. doi:10.1029/2003JB002830. O'Hara, T.F., Jacobson, J.P., 1991. Standardization of the Cumulative Absolute Velocity. Report No. EPRI-TR-10082. Electric Power Research Institute. Palo Alto, CA.

277

Ordaz, M., Reyes, C., 1999. Earthquake hazard in Mexico City: observations versus computations. Bulletin Seismological Society of America 89, 1379–1383. Papadimitriou, P., Kaviris, G., Makropoulos, K., 2006. The Mw = 6.3 2003 Lefkada earthquake (Greece) and induced stress transfer changes. Tectonophysics 423, 73–82. Papaioannou, C.A., Papazachos, B.C., 2000. Time-independent and time-dependent seismic hazard in Greece based on seismogenic sources. Bulletin Seismological Society of America 90, 22–33. Papazachos, B.C., 1990. Seismicity of the Aegean and surrounding area. Tectonophysics 178, 287–308. Papazachos, C.B., 1999a. Seismological and GPS evidence for the Aegean–Anatolia interaction. Geophysical Research Letters 26, 2653–2656. Papazachos, C.B., 1999b. An alternative method for a reliable estimation of seismicity with an application in Greece and the surrounding area. Bulletin Seismological Society of America 89, 111–119. Papazachos, B.C., Papazachou, C., 1997. The earthquakes of Greece, 1st edition. P. Ziti & co., Thessaloniki, Greece. Papazachos, B.C., Papaioannou, C.A., Papazachos, C.B., Savvidis, A.S., 1999. Rupture zones in the Aegean region. Tectonophysics 308, 205–221. Reasenberg, P., 1985. Second-order moment of Central California seismicity, 1969–1982. Journal of Geophysical Research 90, 5479–5495. Reilinger, R., McClusky, S., Vernant, P., Lawrence, S., Ergintav, S., Cakmak, R., Ozener, H., Kadirov, F., Guliev, I., Stepanyan, R., Nadariya, M., Hahubia, G., Mahmoud, S., Sakr, K., ArRajehi, A., Paradissis, D., Al-Aydrus, A., Prilepin, M., Guseva, T., Evren, E., Dmittrotsa, A., Filikov, S.V., Gomez, F., Al-Ghazzi, R., Karam, G., 2006. GPS constraints on continental deformation in the Africa-Arabia–Eurasia continental collision zone and implications for the dynamics of plate interactions. Journal of Geophysical Research 111. doi:10.1029/2005JB004051. Reiter, L., 1990. Earthquake hazard analysis: issues and insights, 1st edition. Columbia University Press, New York. Romeo, R., Paciello, A., Rinaldis, D., 2000. Seismic hazard maps of Italy including site effects. Soil Dynamics & Earthquake Engineering 20, 85–92. Scherbaum, F., Cotton, F., Smit, P., 2004. On the use of response spectral-reference data for the selection and ranking of ground-motion models for seismic hazard analysis in regions of moderate seismicity: the case of rock motion. Bulletin Seismological Society of America 94, 2164–2185. Shapira, A., 1983. Potential earthquake risk estimations by application of a simulation process. Tectonophysics 95, 75–89. Sinadinovski, C., Edwards, M., Corby, N., Milne, M., Dale, K., Dhu, T., Jones, A., McPherson, A., Gray, D., Robinson, D., White, J., 2005. Earthquake Risk. In: Jones, T., Middleman, M., Corby, N. (Eds.), Natural Hazard and Risk in Perth, WA. Geoscience Australia, Canberra, Australia. Skarlatoudis, A.A., Papazachos, C.B., Margaris, B.N., Theodolidis, N., Papaioannou, Ch., Kalogeras, I., Scordilis, E.M., Karakostas, V., 2003. Empirical peak ground motion predictive relationships for shallow earthquakes in Greece. Bulletin Seismological Society of America 93, 2591–2603. Slejko, D., Rebez, A., 2002. Probabilistic seismic hazard assessment and deterministic ground shaking scenarios for Vittorio Veneto (N.E. Italy). Bolletino Geofisica Teorica e Applicata 43, 263–280. Smith, W.D., 2003. Earthquake hazard and risk assessment in New Zealand by Monte Carlo methods. Seismological Research Letters 74, 298–304. Stafford, P.J., Strasser, F.O., Bommer, J.J., 2008. An evaluation of the applicability of the NGA models to ground-motion prediction in the Euro-Mediterranean region. Bulletin of Earthquake Engineering 6, 149–177. Stepp, J.C., 1971. An investigation of earthquake risk in the Puget Sound area by use of the type I distribution of largest extremes., Ph.D. Thesis, Pennsylvania State Univeristy, Pennsylvania. Stirling, M.W., Anooshehpoor, A., 2006. Constraints on probabilistic seismic hazard models from unstable landform features in New Zealand. Bulletin Seismological Society of America 96, 404–414. Strasser, F.O., Bommer, J.J., Abrahamson, N.A., 2008. Truncation of the distribution of ground motion residuals. Journal of Seismology 12, 79–105. Theodulidis, N.P., Papazachos, B.C., 1992. Dependence of strong ground motion on magnitude-distance, site geology and macroseismic intensity for shallow earthquakes in Greece: I, Peak horizontal acceleration, velocity and displacment. Soil Dynamics and Earthquake Engineering 11, 387–402. Tranos, M.D., Papadimitrou, E.E., Kilias, A.A., 2003. Thessaloniki-Gerakarou Fault Zone (TGFZ): the western extension of the 1978 Thessaloniki earthquake fault (Northern Greece) and seismic hazard assessment. Journal of Structural Geology 25, 2109–2123. Tselentis, G.-A., Danciu, L., 2008. Empirical relationships between modified Mercalli intensity and engineering ground motion parameters in Greece. Bulletin Seismological Society of America 98 (4), 1863–1875. Tselentis, G.-A., Danciu, L., 2010a. Probabilistic seismic hazard assessment in Greece — Part I: Engineering ground motion parameters. Natural Hazards and Earth System Sciences 10, 25–39. Tselentis, G.–.A., Danciu, L., 2010b. Probabilistic seismic hazard assessment in Greece– Part III: deaggregation. Natural Hazards and Earth System Sciences 10, 51–59. Tselentis, G.-A., Danciu, L., Sokos, E., 2010. Probabilistic seismic hazard assessment in Greece Part II: acceleration response spectra and elastic input energy spectra. Natural Hazards and Earth System Sciences 10, 41–49. Utsu, T., 1999. Representation and analysis of the earthquake size distribution: a historical review and some new approaches. Pure and Applied Geophysics 155, 509–535. Vannucci, G., Pondrelli, S., Argnani, A., Morelli, A., Gasperini, P., Boschi, E., 2004. An atlas of Mediterranean seismicity. Annals of Geophysics 47, 247–306 Supplement to.

278

G. Weatherill, P.W. Burton / Tectonophysics 492 (2010) 253–278

Wald, D.J., Allen, T.I., 2007. Topographic slope as a proxy for seismic site conditions and amplification. Bulletin Seismological Society of America 97, 1379–1395. Weatherill, G.A., 2009. A Monte Carlo Approach to Probabilistic Seismic Hazard Analysis in the Aegean Region. PhD Thesis. University of East Anglia, Norwich, United Kingdom. Weatherill, G.A., Burton, P.W., 2009. Delineation of shallow seismic sources using Kmeans cluster analysis, with specific application to the Aegean region. Geophysical Journal International 176, 565–588. Wiemer, S., 2001. A software package to analyze seismicity: ZMAP. Seismological Research Letters 72, 374–383.

Wiemer, S., Wyss, M., 2000. Minimum magnitude of completeness in earthquake catalogs: examples from Alaska, the Western United States, and Japan. Bulletin Seismological Society of America 90, 859–869. Wiemer, S., Giardini, D., Fah, D., Deichmann, N., Sellami, S., 2009. Probabilistic seismic hazard assessment of Switzerland: best estimates and uncertainties. Journal of Seismology 13 (4), 449–478. Woessner, J., Wiemer, S., 2005. Assessing the quality of earthquake catalogues: estimating the magnitude of completeness and its uncertainty. Bulletin Seismological Society of America 95, 684–698. Youngs, R.R., Chiou, S.-J., Silva, W.J., Humphrey, J.R., 1997. Strong ground motion attenuation relationships for subduction zone earthquakes. Seismological Research Letters 68, 58–73.