The 2013 European Seismic Hazard Model: key components and results

The 2013 European Seismic Hazard Model (ESHM13) results from a community-based probabilistic seismic hazard assessment supported by the EU-FP7 project “Seismic Hazard Harmonization in Europe” (SHARE, 2009–2013). The ESHM13 is a consistent seismic hazard model for Europe and Turkey which overcomes the limitation of national borders and includes a through quantification of the uncertainties. It is the first completed regional effort contributing to the “Global Earthquake Model” initiative. It might serve as a reference model for various applications, from earthquake preparedness to earthquake risk mitigation strategies, including the update of the European seismic regulations for building design (Eurocode 8), and thus it is useful for future safety assessment and improvement of private and public buildings. Although its results constitute a reference for Europe, they do not replace the existing national design regulations that are in place for seismic design and construction of buildings. The ESHM13 represents a significant improvement compared to previous efforts as it is based on (1) the compilation of updated and harmonised versions of the databases required for probabilistic seismic hazard assessment, (2) the adoption of standard procedures and robust methods, especially for expert elicitation and consensus building among hundreds of European experts, (3) the multi-disciplinary input from all branches of earthquake science and engineering, (4) the direct involvement of the CEN/TC250/SC8 committee in defining output specifications relevant for Eurocode 8 and (5) the accounting for epistemic uncertainties of model components and hazard results. Furthermore, enormous effort was devoted to transparently document and ensure open availability of all data, results and methods through the European Facility for Earthquake Hazard and Risk (www.efehr.org).

1 Introduction acceleration (PGA), exceeds a threshold level in a given period. A standard output of a PSHA is a map displaying the PGA level that has a 10 % exceedance probability in 50 years-an input parameter currently required by almost all National Annexes of Eurocode 8 (CEN 2004) as an anchor to Types I and II elastic design spectral shapes (e.g. Bisch et al. 2011). Despite the current importance of PGA in Eurocode 8, there are a number of other outputs of a PSHA that can and should be used in seismic design, including the uniform hazard spectra (UHS) for a range of return periods and the disaggregation of spectral ordinates in terms of the magnitude, distance and epsilon of the controlling earthquakes. Due to this direct connection to seismic design codes, PSHA has gained a central role in engineering applications, seismic risk analysis and the definition of the design basis ground-motions.
Modern PSHAs generate a wide range of results due to the growing understanding of the earthquake source processes, to the improved abilities for modelling ground motion attenuation, to an improved characterisation of uncertainties, and to enhanced computational capabilities (Pagani et al. 2014). Thus, the results of this comprehensive PSHA are hosted in a database accessible online via a web-portal. For the first time data for the entire European continent can be accessed interactively and used for applications in science, engineering, disaster mitigation, and in the insurance industry (Giardini et al. 2013a(Giardini et al. , 2014. PSHAs are performed to estimate hazard for a specific site and are extended to national or continental scales considering multiple sites. In Europe, the most recent efforts to provide a combined PSHA at the continental scale were completed over a decade ago (Giardini 1999;Jimenez et al. 2001), and since then most national hazard models have been updated at least once. National models or those covering several countries, however, are traditionally based on similar yet different procedures that are not harmonised and can result in considerable differences at country borders (e.g. Grünthal et al. 1998;Grünthal and Wahlström 2000;Wiemer et al. 2009b;Stucchi et al. 2011;Dominique and Andre 2012;Musson 2012). The project "Seismic Hazard Harmonization in Europe" (SHARE) was planned to overcome such differences and to benefit from the latest improvements in input data sets, from new approaches in modelling earthquake sources, and from a new generation of ground motion measurements and attenuation models. This project was conceived in 2007, submitted to the European Community in 2008 and funded in 2009 for a period of 3.5 years. Its main goal was to build a reference hazard model for Europe and Turkey, based on data compiled homogeneously without being constrained by national boundaries, and on full scientific integration of all the disciplines involved. In addition, SHARE aimed at providing the relevant input for the current version and planned update of the seismic provisions of the European building code (Eurocode 8).
SHARE brought together leading scientists from 18 European research institutions belonging to 12 countries. More than 50 researchers-seismologists, geologists, geodesists, historians, earthquake engineers, computer scientists, statisticians, and outreach specialists -formed the core team (see SHARE consortium members in the Acknowledgements). In addition, the contributions of more than 250 additional key European and worldwide experts were integrated during numerous meetings. These meetings were aimed at assembling high-quality datasets, constructing earthquake source models, discussing future ground motion expectations and their modelling, and defining procedural standards. An extensive, iterative review process was conducted with four major review meetings, aimed at discussing approaches and results with consortium members and key external experts (Fig. 1).
The first 2 years of the project were dedicated to compiling data, preparing prospective models and approaches, understanding the engineering requirements and preparing the computational workflow. The last 14 months focused on hazard calculations through a number of review meetings during which all aspects were thoroughly discussed, together with key experts from the wider engineering and seismological communities. Owing to the involvement of such a large number of European experts, these procedures are now setting a standard not only for Europe but also for future national efforts worldwide. The final seismic hazard maps are computed from Iceland, Great Britain and Portugal including the Canaries in the west to Finland, the Baltic countries, Poland, Romania and Turkey in the east.
SHARE delivered a comprehensive, community-based seismic hazard model that in the following we refer to as the 2013 European Seismic Hazard Model. The ESHM13 is the first regional effort to be completed since the end of the Global Seismic Hazard Assessment Program (GSHAP, Giardini 1999) and the follow-up project SESAME (Jimenez et al. 2001). At the same time, it is the first project to contribute its results to the Global Earthquake Model (GEM, www.globalquakemodel.org) initiative. Within the GEM framework, SHARE partners have closely collaborated with the leaders and scientists of the Earthquake Model for the Middle East (EMME, www.emme-gem.org) to discuss data and methods and to harmonise the hazard results across the national boundaries of south-eastern Turkey.
Any PSHA combines two essential ingredients: an earthquake source model that characterises the seismicity of the area of interest, and a ground motion model, that deals  Fig. 1 Workflow and timeline of the elicitation procedure within the SHARE project. Data, source model and GMPE Preparation workshops were held in the first period of the project, then joint meetings were held to discuss all issues of the hazard model. Key external experts were invited to these meetings (see Acknowledgements) with attenuation of the ground motion amplitudes through a set of ground motion prediction equations (GMPEs). Both are conceptually explained in the following. This document supplies a comprehensive scientific overview for the ESHM13. The reader may refer to further documentation on the SHARE/EFEHR websites (see ESHM13 data and resources) for further details such as the computational implementation of the model in OpenQuake (Pagani et al. 2014).

Data preparation and modelling approach
The ESHM13 is a time-independent hazard model, i.e. a model built under the assumption that earthquakes occur with a constant average frequency, but independently of each other (e.g. Cornell 1968;Reiter 1990). In other words, the model assumes that earthquakes have no memory of when and where other earthquakes occurred.
The area considered by SHARE extends from the Mid-Atlantic Ridge to the East European Platform, from the subduction zones in the Mediterranean and the large transform faults in Turkey to the Baltic shield. The study area thus encompasses highly diverse tectonic regimes. Consequently, the compilation of the aforementioned ingredients for seismic hazard assessment poses major challenges, both because of the seismotectonic diversity, reflected in the heterogeneous seismicity rates and seismic wave propagation properties, and because of the variable quantity and quality of the available data. Performing a regional hazard assessment requires compromises, and larger scale analyses inevitably carry the drawback of coarser characterisation of some areas. We first devised a concept of large areas, which we called superzones, to initially assess the spatial distribution of the following six model characteristics: prevalent tectonic regime, catalogue completeness, prior b-value, maximum magnitude, focal depth distribution and style of faulting. Knowing that these characteristics respond to different physical and tectonic constraints we developed a different set of superzones for each of them. Each superzone geometry is thus different. For example, the definition of a completeness-time history depends on the availability of historical documentation and instrumental recordings of any given region. Assessing maximum magnitudes and/or selecting ground motion prediction equations depends on the definition of the tectonic regime in which an event occurs, and hence on geologic and wave propagation properties as considered when defining the tectonic regimes (Delavaud et al. 2012). All these initial assessments were checked and revised if necessary for smaller scale features of the three final earthquake source models. For technical ease, the polygons enclosing the superzones were drawn to match the envelopes of the zonation models.
We developed three alternative earthquake source models, relying on different assumptions and providing a diverse description of the earthquake activity (details in Sect. 3): 1. An Area Source model (hereinafter AS-model) based on the definition of areal sources for which earthquake activity is defined individually. 2. A kernel-smoothed, zonation-free stochastic earthquake rate model that considers SEIsmicity and accumulated FAult moment (hereinafter SEIFA-model). Activity rates are based on the frequency-magnitude distribution model of the SHARE European Earthquake Catalogue (SHEEC, Grünthal and Wahlström 2012;Stucchi et al. 2012;Giardini et al. 2013a;Grünthal et al. 2013), while the spatial distribution of model rates depends on the density distributions of earthquakes and fault slip rates.

A Fault Source and
BackGround model (hereinafter FSBG-model), based on the identification of large seismogenic sources using tectonic and geophysical evidence. Activity rates are based on documented fault slip-rates.
We followed a hierarchical modelling approach. The set of superzone layers define large geographical regions providing an initial assessment essential for a homogeneous definition of model properties. These properties are then mapped onto the earthquake source models. The tectonic regionalisation is used for the appropriate regional selection of GMPEs.

Compiling homogeneous datasets
A milestone of the SHARE project was to compile harmonised datasets suitable for assessing the seismic hazard of the entire Euro-Mediterranean region. Harmonised datasets are the best way to handle the inherent heterogeneity in the tectonic environments and data availability as they help remove differences due to country-specific approaches and emphasize the differences due to actual seismotectonic diversity. Such datasets were compiled during the initial phases of the project (Fig. 1) and are described in the following sections. All data is available on-line through the European Facility for Earthquake Hazard and Risk website (http://www.efehr.org/).

The SHARE European Earthquake Catalogue
For the scope of the project we compiled a new homogeneous earthquake catalogue (SHEEC-the "SHARE European Earthquake Catalogue"; Stucchi et al. 2012;Grünthal et al. 2013 Data (AHEAD) (Locati et al. 2014) after reassessing all earthquake parameters, based either on raw macroseismic data or on existing, selected regional catalogues (Gomez-Capera et al. 2014);and(b) 1900-2006, harmonised from existing national catalogues and numerous special studies as described in Grünthal and Wahlström (2012) and Grünthal et al. (2013). The catalogue was extended to cover Central and Eastern Turkey for a more homogeneous hazard analysis. The assessment is documented as the SHARE earthquake catalogue for Central and Eastern Turkey (SHARE-CET, available at http://www.emidius. eu/SHEEC/docs/SHARE_CET.pdf). For the time-window 1000-2006 the catalogue lists over 30,000 earthquakes in the magnitude range 1.7 ≤ M W ≤ 8.5 (Fig. 2). The completeness of the SHEEC catalogue was assessed from a historical perspective within 22 completeness superzones. The superzones were defined as groups of single area sources within the AS-model, covering zones that exhibit a fairly homogeneous capacity of documenting the earthquake occurrence through history. The completeness time intervals defined at the scale of each superzone were reviewed for individual area sources to take into account local variations in seismicity and data availability. The results of the completeness assessment and the differences in the superzones are described in Appendix #3 of Stucchi et al. (2012). Superzones are groups of Area Sources (unrelated to country borders) also defined taking into account data homogeneity.
The catalogue was declustered with the windowing-approach using the window paramters of Burkhard and Grünthal (2009). Slightly more than 8000 of the original 30,000 events survived after declustering and completeness cut-off, implying specific challenges in some low-seismicity areas.

The European Database of Seismogenic Faults (EDSF)
The European Database of Seismogenic Faults (EDSF; Basili et al. 2013;Fig. 3), was specifically designed to fulfill the project requirements in terms of exploiting the wealth of Faults are colour-coded by slip rate (mm/y), thus illustrating the expected crustal deformation rates. Subduction zones are delineated with black lines. The gray background illustrates a preliminary model of strain rates in the Earth's crust inferred from geologic and geodetic data (see SHARE deliverable D3.5) Bull Earthquake Eng (2015) 13:3553-3596 3559 geologic, tectonic, and paleoseismic knowledge available in the scientific community. The database includes faults that are deemed capable of generating earthquakes of M W ≥ 5.5 and ensures homogenous input for use in probabilistic ground-shaking hazard assessment in the Euro-Mediterranean area. The chosen minimum magnitude threshold corresponds to a fault size of about 5 km and is assumed to define the resolution limit for identifying a seismogenic fault through geological investigations. The EDSF hosts two types of formalized seismogenic faults stored in two separate GIS layers: crustal faults and subduction zones.
For the parameterisation of crustal faults the EDSF adopts the Composite Seismogenic Source (CSS) model defined in previous works (Basili et al. 2008(Basili et al. , 2009 with standards used in other seismogenic fault databases (Haller and Basili 2011). The CSS is a simplified 3D model of a fault or fault system containing an unspecified number of aligned earthquake ruptures that could not-or simply were not-singled out. The fault is represented with a geo-referenced plane having variable strike but fixed dip along its entire length. Its surface projection is confined between two horizontal lines. In addition to this geographic feature, each fault is characterised by geometric parameters (minimum and maximum depth, strike, and dip) and behaviour parameters (rake, slip rate, and maximum earthquake magnitude).
For the parameterisation of subduction zones (SUBD) the EDSF adopts a newly defined model, partly following common definitions established in recent global efforts (Berryman et al. 2013). SUBD is a simplified 3D model of a subducting slab; similarly to the CSS model, the number and size of potential earthquake ruptures is unspecified. The slab is a complex surface represented by a georeferenced set of depth contours. In addition to its geographic attributes, each subduction zone is characterised by geometric parameters (minimum and maximum depth of the seismic interface and prevalent slab dip direction) and behaviour parameters (convergence azimuth, convergence rate, and maximum earthquake magnitude).
In order to capture the natural variability of the data (aleatory uncertainty) for both, crustal faults and subduction zones, all parameters are given as a min-max range. For each source the compilers also provide short explanations and qualified evidence on how each parameter was estimated. A list of pertinent references to scientific literature is also annexed to each seismogenic fault. For simplicity, we will in the following refer to seismogenic faults as fault sources, as opposed to area sources.
The EDSF was compiled through successive stages of update, revision, and homogenisation between June 2009 and May 2012, thanks to the contributions of over one hundred scientists from several European research institutions. It incorporated data from regional initiatives such as the well established DISS for Italy (DISS Working Group 2010) and surrounding areas, and the more recent GREDASS for Greece (Pavlides et al. 2010;Caputo et al. 2012Caputo et al. , 2013 and QAFI for Spain and Portugal (Garcia-Mayordomo et al. 2012). It also incorporated data from the sibling project EMME (Earthquake Model of the Middle East) for Turkey and from a number of local and regional studies. The final version of the database contains 1128 records for 63,775 km of crustal faults, from the Iberain peninsula to Turkey, and from three subduction zones known as Calabrian Arc, Hellenic Arc, and Cyprus Arc in the central and eastern Mediterranean Sea.

Modelling maximum magnitudes
For many hazard practitioners M max is best handled by using multiple values in a logic-tree with weights decreasing with increasing magnitude. Unfortunately it is impossible to test rigorously this parameter from a statistical point of view as this would require observation periods that extend far beyond human lifetime (Holschneider et al. 2014), thus M max remains a matter of choice, depending on the applied methodology.
For the ESHM13 M max is defined as the time-independent, ultimately largest earthquake that can occur in a region. The methods to assess M max are different for each model and explained in the following. In short, standard methods to assess the maximum magnitude were used for the AS-model (Melletti et al. 2009;Wheeler 2009). The SEIFA-model includes a single M max separately for crustal and subduction zone seismicity, and allows the largest event to occur anywhere geographically, however, with a model-driven cut-off leading to an actual spatial variation of M max . For the FSBG-model, M max is determined from the mean of multiple scaling relations that convert the rupture area to moment magnitude (Fig. 4). For the AS-model, the general approach was to give priority to actual data everywhere they are available and reliable, rather than to models. We adopted different approaches in low/moderate and in high seismicity regions. To this end we first defined four groups of area sources: (1) stable continental regions (SCR), (2) oceanic crust and mid-Atlantic ridge (OC), (3) Azores-Gibraltar zone (AG) and (4) generic active zones. For each source zone we defined four values of maximum magnitude which were then entered in a logic tree with decreasing weight values. The smallest and the largest values are displayed in Fig. 4a, b.
In low-to-moderate seismicity regions (mainly those in the stable continental regions) a single distribution was assumed, in analogy with the global analog approach (Wheeler 2009(Wheeler , 2011: the magnitude of the largest observed earthquake, with proper consideration of its uncertainty, was taken as the lower value for the distribution of maximum magnitude, i.e. M W = 6.5, whereas the other values were obtained by 0.2 increments (i.e. 6.7, 6.9, 7.1).
In the remaining seismogenic areas, normally characterised by high seismicity level and by a better understanding of historical seismicity and seismogenic sources, the distribution of M max was anchored to the larger value between the largest earthquake reported in the catalogue and the maximum magnitude expected based on the fault database, again with consideration of its uncertainty; this procedure was repeated in each AS in the active regions or in each superzone for oceanic crust, mid-Atlantic ridge and Azores-Gibraltar zones. Larger values of the M max distribution were obtained by 0.2 increments. The values obtained for the AS-model are shown in Fig. 5 together with their weights.
For the SEIFA-model, we defined M max = 8.6 for crustal seismicity as the truncation magnitude accounting for the largest moment magnitude earthquake listed in the SHEEC catalogue, i.e. the 1755 M W 8.5 Lisbon earthquake. For subduction zones we defined M max = 9.0 based on historically documented earthquakes (i.e. the 365AD, M8.5+ event in the Hellenic Arc; Stiros 2010), the overall size of the Hellenic Arc (e.g. Strasser et al. 2010), and subduction zone characteristics in the Mediterranean. However, it is important to note that magnitude bins with an annual rate smaller than 10 −6 are not considered in the hazard calculation. This leads to an apparent spatial variation of the maximum magnitude which is model-driven and not defined by a separate analysis, and thus constitutes a novel philosophy (Fig. 4c).
For the fault sources of the FSBG-model ( Fig. 4d) M max is calculated as the arithmetic mean of the estimates obtained from multiple fault scaling relations (Wells and Coppersmith 1994;Hanks and Bakun 2002;Kagan 2002;Leonard 2010). This additional assessment of M max , with respect to that given by compilers of the EDSF, was necessary to further homogenize the fault source dataset.
All faults in EDSF were classified according to their primary style of faulting and then the relevant fault scaling relations were applied. Considering that seismogenic faults in the EDSF may span more than one rupture in length, the maximum length of a single rupture was constrained based on aspect ratio for dip-slip faulting; conversely, no limitation was applied for strike-slip faulting. Given the fault classification and using the expected value for each fault scaling law (i.e. incorporating the epistemic uncertainty between relations), we obtained two to seven different moment magnitude estimates for each fault source. To capture also the internal variability of relations, being interested mainly in the positive side of the uncertainty, we augmented the mean M max with an uncertainty of ΔM W = 0.3 that represents a reasonable scatter of the relevant regressions. Colours depict the branching levels for the earthquake source models (yellow), maximum magnitude models (green) and ground motion models (red). Values below the black lines indicate the weights, for the source model these indicate the weighting scheme for different return periods as in Table 1. Tectonic regionalisation is not a branching level (grey) of the model, however, it defines the GMPEs to be used Bull Earthquake Eng (2015) 13:3553-3596 3563

Earthquake source models
A key step in our seismic hazard assessment is modeling epistemic uncertainty through a comprehensive logic-tree approach (Kulkarni et al. 1984;Budnitz et al. 1997;Bommer and Scherbaum 2008). The ESHM13 is based on three branching levels: (1) the earthquake source models, (2) the maximum magnitude, and (3) the ground motion models (Fig. 5); tectonic regionalisation is not a branching level but it is used to define the GMPEs set to be used. The earthquake source models (AS, FSBG, SEIFA) represent established approaches to model earthquake occurrences based on seismological, geological, tectonic and geodetic information, with a varying degree of importance represented in the source typologies (Fig. 6). As we discuss in better detail in Sect. 3.3, we assume full seismic coupling over the entire study region, i.e. that all tectonic and geodetic moment goes into seismic moment. We are aware that this is a simplification and may not be true everywhere (e.g. Ergintav et al. 2014); the coupling coefficient is indeed a crucial parameter, yet it remains rather controversial and very hard to assess homogenously at the scale of the entire continent. The three earthquake source models are characterised by alternative options to calculate and spatially distribute future seismic activity (Fig. 5). They were used to model crustal seismicity with depth ≤40 km. The seismicity of subduction zones and of the Vrancea region was modelled separately. Fig. 6 Schematic illustration of the four contributing source typologies used in the ESHM13: area sources with homogeneous distribution of rates; fault and background sources with M W ≥ 6.5 constrained to occur on fault sources and M W \ 6.5 events throughout the background zone; Kernel-smoothed model; and subduction zones model with seismicity on the interface modeled as complex fault, while in-slab seismicity is modeled as volumes at depth

Area source model
For several years the evaluation of seismic hazard has largely relied on area source models (Grünthal et al. 1998;Sleijko et al. 1998;Grünthal and GSHAP Working Group Region 3 1999;Musson 1999;Sleijko et al. 1999;Grünthal and Wahlström 2000;Musson 2000;Jimenez et al. 2001;Demircioglu et al. 2004;Meletti et al. 2008;Wiemer et al. 2009b;Papaioannou and Papazachos 2010), both in Europe and worldwide (e.g. Giardini 1999). Area source zones are designed with the assumption that seismicity may occur anywhere within each zone. The delineation considers seismicity, tectonics, geology and geodesy upon expert evaluation, generally without assigning uncertainties to the zone boundaries and only infrequently documenting the importance of the data used for the zonation (Burkhard and Grünthal 2009;Musson et al. 2009;Schmid and Slejko 2009;Wiemer et al. 2009a;Vilanova et al. 2014). Our zonation approach, based on the ESHM13 AS-model previously described, was applied to crustal seismicity down to a depth of 40 km.
The ESHM13 AS-model is based on the latest national area source models and on the former Euro-Mediterranean model (Jimenez et al. 2001) that have been merged and harmonised at national borders (Arvidsson and Grünthal 2010). Besides a number of political issues, challenges in the construction phase included the variability of seismotectonic environments, the differences in the nature of documentation, data availability and philosophies of the former models, and the need to harmonise the AS-model with the other datasets being compiled simultaneously within SHARE. In a final harmonisation step, small area sources that turned out to be devoid of seismicity were removed for consistency with the final datasets compiled within the project. The final area model consists of 432 area sources depicted as black lines in Fig. 7b.
The considerable differences in the tectonic environments across Europe are reflected in the final delineation of the area sources. In the Mediterranean area, and especially in the (a) compression-dominated areas including thrust or reverse faulting, associated transcurrent faulting (e.g., tear faults), and contractional structures in the upper plate of subduction zones (e.g., accretionary wedges), and (b) extension-dominated areas including associated transcurrent faulting, (c) major strike-slip faults and transforms, and (d) mid oceanic ridges, 4 subduction zones with contours of the depth of the dipping slab at 50 km intervals, 5 areas of deep-focus, nonsubduction earthquakes and 6 active volcanoes and other thermal/magmatic features. B Mapped tectonic regionalisation used for the selection of GMPEs. Area source zones are depicted as thin black lines and the tectonic regime is then used as a parameter for each area source zone Balkans, Greece, Turkey, and Italy, mapped active faults played a major role in defining the area sources. The knowledge on active faults of southern Spain, which was improved during the SHARE project, allowed for a redefinition of existing source models. In contrast, very few active structures are known in the continental shield of northern Europe. In these areas, seismicity was the main basis for the designing the area source model. Only a few major tectonic structures, such as the Rhine Graben in Central Europe, are represented in the zonation with a higher level of detail than large structures like the Baltic Shield. The AS-model was matched with the tectonic regionalisation (Fig. 7a, b) in order to define, for each area source, the prevailing style-of-faulting, the upper and lower bounds of the characteristic seismogenic depth, and the distribution of hypocentral depth.
As the declustered catalogue contains about 8000 events, a large number of the area sources contain only a few data; for instance, 196 zones contain 10 events or less. Those that did not contain any seismicity were assigned an activity rate by regional experts, or a water-level annual activity for M W ≥ 4.5 of 5 9 10 −8 events per km 2 .
We estimated activity rate parameters using a Bayesian penalized maximum likelihood (PML) approach (Johnston et al. 1994;Coppersmith et al. 2012) that updates the prior bvalue only if it is supported by local data. The prior b-values are estimated in the specific bvalue superzones with the method by Weichert (1980). The PML-estimates for the a-bvalue pairs and their associated uncertainty recover well the distribution of seismicity, also for sources with little data. The estimates are, however, guided by the more abundant small magnitude data. Two observations were made based on the PML-results as shown in the examples in Fig. 8-an apparent underestimation of rates in the moderate to large magnitude events in the range 5.5 ≤ M W ≤ 7.0, and a relatively large variation in b-values.
Like any maximum likelihood estimator, the PML-estimator is largely driven by the more abundant small magnitude seismicity. Thus, underestimating the rate of moderate to large events may result from how the catalogue completeness is assessed, from the variability of magnitude determinations for different time periods (e.g. historical vs. instrumental) and, more generally, from spatio-temporal variations of the observed seismicity. As such, this problem may be rooted both into the data and into the selected Uncertainties are displayed as confidence levels and single solution of the PML-approach. Black line indicates the frequency-magnitude distribution with a-/b-values according to expert judgement used in the hazard computation. The latter is the weighted cumulative sum of the frequency-magnitude distributions considering the four maximum magnitudes. The combination causes the stair-step like decay at the tail of the distribution methods of data analysis. In many regions the rates of small magnitude events are mainly governed by recent instrumental data, whereas the larger magnitude events are mainly historical earthquakes. Due to their origin from intensity assessments, the magnitudes of historical earthquakes can only assume a limited number of discrete values; as such they are binned differently from instrumental earthquakes, which may cause artificial perturbations in the seismicity rates.
We carefully analyzed each individual fit to the data and replaced the PML-results by subjective judgment in case b-values deviated more than Δb = 0.4 from the global estimate of the SHEEC catalogue (b = 0.9: Hiemer et al. 2014). The b-values obtained using the PML-approach for area sources characterised by very few data were constrained in the range 0.8-1.2, with the only exception of volcanic and swarm-type seismicity area sources, for which the b-value was allowed to increase up to 1.5 based on observed data.
The a-b-value pairs, with the a-value defined at M W = 0, were then used to calculate a weighted mean frequency-magnitude distribution (FMD) considering the weights of the Fig. 9 Parameter values of the crustal AS-model. a Annual rate of earthquakes for log 10 (λ(M W ≥ 4.5/km 2 )), b annual activity rate a-value: a(M W = 0)) and c b-value Bull Earthquake Eng (2015) 13:3553-3596 3567 maximum magnitude distribution (Fig. 8, solid black line). We calculated the incremental FMDs truncated at the different maximum magnitudes and then stacked them with the weights of the M max . Uncertainties are plotted for the 2.5 %/15 %/85 %/97.5 % percentiles of the distribution (Fig. 8, red and green dashed lines). We illustrate our results for two area sources, one in southwestern Switzerland (CHAS089) and one in the Central Apennines (ITAS308), containing 8 and 52 events above completeness, respectively (Fig. 8). For CHAS089, the difference between the two estimates is substantial (a = 2.4, b = 0.9; a pml = 2.71, b pml = 1.03; Fig. 8a). In contrast, for ITAS308 the difference between the estimates taken from the current Italian Seismic Hazard Model (a = 4.2, b = 1; Stucchi et al. 2011) and the PML (a pml = 3.92, b pml = 0.98) is negligible. Figure 8 also illustrates the difference in the uncertainties when estimating rates for magnitudes in the small and in the large magnitude ranges: for CHAS089 these uncertainties are much larger in the upper magnitude range than for ITAS308 due to the limited amount of data.
The AS-model activity rates are fully defined by the a-b-values and by the distribution of M max . The time-independent annual rate forecast for M W ≥ 4.5 varies across the continent and highlights the higher activity rates along the active fault systems of Southern Europe, Turkey and Iceland in comparison to the lower activity rates observed in stable continental regions (Fig. 9a). The rates are computed from the productivity parameter (avalue), that varies strongly across the continent, and from the b-value, that by construction varies only in a limited range (Fig. 9b, c).

Kernel-smoothed stochastic rate model considering seismicity and fault moment release (SEIFA-model)
In contrast to the other source models, this refined smoothed seismicity model uses estimates of the total productivity and frequency-magnitude distribution by fitting the entire declustered SHEEC catalogue using the catalogue completeness levels of the superzones simultaneously first. In a second step, earthquakes rates are spatially distributed according to a weighted combination of two spatial probability density functions (Hiemer et al. 2014). These are estimated from past seismicity and from accumulated moment release along faults, inferred from their slip rates and geometry. The relative weights of these two spatial density functions depend on the earthquake magnitude and are determined by pseudo-prospective forecasting experiments. The model, initially developed for California (Hiemer et al. 2013), is applied at the European scale: we generated a crustal model using earthquakes shallower than 40 km depth and active faults of the European Database of Seismogenic Faults (EDSF; Basili et al. 2013). We separately modelled seismicity in subduction zones in combination with the models of the Calabrian, Hellenic and Cyprus arcs with their complex 3D geometry. We treat active faults and subduction zones in the same way, as they are all characterised by their size, geometry and slip rate converted to moment release per unit area (discretized to 5 km 2 patches). The following sections provide a general description of the model; full details are given in Hiemer et al. (2014). The SEIFA-model assumes that (1) the frequency-magnitude distribution of past earthquakes is the best estimate of productivity and size distribution, (2) future earthquakes will occur in the vicinity of past events, (3) larger earthquakes occur more likely on mapped faults, and more likely along fast-slipping than along slow-slipping faults. These assumptions are parameterised to estimate the annual earthquake rate as a function of space and magnitude. The total productivity and the Gutenberg-Richter relations are estimated from the entire catalogue considering the space-time completeness histories as defined for the superzones (Hiemer et al. 2014). The cumulative annual earthquake rate forecast log 10 λ(M W ≥ 4.5) for crustal seismicity depicts a relatively smooth image with higher rates in currently active regions and along active fault systems (Fig. 10) as compared to the stable continental areas. This earthquake rate forecast shows more smaller scale variability in regions of frequent seismicity, a feature that can also be observed in the hazard map calculated only from this model (Fig. 16b).
The model allows for a maximum magnitude of M W = 8.6 for crustal events throughout Europe and Turkey. As the annual rates of larger magnitude earthquakes are locally very small, however, the actual maximum magnitude that entered the hazard calculations varies as we disregard magnitude bins for which the inferred annual rates are smaller than 10 −6 (Fig. 5c). According to our sensitivity tests, the contribution of these very low rate events to PGA hazard up to return periods of 5000 years is negligible, thus they were removed.

Fault sources and background (FSBG) model
In the FSBG-model, activity rates of the mapped fault sources are determined from the available deformation data as inferred from geodetic and geological methods (Basili et al. 2008(Basili et al. , 2009Haller and Basili 2011;Basili et al. 2013). The activity rates in this model are primarily dependent on the slip rate and the maximum magnitude of the fault sources. As crustal seismogenic faults are not known over the entire area covered by SHARE, and more specifically lack in regions of low tectonic deformation far from plate boundaries, the sources comprising the AS-model are used and adjusted in size and productivity at the edges of the background sources. Background sources were delineated with consideration for the distribution of fault sources. A background source polygon includes an entire fault system and does not cut across any fault source. A background zone is assumed to include a complete set of fault sources, i.e. one or several sources are embedded in one background source. Within one background zone we assume that all fault sources capable of generating M W ≥ 6.5 earthquakes are known and dominate the seismic moment release. Events with magnitudes M W ≥ 6.5 are modelled on the 3D geometry of the fault sources; their exact location varies depending on the modelled hypocentral location (Pagani et al. 2014), magnitude, and determined area, however, never extends beyond the mapped structures. Earthquakes in the magnitude range 4.5 ≤ M W ≤ 6.4 can occur throughout the background source and are homogeneously distributed.
The earthquake potential of each fault source is evaluated by converting the long-term tectonic moment rate into an earthquake activity rate. The tectonic moment rate is obtained from Ṁ 0 = µAḊc, where µ is rigidity, A is the area of the fault source, Ḋ is the slip rate and c is the seismic coupling, which indicates how much of the tectonic moment goes into the seismic moment. The earthquake rates are calculated assuming that all the accumulated tectonic moment is released seismically, thus the seismic coupling was set to c = 1 everywhere. We also adopted the maximum slip rate assigned to fault sources, thereby providing an upper bound of the cumulative tectonic moment.
After considering several methods for obtaining earthquake activity rates from slip rates (e.g. Anderson and Luco 1983;Youngs and Coppersmith 1985;Frankel et al. 2002;Stirling et al. 2002;Field et al. 2009;Stirling et al. 2012) we opted for the original contribution proposed by Anderson and Luco (1983), which generates activity rates that conform to a double-truncated Gutenberg-Richter distribution (Anderson and Luco 1983;Youngs and Coppersmith 1985;Bungum 2007).
The scale-invariant, power-law decay of the distribution is controlled by the b-value. Within a background zone the b-value is assumed to be the same everywhere, though it varies between background zones. Similarly to the AS-model, it is computed with the PML-method (Johnston et al. 1994;Coppersmith et al. 2012). The right-hand tail of the distribution, instead, is controlled by the maximum magnitude of the fault source. The total productivity within one background zone is the sum of the contribution of all fault sources. As a whole, the distribution complies with the moment conservation principle by summing up to the calculated seismic moment rate.
The model is illustrated for a background zone covering the Northern Apennines, Italy, that contains eight fault sources, all contributing to the overall productivity with their maximum slip rate (Fig. 11, source ITBG068, dashed grey lines). The b-value is estimated using the penalised maximum likelihood approach and equals b pml = 1.3. The summed activity using the maximum slip rates (Fig. 11, red line) defines well an upper bound of the cumulative frequency-magnitude distribution of the observed seismicity. Summing the productivity from the minimum slip rates matches the observed seismicity (Fig. 11, dashed green line) and overlies the PML-estimate. The tail of the distribution resembles a staircase function emerging through the summation of the activity rates calculated for each fault source that all have individual maximum magnitudes (Fig. 11, red line). We consistently used the maximum slip rates to estimate activity rates on fault sources which, as illustrated, often defined an upper bound of activity compared to the data observed in the catalogue. The FSBG-model covers mostly southern Europe and Turkey and some structures in Northern Europe such as the Lower Rhine Embayment (Fig. 12; Vanneste et al. 2013). Not all faults mapped and included in the EDSF could be used, however. For example, we excluded all mapped faults that represent only a small portion of a large background zone, thereby assuming that the fault coverage is not complete. The resulting earthquake rate forecast model for the cumulative annual rate [log 10 λ(M W ≥ 4.5)/km 2 ] is shown in Fig. 12 together with the tectonic association of the background zones and maximum slip rates of the fault sources.

Comparing the total rate forecasts for crustal seismicity
To obtain a global overview of the three crustal seismicity rate forecasts we calculated the total expected average annual earthquake rates and compared them with the cumulative frequency-magnitude distribution of the declustered earthquake catalogue using the superzone catalogue completeness (Fig. 13). For the AS-and the FSBG-models this involves summing up the rates of all individual zones. In contrast, by construction the SEIFA-model uses the combined maximum likelihood fit to the entire declustered catalogue (Hiemer et al. 2014).
Compared to the observed seismicity, Fig. 13 illustrates that the AS-model (blue) fits well the observed seismicity up to magnitude M W = 7, shows a slight overestimate up to M W = 8, and falls in between the rates of the largest magnitudes. The SEIFA-model Fig. 11 Activity rate plot for the source Northern Apennines (Source Id: ITBG068) that contains eight faults. Each fault contribution to the overall activity is indicated calculated with the specific maximum slip rate of the individual fault (dashed-grey lines). The total activity rate is the sum of the activity of the single faults (red line). Using the minimum slip rate on each fault adds to the dashed green distribution. Arrows indicate the ranges of magnitudes for earthquakes modelled as fault sources (M ≥ 6.5) or distributed in the entire background zone. Observed seismicity is display as cumulative and incremental frequency-magnitude distribution (triangles and circles, respectively) Bull Earthquake Eng (2015)  forecasts a smaller total earthquake rate compared to the observed seismicity in the magnitude range M W = 4.5-7.5, while it is qualitatively in good agreement with the observed rates in the uppermost magnitude range (dark yellow). In contrast, the FSBGmodel forecasts higher rates from M W = 4.5-7.2 and lower rates for the magnitude range above 7.5. The variability of the total rate forecasts illustrates the epistemic uncertainty in the source models and has several implications: 1. The frequency-magnitude distribution of the declustered earthquake catalogue (cut at the completeness levels) shows increased rates starting at M W = 6.5 as compared to Fig. 12 Fault sources and background model. a Annual rate of earthquakes for log 10 (λ(M W ≥ 4.5/km 2 )), and b fault sources, colour-coded by maximum slip rate (mm/y); background sources are colour-coded by tectonic regime, i.e. stable continental (blue) and active (brown) what is expected when extrapolating from smaller magnitudes. Reasons for this increased rate may arise from the completeness intervals, the declustering procedure and the different approaches used to determine magnitudes-or a combined effect. 2. The activity rates calculated from slip rates on fault sources (FSBG-model) show an underestimation of M W ≥ 7.5 events, while events of M W ≤ 7.0 are overestimated; the range in between is matched qualitatively well. We assume a Gutenberg-Richter type behaviour of faults constraining the maximum magnitude based on a weighted combination of scaling relations-the weighted M max are smaller compared to those used in the area source model. From an overall perspective this leads to an overestimation of small events, however, with considerable heterogeneity between slow-slipping and fast-slipping faults. The overestimation of smaller events could be reduced by considering a larger M max , which is however not supported by the current data. An alternative could be to consider characteristic type earthquake rate distributions such as in Schwartz and Coppersmith (1984). While it remains to be tested whether this would lead to a better representation of observed earthquake rates (e.g. Kagan et al. 2012), a combination of a Gutenberg-Richter type and a characteristic model might comprise a viable alternative. 3. The SEIFA-model forecasts a slightly lower activity rate compared to the observed data ( Fig. 13, dark yellow line). The power-law fit considers the space-time history for catalogue completeness from all completeness superzones (Stucchi et al. 2012;Hiemer et al. 2014). The visual comparison in the upper magnitude range is qualitatively satisfying. The forecasts are consistent with the SHEEC catalogue as well as with data from the Global CMT catalogue as shown in quantitative retrospective and pseudoprospective tests (Hiemer et al. 2014). 4. The mean cumulative rate forecast of the AS-model (Fig. 13, blue line) obtained from individual area sources shows a good agreement with the observed data for the magnitude range 4.5 ≤ M W ≤ 6.5 and slightly higher rates in the magnitude range 6.7 ≤ M W ≤ 7.7. Qualititatively the model matches well the uppermost range of magnitudes.
Combining the three source models within the logic-tree according to the assigned weights (0.5 for the AS-model, 0.3 for the SEIFA-model and 0.2 for the FSBG-model: Fig. 6; Table 1) we obtain a good qualitative match with the observed data (Fig. 13, red line). The three source models therefore generate forecasts that describe well the epistemic uncertainty of earthquake activity for crustal seismicity in Europe. Spatial variability in matching observed seismicity with the different models obviously exists across the study region. However, these are not established well enough to modify weights on a regional basis.

Modelling subduction zones and deep seismicity in Vrancea
We separately modelled the deep seismicity of the Calabrian, Hellenic, and Cyprus Arcs subduction zones as well as the non-subduction deep seismicity of the Vrancea Region (Romania). The activity rates of the three subduction zones has been jointly computed for interface and in-slab seismicity as depth uncertainties did not allow for a separation.
The subduction interface is modelled as a complex fault as defined by depth contours in the EDSF. Seismicity on the interface is concentrated between 20 and 50 km depth for the Hellenic and 10-50 km depth for the Cyprus Arcs. The surface projection of the subduction interface is shown in Fig. 14, coloured according to the assigned maximum magnitudes. Activity rates are taken to be uniform across the complex structure. In-slab seismicity is modelled in depth volumes located at 40-160 km depth in all subduction zones. In both instances seismicity is homogeneously distributed throughout the area. Maximum magnitudes as well as the b-value are the same for subduction interface and in-slab seismicity. Maximum magnitudes are estimated to range between M W = 8.3 and M W = 8.9 with  increments of 0.2 and decreasing weights (0.5, 0.2, 0.2, 0.1). The estimated b-value is 0.9 with no difference for in-slab and interface events. The SEIFA-model uses the same structure for interface and in-slab events. However, in contrast to the model used in combination with the crustal AS-and FSBG-model, the SEIFA-model uses a heterogeneous distribution of the activity rate parameters determined from a combined probability density function determined from smoothing the earthquake density and the moment release on the subduction interface (see Hiemer et al. 2014 for details). Within this model, both the maximum magnitude (M W = 9.0) and the b-value (b = 0.95) are the same for subduction interface and in-slab seismicity.
In the Vrancea region, situated beneath the southern Carpathian Arc, the most hazardous seismicity falls in the depth range 70-150 km with earthquakes occurring in the intervals 70-90 and 110-150 km. Lacking a well-defined structural description of the seismogenic processes of the area, we model the deep seismicity with four volumes at depth. The shape differs from the delineation of the crustal area sources (Fig. 14). This specific model is used in combination with all crustal seismicity models.
Earthquake density, depth levels and maximum magnitudes vary within these zones as well as the parameterisation of the preferred rupture orientation. The maximum magnitude expected for Vrancea earthquakes was assigned based on observed seismicity. We allowed for larger magnitudes by adding four magnitude increments with decreasing weights as we did for crustal sources. The M max values in the 70-90 km depth range are 7.5 and 7.8; the seismicity between 110 and 150 km is limited to M max values of 7.8-8.1 (Fig. 14).
Attenuation of seismic wave amplitudes exhibits azimuthal dependence-a characteristic that is likely due to the rupture characteristics of the deep events and wave propagation effects rather than amplification effects. This characteristic has been implemented in intensity attenuation relations (Sokolov et al. 2009;Sørensen et al. 2010) but is not considered in GMPEs that consider several spectral ordinates, as required for a harmonised European assessment. To overcome this limitation the deep seismicity is modelled to mimic the azimuthal dependence with a set of sources at depth, a scheme similar to that adopted by Sokolov et al. (2009).

Ground motion model
The logic-tree for GMPEs was developed using a new approach that combines the elicitation of key European experts with objective data testing through an iterative procedure that attempts to capture the full distribution of possible ground motion expected from future earthquakes (Delavaud et al. 2012). The data testing phase involved the tectonic regionalisation of Europe based on geological grounds and the compilation of a comprehensive strong-motion waveform data set, ultimately leading to the selection of fourteen GMPEs. These characterise the expected ground motions for all magnitude, depth and distance ranges and for all geological conditions of the source volume that exist in Europe.
Data-driven guidance provided valuable information about the ability of GMPEs to predict the ground motion in different regions based on the strong ground motion database compiled in the SHARE project ( Fig. 15; Yenier et al. 2010). The time-histories comprising the database were recorded around the world, including Europe, Middle East, California and Japan. This database (13,500 records from 2268 events recorded at 3708 stations) carefully considered the removal of duplicated entries in the event, station and waveform information through a hierarchical approach.
The strategy to select the GMPEs involved six steps: 1. Pre-selection of GMPEs suitable for the adopted tectonic regionalisation (Cotton et al. 2006;Bommer et al. 2010). 2. Expert elicitation guided by a set of predefined rules, aimed at evaluating each expert's confidence in the ability of the candidate GMPEs to predict earthquake ground motions in any given tectonic regime. 3. Objective testing of candidate GMPEs against aggregated data following the procedures described by Scherbaum et al. (2009). 4. Combination of the results of the independent steps 2 and 3 into a proposed weighting scheme. 5. Sensitivity assessment of the proposed weighting schemes of the logic tree per tectonic regime and revision of the proposed weights, if needed. 6. Definition of the final weights.
The entire procedure involved the elicitation of consortium members and key European experts on ground motion modelling (Fig. 1). The final weights to be used in the logic tree are based on expert opinion and sensitivity analyses and reflect a weighting scheme that all experts could agree on for the various tectonic regimes. The logic-tree for GMPEs is illustrated in Fig. 6.
We finally selected the following GMPEs: four for active shallow and oceanic crust (Zhao 2006;Cauzzi and Faccioli 2008;Chiou and Youngs 2008;Akkar and Bommer 2010) (Campbell 2003;Cauzzi and Faccioli 2008;Spudich and Chiou 2008;Akkar and Bommer 2010); the model developed by Toro (2002, unpublished) was selected as it is an update of that by Toro et al. (1997). For in-slab and interface subduction earthquakes we selected four models (Youngs et al. 1997;Atkinson and Boore 2003;Zhao 2006;Lin and Lee 2008), and one for volcanic and swarm type areas (Faccioli et al. 2010). The models by Toro (2002, unpublished) and Campbell (2003) were adjusted for the generic European rock site conditions established within the SHARE project, described with a shear wave velocity v rock = 800 m/s and a kappa value κ = 0.03 following Van Houtte et al. (2011). Based on a further sensitivity analysis, we modified the logic-tree proposed by Delavaud et al. (2012) and used only two ground motion models for the seismicity of the deep sources in the Vrancea region (Youngs et al. 1997;Lin and Lee 2008).

Hazard contributions of individual earthquake source models
In principle each earthquake source model, combined with the full ground motion model, represents a viable seismic hazard assessment scheme for the European region. However, each model is based on different assumptions, treats the data in different ways and reflects uncertainties differently. Delivering a full and reliable portrait of the ground shaking hazard expected from future earthquakes hence requires a well balanced combination of all models that takes into account the assumptions, the strengths and the weaknesses of each data set.

Properties of the source models
The PGA values calculated from individual source models for a 10 % exceedance probability in 50 years highlight the impact of the different assumptions underlying each model (Fig. 16). The AS-model displays relatively smooth transitions across the full range of hazard values, showing the highest values along the main fault systems of Turkey, along the coasts of Greece, along the entire length of the Italian Apennines and across Iceland. Some high hazard spots are spread across Europe, such as in the Lisbon area, in the western Pyrenees and in the Swabian Alb (Fig. 16a). As a consequence of the zonation approach, some small zones delineate hazard regions that experienced moderate-to-large historic events with relatively sharp transitions. The SEIFA-model concentrates hazard where seismicity occurred in the past, a consequence of the model assumptions (Fig. 16b). By using an adaptive kernel and considering the accumulated fault moment for spatially distributing seismicity rates, the SEIFA-model achieves a smoother representation of seismic hazard in comparison to smoothed seismicity models that use fixed kernel widths (e.g. Frankel 1995). The width of the adaptive kernels in the SEIFA-model have been optimized through pseudo-prospective likelihood tests for 5 and 10 years intervals, and hence they are more representative for these periods. Using longer optimisation periods was not possible due to the scarcity of earthquake data. Due to clustering of seismicity around some rare historical earthquakes, a few localized features remain in regions with relatively low seismic activity. The hazard computed with this model exhibits larger variations and also sharp transitions. An interesting feature is that the ground-shaking hazard varies along known fault systems such as the North Anatolian Fault Zone-a feature that may convey additional information on exceedance probabilities far larger than 10 % in 50 years. The FSBG-model concentrates its contribution to the northern Mediterranean area and the Lower Rhine Embayment. This model predicts the highest ground shaking hazard along the fast-slipping fault systems of Turkey, whereas the slower-slipping faults of south-eastern Europe and Italy return moderate-to-low hazard estimates. Hazard transitions along the slow-slipping faults are smooth; stronger gradients are found along the North Anatolian Fault system. These regional differences are mainly a consequence of the rate modelling procedure: while earthquake rates calculated for the fastest strike-slip faults may exceed the observed seismicity by a factor of up to ten, the rates calculated for slower faults-such as the normal faults of the Appennines-better match the observed rates. The comparison of the overall forecast activity (Figs. 8,9,12) shows clearly that the expected rates vary regionally and that the fault model for the Appennines forecasts smaller rates, and thus smaller hazard.

Assigning weights to individual source models
The modelling team and the consortium members agreed on weights for the individual models by judging their assumptions and main ingredients. The weighting scheme varies with the probability of exceedance levels ( Fig. 6; Table 1). Note that the combination of weights 0.5/0.3/0.2 (AS/SEIFA/FSBG) is the preferred weighting scheme for exceedance probabilities between 2 and 10 % in 50 years. This weighting scheme is used for aggregating the seismic hazard curves that are released on the website (www.efehr.org), while maps are aggregated as appropriate according to the scheme in Table 1. The rationale of using different weighting schemes originates from the assumptions of the methods applied to calculate activity rates. Using slip rates to calculate activity rates in the FSBG-model is The AS-model received the largest weight for all exceedance probability levels. Its weight increases slightly for increasing return periods, due to its qualitatively good match to the observed data and to the complex elaboration and elicitation procedure. This included the largest community-based effort aimed at combining local knowledge and expertise from all over the continent. The model incorporates information from tectonic and geological features, reflects the observed seismicity and includes specific local knowledge not present in the other models in such detail. It also represents the most widely discussed and reviewed part of the model.
In contrast, the weight of the SEIFA-model decreases for increasing return periods. This is motivated by the fact that the kernels were optimized for periods of 5 and 10 years in pseudo-retrospective likelihood-based tests (Hiemer et al. 2014). The distribution of events from the instrumental earthquake record plays a major role here. The SEIFA-model performs significantly better in the likelihood-based consistency testing suite of the Collaboratory for the Study of Earthquake Predictability (CSEP) compared to the ASmodel in 5-and 10-year pseudo-prospective tests; thus we believe it is justified to vary the weights of the SEIFA, starting with higher ones for the shorter return periods. This is done in recognition of the fact that optimizing and testing is performed for intervals in the order of 10 or less years; in other words, the implications of the 5-and 10-year pseudoprospective testing period may convey only limited information content, whereas the ASmodel is designed to deliver reliable information over much longer intervals.
Finally, the weight of the FSBG-model increases with respect to the other models for increasing return periods because the earthquake rates are estimated from geological slip rates that are likely to be more representative of long-term activity. The activity rates of the FSBG-model represent an upper limit of what is currently expected as we use the maximum slip rate for each fault source and we assume all energy to be released seismically.
Thus it is expected that the FSBG-model predicts seismicity rates that are higher than those observed in the historical record.

Seismic hazard results
PGA is one of the most common parameters used for engineering design purposes. We illustrate the spatial distribution of the PGA-values aggregated for the preferred weighting scheme for a 10 and 2 % probability of exceedance in 50 years, corresponding to return periods of 475 and 2475 years (Figs. 17, 18).

Results for the 10 % exceedance probability in 50 years
The ESHM13 map for the mean PGA at the 10 % probability of exceedance in 50 years ( Figure 19 shows the quantile estimates (5, 15, 85 and 95 %) of the hazard for the 10 % exceedance probability in 50 years, illustrating the range of uncertainties resulting from our approach (Fig. 19).

Site-specific hazard results
Hazard maps are aggregated from site-specific calculations. For each site we computed hazard curves for spectral periods (from PGA up to 10 s) and the associated uncertainties. This results in the full range of information needed for engineering applications at each site. In the following we show a few examples from cities across Europe: the reader is encouraged to refer to the online resource (http://www.efehr.org/) to explore more sites.
We display hazard curves and UHS for the cities of Rhodes, L'Aquila, Lisbon and Cologne (Fig. 20, 21; exact locations listed in Table 2). These cities lie in different tectonic settings, their activity rates are on quite different and the estimated ground shaking levels originate from varying combinations of ground-motion prediction equations. The city of Rhodes is the capital of the largest Dodecanese island, located northeast of Crete atop the Hellenic arc. L'Aquila lies in the central Appennines next to large and very active normal faults. Lisbon, located on the Atlantic shore of Portugal, suffered from the 1 November 1755 earthquake, the largest in the SHEEC catalogue (M W = 8.5 ± 0.3); in spite of its size, however, the exact location of its causative fault is still disputed (Fonseca 2006). Cologne is situated in the Lower Rhine Embayment and experiences ground shaking only from a set of well-identified active normal faults (Vanneste et al. 2013). Figures 20 and 21 illustrate the variability of the expected ground motions, showing the mean hazard level (black line) and the uncertainties of four quantiles (5, 15, 85, 95 %). We prefer to show quantiles instead of standard deviations as the underlying distributions are not necessarily Gaussian.
UHS have been computed for a range of periods of vibration (PGA to 4 s) and for probabilities of exceedance ranging from 1 to 50 % in 50 years (corresponding to return periods from 72 to 4975 years). Figure 21 shows the mean, 15th and 85th percentile UHS for a 10 % probability of exceedance in 50 years for the four selected cities. The Eurocode 8 spectral shapes applied to the locations in these cities, anchored to the SHARE mean PGA-values, are also plotted for comparison (dashed-dotted lines). The Type II spectral shape, currently used in areas where the hazard is controlled by earthquakes smaller than M W = 5.5, is seen to match fairly well the mean UHS in Cologne, but in the other cities the mean UHS is seen to be much lower than the Eurocode 8 Type I design spectrum. This difference is believed to be at least partly due to significant improvements in groundmotion prediction equations since the Eurocode spectral shapes were first proposed in the 1990s. Hence, wherever a comparison of ESHM13 maps with previous hazard studies shows an increase in the level of PGA, it is important to consider that at higher periods of vibration the ground motions from the ESHM13 UHS might actually be lower than those stipulated in Eurocode 8. The reader may refer to Weatherill et al. (2013) for further discussion on the variation of spectral shapes from ESHM13 across Europe.
Another interesting feature of the ESHM13 is the impact of epistemic uncertainty on the UHS. As shown by Figs. 20 and 21, this uncertainty in the hazard has a large influence on b Fig. 22 Uniform hazard spectra for multiple locations with varying mean PGA hazard levels (10 % EP in 50 years). Sites are selected for a a low hazard value of PGA ≈ 0.1 g (Cologne (Germany), Murcia (Spain), Haugesund (Norway), Vogtland (Germany)), b a medium hazard value of PGA ≈ 0. the spectra. Despite this influence, the modelling of epistemic uncertainty in seismic design has not received much attention to date. The authors feel that the possibility of using the epistemic uncertainty at a given return period to constrain the ground motions for more important structures such as hospitals and schools deserves further consideration, potentially by calibrating importance factors using the aforementioned epistemic uncertainty. Figure 22 shows the UHS for a number of locations with low (PGA = 0.1 g), medium (PGA = 0.2 g) and high (PGA = 0.32 g) levels of hazard. The spectral shapes are surprisingly similar, with the only significant difference found in the spectral shape of Bucharest in the high hazard plot. Differences in spectral shapes are expected when scenarios with varying earthquake magnitude are found to control the hazard at each spectral ordinate. Therefore, in order to explain the similarities and differences in the spectral shapes it would be necessary to disaggregate the hazard at each location and at each spectral ordinate, a yet outstanding task that is also necessary for the selection of appropriate accelerograms for seismic design and assessment.
Changes in the shape of uniform hazard spectra are also expected for different return periods. Figure 23 shows the UHS for a number of return periods for sites in three cities in Europe. The spectral shape clearly varies with the return period, and as expected this variation is location-specific. The current assumption in Eurocode 8 that a given spectral ordinate can be scaled to another return period of interest using a single fixed factor is thus not reasonable. Weatherill et al. (2013) have shown how this factor varies for different spectral ordinates across Europe.

Discussion and outlook
The ESHM13 delivered by the SHARE project is the first regional and community-based effort that was completed since the conclusion of the Global Seismic Hazard Assessment Program and of the SESAME program (GSHAP, Giardini 1999;Jimenez et al. 2001). The ESHM13 delivers a complete set of harmonised seismic hazard results and associated uncertainties. It may serve as a reference for updating seismic regulations at national and regional scale in Europe and may be used as a template for similar applications worldwide.
The ESHM13 represents a significant improvement in comparison to previous efforts thanks to (1) the compilation of harmonised databases of all parameters required for PSHA, (2) the adoption of rigorous, standardized procedures in all steps of the process, especially in eliciting opinions and building consensus from hundreds of involved European experts, (3) the multi-disciplinary input from all branches of earthquake science and engineering, (4) the direct involvement of the CEN/TC250/SC8 committee in defining output specifications relevant for Eurocode 8, (5) the accounting of epistemic uncertainties for both model components and hazard results, and (6) the efforts to transparent documentation and open availability of all data, results and methods.
Rigorous statistical validation of shaking hazard estimates is not possible to date (e.g. Beauval et al. 2008); similarly, it is impossible to test for maximum magnitude (Holschneider et al. 2014). However, testing of model components was performed for the source model earthquake rate forecasts (Hiemer et al. 2014). Further prospective experiments in the framework of the GEM testing facility are necessary to understand the longterm performance of the ESHM13 predictions. We maintain that such component testing is a suitable tool for quality assurance in the construction of probabilistic seismic hazard models.
Based on the current model and on the experience gained during its development, we formulate some general recommendations that might be considered in future projects.
The variety of seismotectonic environments and the combination of previous national and regional zonation models were major challenges for the construction of the area source model (AS-model). More stringent criteria to delineate area sources with a clear quality ranking of the basic data would be an additional step towards harmonisation (e.g. Wiemer et al. 2009a;Vilanova et al. 2014). Despite the effort made, the zonation model still reflects the inherent heterogeneity of the national models because the zonation model was built simultaneously with the compilation of the earthquake and fault data sets, and only revised afterwards.
We used the concept of superzones to initially assess parameter values, which were then assigned to the three earthquake source models. We specifically developed a tectonic regionalisation scheme for the entire area and used it for multiple purposes, including the homogeneous assessment of maximum magnitudes and the selection of GMPEs and fault scaling relations. While a selection of a GMPE based on the hypocentral location is natural, future work may incorporate more detailed information as suggested by .
Determining activity rates from sparse data is another major challenge. We chose to assess them using the PML-approach as this method has been developed for low-activity regions using a Bayesian approach (Johnston et al. 1994;Coppersmith et al. 2012). Besides supporting further initiatives to improve the assessment of historical earthquakes, it is similarly valuable to use small magnitude events from high quality networks for understanding the spatial and temporal variability of activity rate parameters. The depth distribution of earthquakes can have a strong influence on the resulting hazard, particularly when differentiating between processes in subduction zones or the type of faulting for crustal events. The ongoing technological advancements in seismic monitoring are expected to improve the understanding of earthquake depth distributions, thus leading to crucial improvements in the hazard calculations.
The EDSF is the largest compilation of its kind in Europe and represents a milestone for the next generation of hazard models. SHARE promoted or reinforced several regional initiatives to improve the mapping of seismogenic faults. However, most well known seismogenic faults are those located in the more active regions; those located in the continental interiors are often poorly known with only a few exceptions such as southeastern Iberia, the Lower Rhine Embayment, and the northern Balkans. Extending the mapping of seismogenic faults to the rest of Europe is a challenging task that will definitely improve future hazard models (England and Jackson 2011). Fault source parameterisation still requires an effort toward a full standardisation. A continuous update of the seismogenic fault database will contribute to reducing and understanding the incompleteness of the current release and the overall understanding of the fault network. Improving the geometric characterisation of faults, for example, works towards reducing the uncertainty in fault-site distance metrics when applying GMPEs. A more systematic collection of geological slip rates will provide better constraints on recurrence models. In addition, the current recurrence model, like in most other seismic hazard efforts, adopts the conservative assumption of full seismic coupling, thus neglecting observations of tremor, slow slip, aseismic slip and afterslip as locally observed in many parts of Europe and worldwide.
The preparation of the GMPE logic-tree defines a novel procedure that may serve as a guideline for similar projects (Delavaud et al. 2012). The selection of GMPEs, however, bears some shortcomings that should be addressed in future efforts.
Firstly, there is a lack of models for relatively rare tectonic regimes; for instance, data and models for deep seismicity that are the last glimpses on subduction processes would improve ground motion modelling in regions like Vrancea or the Betics. In particular, the azimuthal dependence of ground shaking is difficult to implement for rock hazard and also difficult to resolve when intensities are not being used. The selected ground-motion models have been calibrated using large databases that may mix data from different regions. A major challenge is thus to analyse the regional variations of earthquake properties and ground motion. This challenge is interdisciplinary as seismology, structural geology and tomographic imaging may provide independent evidence for these regional variations. Such analysis is particularly important for the continental part of Europe for which five models-two from Eastern North America and three from active shallow crustal regionswere taken into account to capture the uncertainty of regional variations of ground-motion.
Secondly, all GMPEs feature limited earthquake magnitude and fault-site distance ranges; their extrapolation to calculate hazard starting from lower magnitudes or using inappropriate distances is critical.
Thirdly, the selected GMPEs are in many ways similar as all of them are regression models. Thus it is worth investigating physics-based simulation strategies as a complementary pathway for ground motion modelling, for improving the databases which will be needed to test and validate these models, and for investigating alternative approaches such as the use of stochastic GMPEs (e.g. Edwards and Fäh 2013). Indeed, the GMPE testing activities conducted within SHARE have shown the need to update the metadata information and develop a uniformly European processed accelerometric data archive (Akkar et al. 2013).
Furthermore, the project was limited to cover on-shore regions as GMPEs do not exist for the continental shelf. The interest for these regions is rising quickly as the level of their economic exploitation increases, but due to the lack of calibration of the GMPEs the understanding of seismic hazard and risk of these important areas remains rudimentary.
We focused primarily on the development of a seismic hazard map for an average shear wave rock velocity of v S30 = 800 m/s and κ = 0.03 s. We investigated also the definition of a viable proxy to provide site amplification values for the entire European region, i.e. the correlation between topographic slope and v S30 , but the conclusion is that this approach is not fully applicable (Lemoine et al. 2012). We conclude that future projects on this scale need to investigate this finding in better detail.
The ESHM13 results need to be evaluated in detail from an engineering perspective, however, this is outside the scope of this overview. Nevertheless, the main conclusion that can currently be drawn from this limited analysis is that the availability of hazard curves, UHS and disaggregation results via an interactive web portal will allow design codes to make more direct use of PSHA results, having to rely less on assumptions and simplifications.
Probabilistic seismic hazard assessment is a community effort which requires a welltested expert elicitation procedure. The structure of the project included multiple review procedures, for the data as well as for modelling. It is however desirable that in future projects guidance and formal criteria be defined in advance and followed considering the highest level guidelines (Krammerer and Ake 2011). We focused on transparent communication, open access to data and the use of published approaches. Data and results are available from the European Facility for Earthquake Hazard and Risk at www.efehr.org. Hazard was computed with the GEM OpenQuake engine with the intention that results be and remain reproducible (Pagani et al. 2014). We are confident that this effort will leave a legacy of scientific rigor in PSHA, both at the European level and worldwide. The ESHM13 is a starting point for harmonising hazard as input to Eurocode 8 and we expect that further efforts at the national scale will help improving seismic hazard and risk mitigation all over Europe.

ESHM13 data and resources
The resulting ESHM13 provides more than five hundred maps displaying the ground shaking expected to be reached or exceeded with probabilities of exceedance in the range 1-50 % in 50 years and for a range of periods of ground acceleration (from PGA to 10 s).
Hazard curves, UHS and hazard disaggregation schemes are available for more than 120,000 sites equally spaced at 10 km across Europe and Turkey. Mean, median and quantile results are produced for a reference rock condition (available online), i.e. Eurocode 8 Type A with an average shear wave velocity v s30 = 800 m/s and κ = 0.03 s. Hazard curves are computed for return periods up to 10,000 years, however they must be used with much caution because the inclusion of very low-activity faults that have not entered the SHARE model due to its regional scope, or insufficient uncertainty descriptions in ground motion models may affect the results. We therefore limit the representation of results for hazard maps to 5000 years.
All data, results, references and print material, including the official poster of the reference hazard map (Giardini et al. 2013b), are freely accessible online at http://www. efehr.org, the portal of the European Facility for Earthquake Hazard and Risk, and at the SHARE project website (www.share-eu.org).
The EDSF is made available to the public through the web site (http://diss.rm.ingv.it/ share-edsf/), hosted and maintained by Istituto Nazionale di Geofisica e Vulcanologia (INGV), which features map navigation, tabulated-data browsing, as well as seamless download of data files in various GIS formats.
To date more than 30 peer-reviewed scientific manuscripts have been published as a result of the project; they are listed at www.share-eu.org.
The ESHM13 model was generated through the community effort of over 300 scientists and engineers contributing their data and expertise. It is the first regional contribution to the Global Earthquake Model (GEM, www.globalquakemodel.org) initiative. SHARE adopted the OpenQuake hazard engine ver. 1.0 (Pagani et al. 2014) to compute the ESHM13.