# The 2020 update of the European Seismic Hazard Model - ESHM20: Model Overview

<figure><img src="image_1.png" /></figure>

Laurentiu Danciu, Shyam Nandan, Celso Reyes, Roberto Basili, Graeme Weatherill, Celine Beauval, Andrea Rovida, Susana Vilanova, Karin Sesetyan, Pierre-Yves. Bard, Fabrice Cotton, Stefan Wiemer, Domenico Giardini
# The 2020 update of the European Seismic Hazard Model - ESHM20: Model Overview

Laurentiu Danciu¹, Shyam Nandan¹, Celso Reyes¹, Roberto Basili², Graeme Weatherill³, Céline Beauval⁴, Andrea Rovida⁵, Susana Vilanova⁶, Karin Sesetyan⁷, Pierre-Yves Bard⁴, Fabrice Cotton³, Stefan Wiemer¹, Domenico Giardini⁸

1. Swiss Seismological Service, ETH Zurich, Sonneggstrasse 5, 8092 Zurich, Switzerland,

2. Istituto Nazionale di Geofisica e Vulcanologia (INGV), 00143 Rome, Italy,

3. German Research Centre for Geosciences (GFZ), Section 2.6 Seismic Hazard & Risk Dynamics, 14473 Potsdam, Germany,

4. Institut des Sciences de la Terre (ISTerre), IRD, 38058 Grenoble, France,

5. Istituto Nazionale di Geofisica e Vulcanologia (INGV), 20133 Milan, Italy,

6. Instituto Superior Tecnico (IST), 1049-001 Lisboa, Portugal,

7. Bogazici University, Kandilli Observatory and Earthquake Research Institute, Department of Earthquake Engineering, 34684, Cengelkoy, Istanbul, Turkey,

8. Institute of Geophysics, ETH Zurich, Sonneggstrasse 5, 8092 Zurich, Switzerland
## Rights and permissions

Copyright © 2021 ETH Zurich

Except where otherwise noted, this work is licensed under a Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses/by/4.0/).

This report aims at documenting the 2020 update of the European Seismic Hazard Model. The findings, comments, statements or recommendations expressed in this report are exclusively of the author(s) and not necessarily reflect the views and policies of the institutions listed here: i.e. i.e. Swiss Seismological Service, ETH Zurich, Istituto Nazionale di Geofisica e Vulcanologia (INGV), German Research Centre for Geosciences (GFZ), Institut des Sciences de la Terre (ISTerre), Instituto Superior Tecnico (IST), Bogazici University, Kandilli Observatory and Earthquake Research Institute, Department of Earthquake Engineering, the EFEHR Consortium or the European Union. The core team of ESHM20 has tried to make the information in this product as accurate as possible. However, it does not guarantee that the information is totally accurate or complete. Therefore, you should not solely rely on this information when making a commercial decision. Users of information expressed in this report assume all liability arising from such use. While undertaking to provide practical and accurate information, the authors assume no liability for, nor express or imply any warranty with regard to the information contained hereafter.

EFEHR Consortium is committed to providing web accessible content and results of the ESHM20 via the www.efehr.org platform. If you are having difficulties with accessing this document please email efehr@sed.ethz.ch

## Citation Advice

Danciu L., Nandan S., Reyes C., Basili R., Weatherill G., Beauval C., Rovida A., Vilanova S., Sesetyan K., Bard P-Y., Cotton F., Wiemer S., Giardini D. (2021) - The 2020 update of the European Seismic Hazard Model: Model Overview. EFEHR Technical Report 001, v1.0.0, https://doi.org/10.12686/a15
# Table of contents

<table><tr><td>Acknowledgments</td><td>8</td></tr><tr><td>Executive Summary</td><td>9</td></tr><tr><td>1 Introduction</td><td>10</td></tr><tr><td>1.1 Framework and Model Building Philosophy</td><td>10</td></tr><tr><td>1.2 Outline of Report</td><td>12</td></tr><tr><td>1 Main Input Datasets</td><td>14</td></tr><tr><td>1.1 Historical Earthquake Catalogue</td><td>14</td></tr><tr><td>1.1.1 Background and compilation strategy</td><td>14</td></tr><tr><td>1.1.2 Input data</td><td>14</td></tr><tr><td>1.1.3 Earthquake parameters</td><td>15</td></tr><tr><td>1.1.4 EPICA, the European PreInstrumental earthquake CAtalogue</td><td>17</td></tr><tr><td>1.2 Instrumental Earthquake Catalogue</td><td>18</td></tr><tr><td>1.2.1 Data collection and regionalisation</td><td>19</td></tr><tr><td>1.2.2 Harmonisation</td><td>20</td></tr><tr><td>1.3 Active Crustal Faults and Subduction Zones</td><td>21</td></tr><tr><td>1.3.1 Crustal faults</td><td>21</td></tr><tr><td>1.3.2 Subduction Zones</td><td>27</td></tr><tr><td>1.4 Tectonic Regionalization</td><td>29</td></tr><tr><td>2 Unified Earthquake Catalogue</td><td>31</td></tr><tr><td>2.1 Catalogue Declustering</td><td>32</td></tr><tr><td>2.1.1 Sensitivity Analysis</td><td>33</td></tr><tr><td>2.2 Catalogue Completeness</td><td>35</td></tr><tr><td>2.2.1 Completeness zones (CSZs)</td><td>37</td></tr><tr><td>2.2.2 Improved temporal course of earthquake frequency (TCEF) method</td><td>39</td></tr><tr><td>2.2.3 Time Series of Mc</td><td>40</td></tr><tr><td>2.2.4 Objective quality assessment of time series of Mc</td><td>40</td></tr><tr><td>3 Seismogenic Source Model</td><td>43</td></tr><tr><td>3.1 Unified Area Sources Model</td><td>43</td></tr><tr><td>3.1.1 The use of Super Zones</td><td>44</td></tr><tr><td>3.1.2 Activity Rates Parameters and Magnitude Frequency Models</td><td>45</td></tr><tr><td>3.2 Active faults and background seismicity model</td><td>51</td></tr><tr><td>3.2.1 Active Faults</td><td>51</td></tr><tr><td>3.2.2 Smoothed Seismicity</td><td>53</td></tr><tr><td>3.2.3 Combining the Smoothed Seismicity with the Active Faults</td><td>54</td></tr><tr><td>3.3 Subduction and Deep Seismicity Sources</td><td>56</td></tr><tr><td>3.3.1 Subduction Interface</td><td>56</td></tr><tr><td>3.3.2 Deep Seismicity: Vrancea, Romania</td><td>59</td></tr><tr><td>3.4 Seismogenic Source Model - Logic Tree</td><td>60</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>4</td>
      <td>Ground Motion Model Characterization</td>
      <td>62</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Ground Motion Modelling Approach</td>
      <td>62</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Ground Motion Model Logic Tree for Shallow Seismicity</td>
      <td>62</td>
    </tr>
    <tr>
      <td>4.2.1</td>
      <td>A regionally adjustable backbone Ground Motion Model (GMM) and logic tree</td>
      <td>62</td>
    </tr>
    <tr>
      <td>4.2.2</td>
      <td>Aleatory Uncertainty</td>
      <td>67</td>
    </tr>
    <tr>
      <td>4.2.3</td>
      <td>Site Amplification and Characterisation</td>
      <td>68</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Scaled Backbone GMM Logic Tree for the Stable Craton</td>
      <td>70</td>
    </tr>
    <tr>
      <td>4.3.1</td>
      <td>Defining and delineating the cratonic region of Europe</td>
      <td>70</td>
    </tr>
    <tr>
      <td>4.3.2</td>
      <td>Calibrating the scaled backbone GMM for median ground motions</td>
      <td>72</td>
    </tr>
    <tr>
      <td>4.3.3</td>
      <td>Site Amplification</td>
      <td>73</td>
    </tr>
    <tr>
      <td>4.3.4</td>
      <td>Aleatory variability</td>
      <td>73</td>
    </tr>
    <tr>
      <td>4.3.5</td>
      <td>Complete Logic Tree for GMM in the Stable Cratonic Region of NE Europe</td>
      <td>74</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Subduction and Deep Seismicity</td>
      <td>75</td>
    </tr>
    <tr>
      <td>4.4.1</td>
      <td>Identifying "non-shallow crustal earthquakes"</td>
      <td>75</td>
    </tr>
    <tr>
      <td>4.4.2</td>
      <td>Selecting and modifying a backbone deep/subduction seismicity GMM</td>
      <td>76</td>
    </tr>
    <tr>
      <td>4.4.3</td>
      <td>Quantifying the epistemic uncertainty and constructing the logic tree</td>
      <td>77</td>
    </tr>
    <tr>
      <td>4.4.4</td>
      <td>Modifying the forearc/backarc scaling term</td>
      <td>79</td>
    </tr>
    <tr>
      <td>4.5</td>
      <td>Special cases and adjustments</td>
      <td>80</td>
    </tr>
    <tr>
      <td>4.5.1</td>
      <td>North-western Europe</td>
      <td>81</td>
    </tr>
    <tr>
      <td>4.5.2</td>
      <td>Iceland</td>
      <td>83</td>
    </tr>
    <tr>
      <td>5</td>
      <td>Model Implementation and Hazard Calculation</td>
      <td>85</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Seismogenic Source Model - implementation</td>
      <td>85</td>
    </tr>
    <tr>
      <td>5.1.1</td>
      <td>Area Source Zones</td>
      <td>86</td>
    </tr>
    <tr>
      <td>5.1.2</td>
      <td>Active Faults and Subduction Interface</td>
      <td>86</td>
    </tr>
    <tr>
      <td>5.1.3</td>
      <td>Logic Tree Seismogenic Sources and OpenQuake Input Files</td>
      <td>87</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Ground Motion Model - implementation</td>
      <td>89</td>
    </tr>
    <tr>
      <td>5.2.1</td>
      <td>Implementation of the Ground Motion Models</td>
      <td>89</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>OpenQuake Calculation Settings</td>
      <td>91</td>
    </tr>
    <tr>
      <td>6</td>
      <td>Main Results</td>
      <td>93</td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>Seismic Hazard Maps</td>
      <td>93</td>
    </tr>
    <tr>
      <td>6.1.1</td>
      <td>475yrs ground shaking hazard maps</td>
      <td>93</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Comparison with ESHM13</td>
      <td>94</td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>Hazard Contribution per Source Model</td>
      <td>98</td>
    </tr>
    <tr>
      <td>6.4</td>
      <td>Site Specific Hazard Estimates</td>
      <td>98</td>
    </tr>
    <tr>
      <td>6.4.1</td>
      <td>Uniform Hazard Spectra</td>
      <td>98</td>
    </tr>
    <tr>
      <td>6.4.2</td>
      <td>Hazard Curves</td>
      <td>99</td>
    </tr>
    <tr>
      <td>7</td>
      <td>Data, Products and Resources</td>
      <td>101</td>
    </tr>
    <tr>
      <td>7.1</td>
      <td>Open Data and Licensing</td>
      <td>101</td>
    </tr>
    <tr>
      <td>7.2</td>
      <td>OQ input Files</td>
      <td>101</td>
    </tr>
    <tr>
      <td>7.3</td>
      <td>Online Resources</td>
      <td>102</td>
    </tr>
    <tr>
      <td colspan="2">References</td>
      <td>105</td>
    </tr>
    <tr>
      <td colspan="2">Appendix A: Earthquake Catalogue Completeness - Summary</td>
      <td>111</td>
    </tr>
  </tbody>
</table>
<table><tr><td>Appendix B: ESHM20 Input Datasets: Summary and File Formats</td><td>113</td></tr><tr><td>B.1 Catalogue file:</td><td>114</td></tr><tr><td>B2 Completeness Superzones (CSZ) shape file</td><td>115</td></tr><tr><td>B.3 Tectonic Zone (TECTO) shape files</td><td>116</td></tr><tr><td>B.4 Area Source Zone (ASZ) shape files</td><td>117</td></tr><tr><td>B.5 Active fault source model shapefile</td><td>118</td></tr><tr><td>B.6 Background seismicity SEIFA shapefile</td><td>119</td></tr><tr><td>B.7 Region Filter shapefile</td><td>120</td></tr><tr><td>B.8 Fault buffer shapefile</td><td>120</td></tr></table>
Acknowledgments

We are grateful all researchers, scientist, experts and engineers contributing to the development of the 2020 European Seismic Hazard Model (ESHM20), in various way from data compilation and curation to knowledge transfer during numerous meetings and webinars:

Andrea Antonucci, António Araújo Correia, Sinan Akkar, Kuvvet Atakan, Jure Atanackov, Hamid Sadegh Azar, Stephane Baize, Julien Barriere, Paolo Bazzuro, Maria Belen Benito, Myriam Belvaux, Bjarni Bessason, Dino Bindi, Christian Bosse, Thierry Cameelbeck, Eduardo Cansado Carvalho, Michele Carafa, Fernando Carrilho, Alexandra Carvalho, Carlo Virgilio Cauzzi, Gilles Celli, Eugenio Chioccarelli, Carmen Cioflan, Pasquale Cito, Cécile Cornou, Edward Cushing, Susana Custódio, Snjezana Cvijic-Amulic, Nicola D'Agostino, Vera D'Amico, Nicolas D'Oreye, Atefe Darzi, Mathias Dolsek, John Douglas, Branko Dragicevic, Dejan Dragojević, Stephane Drouet, Blaise Duvernay, João Estevão, Donat Fäh, Ekkehard Fehling, Joao Fonseca, Mariano Garcia-Fernandez, Eulalia Gracia, Gottfried Grünthal, Benedikt Halldorsson, Florian Haslinger, Marian Herak, Jomard Herve, Stefan Hiemer, Fabrice Hollander, Shible Hussein, Iunio Iervolinio, Diethelm Kaiser, Vanja Kastelic, Amir M. Kaynia, Rexhep Koci, Annakaisa Korja, Sreeram Kotha, Svetlana Kovacevic, Milad Kowsari, Olga-Joan Ktenidou, Daniela Kuehn, Kresimir Kuk, Neki Kuka, Pierre Labbe, Steffi Lammers, Giovani Lanzano, Björn Lund, Lucia Luzi, Francesco Maesano, Elena Manea, Päivi Mäntyniemi, Basil Margaris, Christophe Martin, Warner Marzocchi, Frederic Masson, Luís Matias, Carlo Meletti, Nikolaos Melis, Albero Michelini, Jadranka Mihaljevic, Zoran Milutinovic, José Antonio Peláez Montilla, Ilaria Mosca, Roger Musson, Shemsi Mustafa, Adrien Oth, Bruno Pace, Roberto Paolucci, Christos Papaioannou, Florin Pavel, Laura Peruzza, Rui Pinho, Kyriazis Pitilakis, Valerio Poggi, Mircea Radulian, Evi Riga, Philippe Roth, Agathe Roulle, Zafeiria Roumelioti, Radmila Salic, Abdullah Sandikkaya, Antoine Schlupp, Jochen Schwarz, Anselm Smolka, Oona Scotti, Barbara Sket-Motnikar, Vitor Silva, Anne Socquet, Efthimios Sokos, Mathilde Sorensen, Carlos Sousa Oliveira, Thomas Spies, Max Stuchi, Nikolaos Theodoulidis, Mara Tiberti, Paola Traversa, Nino Tsereteli, Radu Vacareanu, Dina Vales, Roberto Vallone, Dimitrios Vamvatsikos, Kris Vanneste, Emmanuel Viallet, Daniele Vigano, Francesco Visini, Stefan Weginger, Polona Zupancic

Special thanks to the SERA project management and logistics: Kauzar Saleh Contell and Romano Meier.

We would like to thank IT team at Swiss Seismological Service, ETH Zurich for their excellent support on maintaining the servers, update the databases and upgrade the EFEHR web-page and the *hazard.EFEHR* web-platform to access the 2020 models: Philipp Kästli, Leandra Eberle, Emil Zylis, Emilia Petronio and Cyrill Bonjour.

Furthermore, we would like to also thank Helen Crowley and Jamal al Dabeek for their efforts and feedback regarding the use of the seismic hazard for risk calculation.

We would like to extend our gratitude to Michèle Marti, Nadja Valenzuela, Irina Dallo and Simone Zaugg for their creativity and coordination of the outreach activities.

In addition, we would like to thank Marco Pagani for scientific discussions and feedback, and Michele Simionato for timely OpenQuake improvements.

Finally, the authors would wish to express their appreciation to the SC8 Ad-hoc Working Groups for the excellent collaboration during the ESHM20 evaluation: Pierre Labbe, Philippe Bisch, António Araújo Correia, Hamid Sadegh Azar, Alexandra Carvalho, Matjaž Dolšek, Blaise Duvernay, João Estevão, Ekkehard Fehling, Max Gündel, Iunio Iervolino, Christophe Martin, Roberto Paolucci, Evangelia Peli, Kyriasis Pitilakis, Jochen Schwarz, Mathilde Sorensen, Evi Riga, Kris Vanneste, Dimitrios Vamvatsikos, Emanuel Viallet, Radu Vacareanu
Executive Summary

The 2020 European Seismic Hazard Model (ESHM20) provides an update of the earthquake hazard
assessment of the Euro-Mediterranean region. ESHM20 has received funding from the European Union's
Horizon 2020 research and innovation programme under grant agreement No.s 730900 of the SERA-
Project (www.sera-eu.org)

The ESHM20 follows the same principles as the ESHM13, with state-of-the art procedures
homogeneously applied for the entire pan-European region, without the country-borders issues. The
model was built upon recently compiled datasets (i.e. earthquake catalogues, active faults, ground
shaking recordings), information (tectonic and geological) and models (seismogenic sources, ground
shaking). A fully probabilistic framework was adopted in the hazard model implementation and all
datasets and inputs are fully cross-border harmonized.

The newly developed seismogenic source model encompass fully harmonized and cross borders
seismogenic sources following the recent national earthquake hazard models. The inherent
uncertainties in characterizing the earthquake rupture forecast are handled by a complex logic tree,
consisting of two main models (branches): an area source-based model and a hybrid fault-smoothed
seismicity model. The ground motion characteristic model is built upon the most complete ground
shaking recordings in Europe and aims at capturing the effects of source and attenuation path of the
expected ground shaking at a site.

The regional variability of the ground shaking as constrain by observations is captured by a novel
approach of spatial clustering and the overall uncertainties are handled in a backbone logic tree. The
master logic tree combines the earthquake rate forecast with the backbone ground motion models into
a computational model for assessing the earthquake ground shaking at across the entire Euro-
Mediterranean region. Full sets of hazard results (i.e. hazard curves and maps, uniform hazard spectra)
for various combinations of return periods and descriptive statistics (mean, median and quantiles) are
available.

The ESHM20 development process involved several regional workshops where the scientists and experts
were consulted and their feedback was acknowledged and considered on the model outcome. Finally,
ESHM20 provides the basis to derive informative hazard maps for two key engineering parameters
defining the anchoring points of the seismic design spectra for the next version of the European Seismic
Design Code (CEN-EC8).

The European Facilities for Earthquake Hazard and Risk (EFEHR) will maintain and further develop this
model in collaboration with the GEM Foundation and the European Plate Observing System (EPOS).

The source data, input models, software and outputs of ESRM20 are thus being openly released with a
Creative Commons CC-BY International 4.0 license (https://creativecommons.org/licenses/by/4.0/).
This license allows re-users to distribute, remix, adapt, and build upon the material in any medium or
format, so long as attribution is given to the creator. The license allows for commercial use. Each product
is released with a clear notice on how it should be cited in order to abide by the license. ESHM20 is online
available at the www.hazard.EFEHR.org.
# 1 Introduction

The 2020 European Seismic Hazard Model (ESHM20) was developed within the EU founded project "Seismology and Earthquake Engineering Research Infrastructure Alliance for Europe" (SERA), under the call INFRAIA-01-2016-2017 Research Infrastructure for Earthquake Hazard in Europe. ESHM20 was initiated within the SERA-Joint Research Activities (JRA3) with the following objectives:

- Revisit the 2013 European Seismic Hazard Model (ESHM13, Wössner et al 2015) to identify the input datasets (earthquake catalogues, active faults, ground motion recordings) and the components to be updated.

- Develop and updated seismic hazard model based on technically sound methods and ensure the quantification of inherent uncertainties.

- Build upon the principles of the ESHM13: consensus and acceptance of the input models, transparency and open access of input datasets, components and output.

- Update the seismogenic source part (catalogues, faults, sources) of ESHM13 in time for the EC 8 revision, so that the 2nd generation codes are based on the most recent and most state-of-the art information.

- Revise the ESHM13 Ground Motion Prediction Equation (GMPE) logic tree in light of the newly available ground motion data.

- Interact with CEN/TC250/SC8 committee as well as with key experts at the national level to ensure the correct information and timely implementation of the ESHM20, and extend the output to serve additional engineering requirements as part of design standard EC 8 update and revision.

- Embed the ESHM20 into the Global Earthquake Mosaic (Pagani et al 2018), ensuring consistency with the overall global model.

- Enable earthquake risk modelling in Europe at regional (ESHRM20, Crowley et al 2021) scale

## 1.1 Framework and Model Building Philosophy

An open and collaborative framework was adopted for the development of the 2020 update of the European seismic hazard assessment. The work was coordinated and conducted by a core team, mainly the authors of this report, with contributions from a large community of earth scientists, seismologists, geologists, earthquake engineers, statisticians, software developers and outreach specialists (see a full list in the Acknowledgements).

The role of the collaborators is critical as they contribute in various ways: provide access to the raw data, provide access to the national hazard models, assess the quality of the datasets and models, review the science and methods, provide technical advancements, etc. Their contribution was integrated in numerous meetings including a) bilateral meetings between the core-team members and the local experts b) regional review meetings and c) community webinars.

The bilateral meetings are an excellent opportunity to interact with the local experts, discuss the scientific methods and the assumption underlying the model development, evaluate the quality and curate the main input datasets, seek consensus on the cross-border harmonization, review the main components of the results, file exchange and provide feedback. This interaction between the various actors involved at regional and national level is essential to the acceptance of the ESHM20 at the national level.

### Regional Meetings

The objectives of the regional meetings were to enhance community understanding and stimulate open discussion on the model development. These meetings were attended in person by about 150 European experts and researchers. Furthermore, the regional meetings provided the opportunity for the ESHM20
core team to present the status update on data compilation, curation and harmonization, development of the seismogenic sources, ground motion models, model implementation, outputs and results spanning across Europe without limitations of the country boundaries. Four regional review workshops were organized in Lisbon (June, 6th to 7th, 2019), Potsdam, June (12th to 13th, 2019), Athens, July (2nd to 3rd, 2019) and Pavia, (October 14th, 2019). An additional workshop was organized in Luxembourg, March 27th, 2019 with focus on the ground motion characterisation model.

## Public Webinars

A public webinar was organized on September 9, 2020 – as a scientific session of the EFEHR Consortium General Assembly. The webinar was attended online by an audience of about 250 registered participants. The feedback received during the online webinar was summarized in a general response and shared with the community. Meeting presentations, supporting materials and results made available during the review and consultation process until October 2021.

Another online webinar was organized as a special session on the 37th General Assembly of the European Seismological Commission, where all the components of ESHM20 and results were presented.

## Interaction with the CEN/TC 250/SC 8 ad-hoc working groups

An important aim of the ESHM20 is to provide the basis for deriving the engineering products for CEN/TC 250/SC 8, i.e., the subcommittee of national experts that is responsible for the development of Eurocode 8 within the European Committee for Standardization (CEN). Two *ad-hoc* working groups i.e. AHG1 with focus on site conditions and AHG2 with focus on seismic hazard evaluation.

The scope of these two *ad-hoc* groups was to analyse the definition of seismic action in EN 1998, interact with the ESHM20 core team, analyse the comments and feedback, with a view to proposing to SC 8 a consensual version of the respective clauses in order to finalise the corresponding informative annex in EN 1998-1-1. Furthermore, the ESHM20 was evaluated with different comparison and testing activities for several countries in Europe (in alphabetic order): France, Greece, Italy, Norway, Portugal, Romania, Slovenia and Switzerland.

To support this evaluation the core team provided the support and access to all materials, input files, results, documentation and presentations. Several bilateral meetings and four plenary meetings took place in the last phase of the model development cycle (2020 to 2021). At the moment when this report is drafted, although a consensus has not yet been reached on some points, on the basis of positions expressed by the (vast) majority of participants, the recommendation of the convener of AHG2 to SC 8 was that $S_\alpha$ and $S_\beta$ maps derived from the median ESHM20 map at 475 years return period are incorporated in EN 1998-1-1 (Doc N1111 of CEN/TC 250/SC 8).

After discussion during the 41st SC 8 meeting, SC 8 approved the inclusion of the ESHM20 maps in Annex A of prEN 1998-1-1:2021, which will now undergo the enquiry and formal voting procedure by CEN member countries.

## Interaction with the CEN SC8 ad-hoc working groups

An important aim of the ESHM20 is to provide the basis for deriving the engineering products for CEN/TC 250/SC 8, i.e., the subcommittee of national experts that is responsible for the development of Eurocode 8 within the European Committee for Standardization (CEN). Two *ad-hoc* working groups i.e. AHG1 with focus on site amplification and AHG2 with focus on seismic hazard evaluation.
The scope of these two ad-hoc groups was to analyse the definition of seismic action in EN 1998, interact with the ESHM20 core team, analyse the comments and feedback, with a view to proposing to SC8 a consensual version of the respective clauses in order to finalise the corresponding informative annex in EN 1998-1-1. Furthermore, the ESHM20 was evaluated with different comparison and testing activities for several countries in Europe (in alphabetic order): France, Greece, Italy, Norway, Portugal, Romania, Slovenia and Switzerland.

To support this evaluation the core team provided the support and access to all materials, input files, results, documentation and presentations. Several bilateral meetings and four plenary meetings took place in the last phase of the model development cycle (2020 to 2021). At the moment when this report is drafted, although a consensus has not yet been reached on some points, on the basis of positions expressed by the (vast) majority of participants, recommendation of the convener to SC8 is that Sα and Sβ maps derived from the median ESHM20 map at 475 years return period are incorporated in EN1998-1-1 Annex G (Pierre Labbe, lead of Sc8-AHG2, personal communication)

**Interaction with the Global Earthquake Model**

During the development phase of the ESHM20, there was an excellent collaboration with the GEM's hazard and software developing team. Model calculation was done with the OpenQuake open-source software (Pagani et al 2014) which was constantly improved during the model development cycle. ESHM20 will replace the ESHM13 with the Global Earthquake Mosaic (Pagani et al 2021).

**ESHM20 release cycle:**

The model was developed between September 2017 and November 2021. In the first 2 years, we focused on compiling and curating the main datasets, designing the blueprint of the models and approaches, software development, and defining the computational workflow. During the following year, the model was implemented iteratively, with calculations undertaken on these preliminary models. In the final year the model was extensively reviewed and quality assurance procedures were established. The model release cycle was:

- ESHM20 v08 – community preview version, released to partners, collaborators and to the hazard and risk community in September 2020

- ESHM20 v12d – beta version, distributed to partners, collaborators in February 2021

- ESHM20 v12e – final version, to be released in December 2021

## 1.2 Outline of Report

In summary, this report provides an overview of the ESHM20 with focus on the main input datasets, main components, calculation framework and results.

Following the background information and scope of the study, Chapter 2 of this report describes the main input datasets from earthquake catalogues, active faults, subduction sources and regional tectonic regionalization. We discuss the methods applied to each dataset, the harmonization strategies and procedures to curate these datasets.

Chapter 3 provides an overview of the unified earthquake catalogue, which is the foundation of any seismic hazard assessment. The newly developed completeness analysis of the unified earthquake catalogue is given also together with insights of the declustering analysis.

Chapter 4 discusses the overall seismogenic sources, reviews the model assumptions, describes the main sources i.e. the area sources, the active faults, background smoothed seismicity and the logic tree.
Chapter 5 addresses the development of the ground motion characteristic model. We provide the insights of the back-bone approach to characterize this uncertainty for seismicity in Europe, incorporating region-to-region source and attenuation variability based on European strong motion data.

Chapter 6 gives guidance on model implementation and hazard calculation settings for the OpenQuake software (Pagani et al 2014).

The main results are summarized in Chapter 7 and Chapter 8 indicates the structure of the inputs with focus on file formats and description of the main attributes. Links to the main repositories and data resources are provided also in this chapter.

The present report does not provide or repeat information on the overall probabilistic seismic hazard assessment (PSHA); rather our purpose is to offer additional, practical insights into the development and use of the ESHM20. A special issue on a peer-reviewed journal is in preparation at the time of writing this report.
# 1 Main Input Datasets

## 1.1 Historical Earthquake Catalogue

### 1.1.1 Background and compilation strategy

The 1000-1899 seismic catalogue compiled for the ESHM20 is the European PreInstrumental earthquake CAtalogue EPICA (https://doi.org/10.13127/epica.1.1; Rovida and Antonucci, 2021a), an update of SHEEC 1000-1899 (https://doi.org/10.6092/ingv.it-share-sheec.1000-1899; Stucchi et al., 2013) based on the same main principles - mostly transparency, repeatability, and continent-wide harmonization of data - and compilation strategies and methods.

As detailed in Stucchi et al. (2013), SHEEC 1000-1899 relied upon the data gathered in the European Archive of Historical Earthquake Data AHEAD (https://doi.org/10.6092/ingv.it-ahead; Albini et al. 2013; Locati et al. 2014; Rovida and Locati, 2015), including both macroseismic intensity data (hereafter MDP, Macroseismic Data Point) supplied by descriptive historical seismological studies and macroseismic databases, and parameters contained in regional catalogues.

For compiling the catalogue, such data were thoroughly analyzed in order to select the most representative of the knowledge of each earthquake, independently from national constraints. From the selected MDPs distributions, earthquake location and magnitude were determined with three methods based on regional attenuation models of macroseismic intensity with magnitude and source distance. In addition, parameters from regional catalogues, selected from AHEAD according to their reliability and public availability, were also considered. The earthquake parameters consisted of locations and magnitude selected among those calculated from MDPs and derived from regional catalogues according to a priority scheme.

The same procedures were applied in the compilation of EPICA, and they are summarized in this report and detailed in dedicated report currently in preparation (Rovida and Antonucci, 2021b). In particular, for the compilation of EPICA version 1.1 as the new version of SHEEC 1000-1899, all the updates concerning both i) input macroseismic datasets and catalogues, and ii) regional calibrations of the methods for parameterizing MDPs, or new methods, were surveyed.

However, since the publication of SHEEC 1000-1899 in 2013, only historical earthquake data have been subject to updates and innovations, whereas data that can significantly improve the calibration have not been published, as detailed in Basili et al., (2018). A new Boxer's calibration has been realized for Italy, only (Rovida et al., 2020a), and no new robust method for deriving location and magnitude from intensity data has been recently proposed.

For these reasons, the compilation of EPICA version 1.1 focused on new input datasets, whereas the parametrization procedures were maintained the same of SHEEC 1000-1899. EPICA version 1.1 contains 5703 earthquakes with either maximum intensity ≥5 or Mw ≥4.0 occurred between 1000 and 1899.

### 1.1.2 Input data

The input data of EPICA were provided by the European Archive of Historical Earthquake Data AHEAD, which collects, systematically organizes and makes available data on European earthquakes in the period 1000-1899 from i) regional and national macroseismic databases, ii) seismological descriptive studies on specific earthquakes, periods or areas, and iii) major modern parametric catalogues.

AHEAD connects multiple datasets referred to the same earthquake, enabling their analysis and comparison for the identification of the most updated and informative one. For both SHEEC 1000-1899
and EPICA, the selection considered the availability of MDPs, the features of the supporting historical seismological research and the update of each dataset, and ignored the national provenance of it. In the end, AHEAD provided EPICA version 1.1 with 160 sources of intensity data and 39 parametric catalogues (see Rovida and Antonucci, 2021b).

The 160 sources of intensity data supply 49852 MDPs related to 3622 earthquakes. These data sources include the nine regional nodes contributing to AHEAD, complemented with the results of several investigations on single or groups of earthquakes, areas and/or periods. As a result, 46 studies providing MDPs to 2014 earthquakes are not among those used for SHEEC 1000-1899, mostly because they are more recent or because they were not considered although already available, e.g. those from the Italian Archive of Historical Earthquake Data ASMI (https://doi.org/10.13127/asmi; Rovida et al., 2017), which did not yet exist at that time.

The 39 regional catalogues EPICA relies upon relate to 5511 earthquakes. For guaranteeing transparency, EPICA exclusively considers published catalogues. Following this principle and for the period of interest (1000-1899), only the following four catalogues were suited for EPICA in addition or substitution of those contributing to SHEEC 1000-1899: CPTI15 (Rovida et al., 2016; 2020a) for Italy, FCAT-17 (Manchuel et al., 2018) for France, Hammerl and Lenhardt (2013) for Lower Austria, and Leydecker (2011) for Germany (although almost identical to the previous version).

### 1.1.3 Earthquake parameters

Parameters in EPICA, as in SHEEC 1000-1899, are assessed from both i) MDPs, with homogeneous and repeatable procedures and ii) regional catalogues, coherently with those from MDPs (Stucchi et al., 2013). The two sets of parameters are then combined in the final parameters, and are both listed in the catalogue file together with the final parameters.

#### 1.1.3.1 Parameters from MDPs

Location with uncertainty, epicentral intensity, and magnitude from MDPs are jointly determined for 3297 earthquakes out of the 3622 with MDPs, because the intensity distribution of the remaining 325 events are not reliable enough to obtain robust parameters.

Macroseismic parameters are assessed with the same three methods as in SHEEC 1000-1899, i.e Boxer (Gasperini et al., 1999; 2010), Meep (Musson and Jimenéz, 2008) and BW (Bakun and Wentworth, 1997). The intensity attenuation models of each of the three methods are calibrated for five attenuation regions and are the same as in SHEEC 1000-1899 (see Gomez Capera et al., 2015) except for the Apennine-Dinarides (APD) region, substituted with the Boxer's calibration updated for CPTI15 (Rovida et al., 2016; 2020a). In addition, the new IO to Mw conversion strategy of CPTI15 is adopted in the newly defined Italian Volcanic Areas (IVA) region.

The choice of the method is mostly in favour of Boxer, with Meep and BW used as exceptions respectively for all the events in the UK and for a few offshore events. In each calibration region, a linear conversion relation from epicentral intensity IO to Mw is derived from the same datasets used for calibrating the three MDPs methods (Gomez Capera et al., 2015; Rovida et al., 2020a).

#### 1.1.3.2 Parameters from regional catalogues

Out of the 39 selected catalogues, 38 provide 5251 earthquakes in EPICA with location and epicentral intensity or magnitude, or both. Locations from the catalogues are adopted without any modifications, whereas magnitude is determined as in SHEEC 1000-1899, i.e.:
* Mw originally provided by catalogues is adopted without any modifications for 3316 earthquakes;

* otherwise, Mw is re-assessed from epicentral intensity I0 according to the regional conversion relations derived from the datasets used for calibrating the MDP methods (1831 events);

* for 56 cases neither Mw nor I0 is available, and Mw is converted from Ms or ML according to Bungum et al. (2003) and Grünthal et al. (2009); for additional 58 earthquakes the unspecified magnitude type from 4 catalogues is assumed as equivalent to Mw.

#### 1.1.3.3 Final parameters

Considering the sets of parameters described above, the final parameters in EPICA are determined as follows.

The epicentral location is selected according to the following criteria:

* When only one location from either MDPs or the selected regional catalogue is available, it is adopted

* When both locations from MDPs and from a regional catalogue, priority is given to the epicentre from macroseismic data.

In EPICA, 3297 (57.8%) epicentres are from MDPs, and 2257 (39.6%) from the selected regional catalogue. In addition, 149 epicentres (2.6%), marked as "preliminary", relate to earthquakes for which the available data do not allow a robust determination of the location. Among the locations from MDPs, 3187 are from Boxer, 82 from Meep, and 28 from BW (see Fig. 1.1). Moment magnitude is determined as follows:

* When Mw determinations from MDPS and from a regional catalogue are both available, they are combined through the weighted mean of the two values, with arbitrary weights of 0.75 and 0.25 attributed to the Mw from MPDs and from the regional catalogue, respectively. In continuity with SHEEC 1000-1899, reverse weights were given to ECOS-09 and CPTI15.

* When only the Mw from MDPs is available, it is adopted

* When only the Mw from a catalogue is available, it is adopted.

Mw is determined as the weighted mean for 3127 earthquakes (55% of the total), it derives from a regional catalogue for 2124 (37%) earthquakes, and from MDPs for 170 (3%) of them (Fig. 1.1). In addition, Mw is not determined for 282 earthquakes because the data they rely upon are not robust enough.

<figure><img src="image_2.png" /><figcaption>Fig. 1.1 Origin of the epicentral coordinates (left) and of the magnitude (right) in EPICA.</figcaption></figure>
1.1.4 EPICA, the European PreInstrumental earthquake CAtalogue

EPICA, the European PreInstrumental earthquake CAtalogue version 1.1 exploits the best knowledge of
the European historical seismicity today available, accounting for the most informative datasets among
those in the European Archive of Historical Earthquake Data AHEAD, and maximizes the harmonization
of earthquake data and parameters at country borders by means of uniform procedures.

EPICA contains 5703 earthquakes with either maximum observed intensity ≥5 or Mw ≥4.0, for the period 1000-1899. As a whole, the 64% of the earthquakes in EPICA is supported by MDPs, and the 97% is reported in a regional catalogue. In addition, for the 60% of the earthquakes both MDPs and a parametric catalogue are available, while the 36% is known from parametric catalogues only (Fig. 1.2). Conversely, the 3% of the earthquakes derive from an historical seismological study and is not included in any published catalogue.

Compared to SHEEC 1000-1899 EPICA contains 1035 new earthquakes, and the number of earthquakes supported by MDPs increased from 2447 to 3622 in EPICA, with 49852 considered MDPs instead of 42581. More than half of the added earthquakes relies on newly published sources of data. In addition, 339 earthquakes from studies or catalogues already considered in SHEEC 1000-1899 are added to EPICA because of the lowered intensity/magnitude threshold or of a revision of the dataset.

Half of the earthquakes in common between EPICA and SHEEC 1000-1899 (i.e. 2332 records) derive from
the same sources of data, and thus present the same parameters, whereas 2336 earthquakes listed in
SHEEC 1000-1899 are included in EPICA with a new or different set of data. In addition, 49 earthquakes
in SHEEC 1000-1899 are not included in EPICA because they have been recognized as fakes, or resulted
to be duplications.

The increased number of earthquakes with MDPs indicates an overall improvement of the knowledge
of European historical seismicity, being the MDPs distributions the results of thorough historical
investigation (see Rovida et al., 2020b).

The wealth of new input data results in an improved reliability and harmonization of earthquake
parameters. However, the availability and nature of historical earthquake data across Europe is still
inhomogeneous, with a large variability from one area to the other and substantial gaps mostly in
Northern and Central-Eastern Europe (Rovida et al. 2020b), which of course cannot but reflect on the
uniformity of any Europe-wide catalogue.
<figure><img src="image_3.png" /><figcaption>**Fig. 1.2** Earthquakes in EPICA by magnitude value and type, which reflects the typology and reliability of supporting data.</figcaption></figure>

## 1.2 Instrumental Earthquake Catalogue

The “instrumental” earthquake catalogue covers the period 1900 to the end of 2014 and builds on the previous European-Mediterranean Earthquake Catalogue (EMEC) from Grünthal & Wahlström (2012). Although it is referred to as an “instrumental catalogue”, the EMEC catalogue and its current update are composite catalogues constructed from local seismicity bulletins, harmonized earthquake catalogues and special studies across various regions of Europe.

Though the majority of events are instrumentally recorded, for much of the early 20th century the catalogues may contain earthquakes that are parameterized from macroseismic data or may contain a mix of macroseismically and/or instrumentally calibrated events depending on the region or data source in question. For many of the data sources, the earthquake magnitudes are calibrated in a mix of scales, including local magnitude ($M_L$) calibrated to different regional conditions, body-wave and surface-wave magnitudes ($m_b$ and $M_s$), moment magnitude ($M_W$) and in some cases maximum intensity ($I_{MAX}$) or epicentral intensity ($I_o$). To use the catalogue in PSHA, in addition to a range of other applications, these must be harmonized into a common scale directly taken from, or equivalent to, moment magnitude $M_W$.

The compilation procedure to update the EMEC catalogue follows closely the approach Grünthal & Wahlström (2012), which was modified slightly for application to the ESHM13 by Grünthal et al. (2013). Both catalogues begin in 1000 CE, but for the ESHM20 only the period from 1900 CE to 2014 CE is considered.

The update incorporates recently published data sets and special studies where available within certain regions and aims to extend the time-period covered from 1000 CE to the end of 2014 CE. To ensure a degree of continuity with respect to the earthquake catalogues adopted by the ESHM13, the same
hierarchy of data sources and the magnitude conversion formulae adopted are retained from the
original work of Grünthal & Wahlström (2012), except where new data sets necessitate a revision to this
process.

1.2.1 Data collection and regionalisation

The construction of the EMEC catalogue is based around a geographical/political regionalization of Europe, reflecting the extent of either political territories of each country or, in some cases, regions of influence that may be treated as a single region in the catalogue compilation.

The regionalization is shown in **Fig. 1.3**. The data collection process identifies new bulletins, compiled catalogue and special studies. The latter may consist of investigations into seismic sequences for which moment tensors were determined and corresponding moment magnitudes defined.

Events from each bulletin (in the form of a date, time, location, magnitude and data source) are added to a master database. Where new representations of existing events in the database are found, these are associated to the existing events based on simple spatio-temporal proximity criteria.

Among the new data sources integrated into the updated EMEC database are recently updated harmonized compilations in Italy (CPTI15, Rovida et al. 2020), France (F-CAT, Manchuel et al. 2018), Turkey (Kadirioglu et al. 2018), Slovenia (Zivcic, 2018) and Romania (INFP, 2018). Bulletin data comprises recent events reported by various seismic networks across Europe. A complete region-by-region description of the data sources and hierarchies is prepared in a separate report.

<figure><img src="image_4.png" /><figcaption>Fig. 1.3 EMEC regions for data compilation</figcaption></figure>
1.2.2 Harmonisation

Data are harmonized using a hierarchical strategy, prioritizing existing M<sub>w</sub> harmonized catalogues and
special studies, then considering moment tensor databases and subsequently local bulletins (without
M<sub>w</sub>). Finally, where no local bulletin can be identified, earthquake information is retrieved from the
International Seismological Centre (ISC).

For local bulletins and ISC events, the magnitude conversions adopted by Grünthal & Wahlström (2012) are retained, with the exception of a small number of cases where updated conversion formulae are identified in the literature. It is important to emphasise that in this approach the local sources of information are treated as the most authoritative sources, particularly when considering existing harmonized bulletins.

This aims to ensure some level of consistency with local scale catalogues that have been used in national seismic hazard modelling projects within Europe, though we note that differences nonetheless inevitably emerge as magnitudes and locations may change for events on the borders of respective countries.

The updated European-Mediterranean earthquake catalogue is shown in **Fig. 1.4**, and contains 55,732 events with M<sub>w</sub> (or equivalent proxy) ≥ 3.5 in the period 1900 to the end of 2014. The density of events per magnitude bin and year is shown in **Fig. 1.5**, alongside the equivalent plot from the Grünthal & Wahlström (2012) catalogue and the difference, i.e. the increase (or decrease) in number of events per cell. The majority of new data comes from the post-2006 period and from dropping the minimum magnitude in southern Europe from M<sub>w</sub> 4.0 to 3.5.

<figure><img src="image_5.png" /><figcaption>**Fig. 1.4** The complete EMEC catalogue for the period 1900 – 2012.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 1.5** Density of events per year and 0.1 M<sub>w</sub> bin width for the original EMEC catalogue of Grünthal & Wahlström (2012) (left), the current version (centre) and the change in number of events per cell (right).</figcaption></figure>

1.3 Active Crustal Faults and Subduction Zones

This section summarizes the main elements of the European Fault-Source Model 2020 (EFSM20),
developed as an update of the European Database of Seismogenic Faults 2013 (EDSF13; Basili et al.,
2013). For this update, we considered primarily the compilations that covered with a consistent
approach significantly large regions, relying on the work of the authors of each compilation for what
regards the accuracy and recentness of the information. We resorted to working on individual studies
or original work only in case of undefined situations, e.g., area of overlap between two regional
compilations or cases where a significant update was available or where the fault information was not
covered by EDSF13 but deemed necessary.

In EFSM20, two main categories of seismogenic faults are considered: 1) crustal faults; and 2) subduction zones. The collated datasets and the main curation aspects are summarized in the next two Sections. The map of collated datasets is shown in **Fig. 1.6**. The model covers an area that encompasses a buffer of 300 km around all target European countries (except for Overseas Countries and Territories, OTCs), and a maximum of 300 km depth for slabs.

1.3.1 Crustal faults

The minimum set of basic fault parameters required for constructing a seismogenic source model refer to Geometry (Location: Lat, Lon, Depth; Size: Length, Width; Orientation: Strike, Dip) and Behavior (Rake and Slip Rate). These are indispensable elements for devising and applying a fault recurrence model to be expressed by a Magnitude-Frequency Distribution (MFD).

Not all fault compilations provide this characterization in full, and strategies needed to be devised to fill
in the missing information. Concerning crustal faults, we identified several regional datasets that vary in
date of latest release, geographical extent, level of fault characterization, and more importantly, data
formats. In addition to those listed below, several others datasets were considered, such as for example
those covering Romania and the northern Black Sea (Diaconescu, Craiu, Moldovan, et al., 2019;
Diaconescu, Craiu, Toma-Danila, et al., 2019; Diaconescu et al., 2021). Although these datasets
represented a significant advancement with respect to EDSF13, they could not be used because did not
fully compliant with the requirements above recalled.

The list below summarizes the main adopted datasets in the various regions (Fig. 1.6).

• Dataset #01

This is the original database EDSF13, compiled in the framework of the SHARE project (Basili et al., 2013; Giardini et al., 2013) which covered Europe and the Mediterranean region. This was the starting dataset
to build the new fault model. The largest regions that remained unmodified are in the Balkans and northern Africa. Most regions were either entirely replaced by new datasets or partly revisited according to individual studies. Additions in regions that were not covered by EDSF13 are in Iceland, France, and the northern Mid-Atlantic plate boundary. The major regional updates are summarized below.

<figure><img src="image_7.png" /><figcaption>**Fig. 1.6** Map of collated fault datasets for the development of the European Fault-Source Model 2020 (EFSM20). See text for the descriptions of the various datasets. From west to east, the subduction systems are: Gibraltar Arc (GiA); Calabrian Arc (CaA); Hellenic Arc (HeA); and Cyprus Arc (CyA).</figcaption></figure>

* Dataset #02

This dataset covers the Mid-Atlantic ridge and transforms. The initial geometry was derived from a global plate-boundary model (Bird, 2003) and the rest of the characterization was based on the oceanic crust age and spreading rate (Müller et al., 2008). For the transform faults the slip rate is directly derived from the spreading rate, aided by more local data for the Gloria fault (Fadil et al., 2006; Koulali et al., 2011). For the normal faults the slip rate is derived by combining the spreading rate with local information about fault spacing and heave (Escartín et al., 1999; MacDonald & Luyendyk, 1977).

* Dataset #03

This dataset, which covers the French region, is derived from BDFA (Jomard et al., 2017). Due to the different strategy of fault mapping used in BDFA, we redrew the fault traces by interpolation and reassigned several parameters based on recent works in the region.

* Dataset #04

This dataset includes a few faults in the Gulf of Corinth. The initial geometry of the faults is based on GREDASS (Caputo & Pavlides, 2013) and the slip rate were updated based on more recent works (Bell et al., 2009; Fernández-Blanco et al., 2019).
• Dataset #05

This dataset covers the offshore regions of the Gulf of Cadiz and the Alboran Sea. In this area we updated the EDSF based on several recent works, either providing updated geometries and slip rates (Koulali et al., 2011; Martínez-Loriente et al., 2018; Martínez-Loriente et al., 2013; Neres et al., 2016).

• Dataset #06

This dataset covers the Italian territory and some surrounding regions. It is mainly based on the most recent version of the DISS (DISS Working Group, 2021).

• Dataset #07

This dataset includes a few faults in the Eastern Betics region. Such faults were modified with respect to the QAFI according to recent works with substantial updates of fault geometries and slip rates (Borque et al., 2019; Gómez-Novell, Chartier, et al., 2020; Gómez-Novell, García-Mayordomo, et al., 2020; Herrero-Barbero et al., 2020).

• Dataset #08

This dataset covers the Aegean region. It is mainly based on the most recent version of GreDaSS (Caputo & Pavlides, 2013).

• Dataset #09

In Iceland we started from the same approach as the Mid-Atlantic ridge and transforms (Dataset #02) and added data and considerations based on local studies (Árnadóttir et al., 2008; Bergerat et al., 1990; Forslund & Gudmundsson, 1991; Garcia et al., 2002; LaFemina et al., 2005; Rögnvaldsson et al., 1998; Rust & Whitworth, 2019).

• Dataset #10

This dataset deals with the Lower Rhine Graben. In this area we started from EDSF13, already based on a local fault model (Vanneste et al., 2013), and updated the slip rate of several faults based on more recent data (Gold et al., 2017).

• Dataset #11

This dataset covers the northwestern African region (Morocco, Algeria, and Tunisia). In this region, we mainly relied on EDSF13 with updates of a few faults in the Moroccan region based on the GEM Global Active Faults Database (Styron & Pagani, 2020) and various other works (Akoglu et al., 2006; Gomez et al., 1996; Pastor et al., 2015; Rigby, 2008; van der Woerd et al., 2014) for refining several fault parameters.

• Dataset #12

The NOAFAULTS database (Ganas, 2021) was used to integrate the dataset in the Aegean region for faults not already included in GreDaSS.

• Dataset #13

This dataset covers the Lower Tagus Valley Fault Zone (LTVFZ), southern Portugal, based on recent works (Canora et al., 2015) and Algarve (Sanz de Galdeano et al., 2020).

• Dataset #14
This dataset covers most of the Iberian region including the Pyrenees. In this region we relied on the Quaternary Faults Database of Iberia (QAFI) database (García-Mayordomo et al., 2012; IGME, 2015). Due to the different strategy of fault mapping used in QAFI, we redrew the fault traces by interpolation.

• Dataset #15

This dataset deals with Slovenia and its surroundings. In this area we relied on the recently published Database of active faults in Slovenia (Atanackov et al., 2021).

• Dataset #16

This dataset covers Anatolia and parts of the Middle East. In this region, we relied on recent data from the project EMME and data from the national update of the Turkish hazard model (Danciu et al., 2018; Demircioğlu et al., 2018; Emre et al., 2018).

To maximize the use of available fault compilations, crustal faults are represented with a down-dip planar geometry. In several compilations, however, crustal faults are represented only by the trace of the fault upper edge. To homogenize the collated dataset, the fault upper edge was resampled at a regular spacing of ~5 km. Whenever a fault was represented by multiple strands, the upper edge was reduced to a single trace by interpolation and smoothing. The fault plane was then extruded downward based on the dip and depth values (Fig. 1.7).

<figure><img src="image_8.png" /><figcaption>Fig. 1.7 Schematic showing the geometric construction of the fault-source model.</figcaption></figure>

For all datasets, we verified the depth datum (local ground surface or mean geoid/spheroid), then checked the possible intersections of the fault plane with the topo-bathymetry using the ETOPO1 Global Relief Model (Amante, 2009; NOAA National Geophysical Data Center, 2009) and with the base of the crust using the European Moho (Grad et al., 2009). Also, the possible intersections at depth between extruded fault planes were verified and, occasionally, this circumstance led to removal of a fault partly or even entirely. The overlaps between different datasets were harmonized based on a prioritization scheme which included considerations on newer data and national data, accuracy and justification of the fault information, coherence with the fault characteristics in the surrounding regions.
Once the fault reconstruction is completed, several geometric parameters, such as total length, end-to-end length, down-dip width, and area are calculated and stored in the GIS file attribute table. For what concern slip rates, we paid attention to harmonize the associated uncertainty. Occasionally, the original datasets we relied upon provided single slip rate values without uncertainty, or the difference between the minimum and maximum value was very small compared with what is usually found in the literature, or did not specified minimum values larger than zero. We thus performed a statistical analysis of the relative uncertainty in all records, considering the different categories of fault type (dip-slip and strike-slip) and tectonic domain (interplate and stable continental regions).

We then adopted the average relative uncertainty (percentage of slip rate), weighted for the number of occurrences in each category, to correct the slip rate variability where it was missing or not reliable. This procedure prevents having apparent small uncertainties for specific faults which could have depended solely on the different strategy of data collection in different datasets or regions.

The maximum earthquake magnitude of the crustal faults is estimated as the magnitude value, in the moment magnitude scale, that corresponds to the largest possible rupture that a fault can host based on its dimensions and magnitude scaling relations (Leonard, 2010, 2014).

The adopted scaling relations are in the form Mw=a+b*Log(S) where S can be the size of any of the following rupture dimensions: end-to-end length, width, area, or displacement, and the parameters "a" and "b" take different values depending on the S, the sense of slip (dip-slip or strike-slip), and the seismotectonic context (interplate or stable continent). The characterization of the faults for adopting the different scaling relations was established by using the provided rake values and following the definition and geographic distribution of the stable continents (Johnston, 1994).

A characterization of the level of activity of each fault is provided by the tectonic moment rate, that can be obtained by the product of fault area multiplied by slip rate multiplied by an assumed rigidity (or shear modulus) of the upper crust.

The parameterization of the crustal faults in EFSM20 is visualized in the maps shown in **Fig. 1.8**. The distribution of occurrences of the different parameters are displayed in the histograms of **Fig. 1.9**.
<figure><img src="image_9.png" /><figcaption>**Fig. 1.8** Maps showing different characteristics of the harmonized fault-source model. From top to bottom and left to right, the four panels show the faulting types (normal, reverse, transcurrent), the maximum magnitude, the slip rates and the moment rates. Slip rates and moment rates are color-coded based on a classification in the logarithmic scale.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 1.9** Histograms showing the frequency of occurrence of the different characteristics of the harmonized fault-source model. From top to bottom and left to right, the four panels show the faulting types (normal, reverse, transcurrent), the maximum magnitude, the slip rates and the moment rates. Slip rates and moment rates are color-coded based on a classification in the logarithmic scale. See correspondence with the geographic distributions in Figure 3.</figcaption></figure>
### 1.3.2 Subduction Zones

Subduction zones are complex systems where different types of potential earthquake sources exist. We rely on the assumption that the earthquakes occurring in the upper plate of a subduction zones are dealt with by the crustal fault model. We thus designed the subduction model to prospectively assess the slab interface seismicity and the intraslab seismicity. To this end the main necessary element is the tri-dimensional geometric reconstruction of the lower plate top surface. In addition, the minimum set of basic parameters required for constructing a seismogenic source model are the crustal thickness of the subducting plate, the upper and lower depths of the seismic interface, and the net convergence direction and rate.

Concerning subduction sources, we identified several datasets that vary in date of the latest release, level of characterization, and data formats. We started with EDSF13, the original database compiled in the framework of the SHARE project (Basili et al., 2013; Giardini et al., 2013) which covered the Mediterranean subduction systems, which include the Calabrian Arc, the Hellenic Arc, and the Cyprus Arc (Figure 1). This was the starting dataset to build the new subduction zone model. In addition to searching the literature for more updated data, we conducted an original work to reconstruct the subduction system of the Gibraltar Arc (Fig. 1.6) which was not included in EDSF13.
The geometry of the Calabrian Arc slab was entirely updated using a large dataset of seismic reflection profiles and the seismicity distribution (Maesano et al., 2017). The geometry of the Hellenic Arc and Cyprus Arc slabs were recently revisited in the framework of a tsunami hazard project (Basili et al., 2021). These datasets have also been re-examined in light of the SLAB 2 model (Hayes et al., 2018) and several others (Ganas & Parsons, 2009; Halpaap et al., 2018, 2019; Sachpazi et al., 2016).

The geometry of the Gibraltar Arc slab was reconstructed using different datasets. For the shallowest depths, we used data of bedrock markers based on the interpretation of multichannel seismic reflection profiles and wide-angle seismic surveys (Gutscher et al., 2009) assuming that the top of the slab coincides with the top of the basement. For the intermediate depths (12-40 km), we used a model of the Moho obtained from a set of diverse datasets using a probabilistic surface reconstruction algorithm (Arroucau et al., 2021) and considering typical values for the old Tethys oceanic crust in the range of 7-9 km (Sallarès et al., 2011). Then, between 40-70 km depth, we obtained the position of the slab by interpolation of seismicity clusters. Within the 140-200 km depth range, the slab was assumed to be vertical based on a tomographic model (Civiero et al., 2018), which shows a nearly vertical high-velocity p-wave anomaly down to 600 km deep.

All these 3D geometries were then resampled at regular spacing and smoothed to ensure the same spatial resolution in the different models. Considering the crustal thickness of the lower plate, the geometry reconstruction is then completed with a lattice of 10x10 km laying on the slab intermediate surface. The various steps of the slab geometric construction are shown in Fig. 1.10, while the resulting geometries of the four slabs are shown by oblique views in Fig. 1.11.

The upper and lower depths of the seismic interface was estimated from data and modeling of the 150°C and 350-450°C isotherm, the seismicity distribution, and the slab intersection with the Moho of the upper plate (Davies, 2013; Grad et al., 2009; Gutscher et al., 2006; Heuret et al., 2011; Syracuse et al., 2010; Thiebot & Gutscher, 2006).

Tectonic rates are derived from geodetic observations on the plate convergence across the subduction zones (Carafa et al., 2018; Devoti et al., 2008; Hollenstein et al., 2008; Howell et al., 2017; Nocquet, 2012; Palano et al., 2015; Reilinger et al., 2006; Stich et al., 2006; Wdowinski et al., 2006).
The maximum earthquake magnitude of the slab interface is estimated as the magnitude value, in the moment magnitude scale, that corresponds to the largest possible rupture that the interface can host
based on its dimensions and magnitude scaling relations specifically developed for this type of earthquake ruptures (Allen & Hayes, 2017).

For what concern the tectonic moment rate, in the case of the slab interface an important aspect to consider is the rigidity (shear modulus) variation with depth that has been observed in subduction zones from around the world (Bilek & Lay, 1999; Sallarès & Ranero, 2019) and already used to simulate earthquake ruptures for tsunami simulations (Geist & Bilek, 2001; Scala et al., 2020).

<figure><img src="image_11.png" /><figcaption>**Fig. 1.10** Schematic showing the geometric construction of the slab model.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 1.11** Oblique views of the tri-dimensional geometry of the four slab models.</figcaption></figure>

**1.4 Tectonic Regionalization**

A tectonic regionalization is a subdivision of the entire source space relevant for the Euro-Mediterranean area into discrete regions that are internally as homogeneous as possible from the standpoint of the dominant tectonic processes.

This regionalization was initially built following basic principles of plate tectonics and building on previous experience of the ESHM13 (Delavaud et al., 2012) and the TSUMAPS-NEAM project (Basili et al., 2021). The following eleven tectonic settings are defined for the crustal level:

1. Active volcano;

2. Back-arc and orogenic collapse;

3. Continental rift;

4. Oceanic rift;

5. Contractional wedge;

6. Accretionary wedge;

7. Conservative plate margin;

8. Transform zones proper;

9. Shield;

10. Stable continental region;

11. Stable oceanic region.

The geographic distribution of the above-defined tectonic settings forms the regionalization shown in **Fig. 1.12**. This construct is a 2D subdivision of the crustal volume. Most of the boundaries between adjacent zones can be assumed to be nearly vertical. However, for some of them this assumption does not hold and further considerations are needed. This is especially the case of the subduction zones,
where the presence of slabs implies a 3D geometry of the regions classified as accretionary wedges. Inclined boundaries exist also elsewhere. In addition, in some cases it should be considered that the transition between one region and another could be gradual and that some regions could have an intermediate classification between those envisaged beforehand.

An important element used to constrain the first draft of the regionalization is represented by the geometry of the main plate boundaries (Bird, 2003). Subdivisions of the plate interiors was mainly based on large-scale geologic maps and tectonic classifications (Asch, 2005; Johnston, 1994; Müller et al., 2008). For the volcanoes, we relied on the global database without modifications (Global Volcanism Program, 2013).

<figure><img src="image_13.png" /><figcaption>**Fig. 1.12** Map of the tectonic regionalization. Each region is color-coded depending on the tectonic setting.</figcaption></figure>
## 2 Unified Earthquake Catalogue

The unified earthquake catalogue consists of two parts: the so-called “instrumental” catalogue (after 1900) based on the updated EMEC catalogue (Weatherill et al. 2021), and the European PreInstrumental Earthquake Catalogue EPICA v. 1.1 (Rovida and Antonucci, 2021; https://doi.org/10.13127/epica.1.1) earthquake catalogue (between the years 1000 CE and 1899). Both catalogues are documented in Chapter 1.

Table 1 summarizes the unified ESHM20 catalogue and its two main components. For comparison, we also include the summary of the unified ESHM13 catalogue. ESHM20 unified earthquake catalogue includes an additional period of 8 years since the end of 2006 and a larger spatial coverage. Thus, it naturally contains more events. Even for the same period (1000-2007), the ESHM20 catalog is overall richer in the number of events at different magnitude thresholds (Fig. 2.1).

However, the ESHM20 catalogue is not homogeneously richer everywhere in space for a given magnitude threshold. For a given magnitude threshold, the difference in the number of events, Δ, between the two catalogs varies widely in space, with some regions featuring more events and others less. We also find that Δ varies substantially with the magnitude threshold.

Especially striking is the Italian region, where the ESHM20 catalogue features substantially more events at a magnitude threshold of 3.5, while at the magnitude threshold of 4.5, the ESHM13 catalog is comparatively richer as a consequence of the updated Italian catalogue CPTI15 (Rovida et al., 2020).

<table><caption>**Table 1** Summary of the updated and unified earthquake catalogue used in ESHM20 and the earthquake catalogue for ESHM13 (SHEEC v3.3).</caption><thead><tr><th></th><th>(1) Historical: EPICA v1.1</th><th>(2) Instrumental: EMEC v20190218</th><th>(1)+(2): Unified ESHM20</th><th>Compare Unified ESHM13</th></tr></thead><tbody><tr><td>Time span</td><td>1000/1899</td><td>1900/31-12-2014</td><td>1000/31-12-2014</td><td>1000/31-12-2006</td></tr><tr><td>Magnitude span</td><td>1.9/8.5</td><td>3.49/8.3</td><td>1.7/8.5</td><td>1.7/8.5</td></tr><tr><td>Longitude span</td><td>-23.5°/32.4°</td><td>-37.0°/51.9°</td><td>-37.0°/51.9°</td><td>-31.65°/45.0°</td></tr><tr><td>Latitude span</td><td>35.0°/69.4°</td><td>26.9°/73.0°</td><td>26.9°/73.0°</td><td>33.2°/73.32°</td></tr><tr><td># Total Number</td><td>5703</td><td>55411</td><td>61127</td><td>30012</td></tr><tr><td># magnitude ≥ 4.5</td><td>2337</td><td>20388</td><td>22728</td><td>13284</td></tr><tr><td># magnitude ≥ 5.0</td><td>1550</td><td>6013</td><td>7565</td><td>5585</td></tr><tr><td># magnitude ≥ 5.5</td><td>885</td><td>1920</td><td>2805</td><td>2066</td></tr><tr><td># NaN depth</td><td>4965</td><td>9451</td><td>14422</td><td>10616</td></tr><tr><td># NaN magnitudes</td><td>282</td><td>0</td><td>282</td><td>303</td></tr><tr><td># depth < 40km</td><td>5355</td><td>49436</td><td>54791</td><td>25666</td></tr></tbody></table>
SHEECv3.3 - SERA; Time Span = 1900-2007

<figure><img src="image_14.png" /><figcaption>**Fig. 2.1** Spatial variation of the difference in number of events, $\Delta$, between the ESHM20 and ESHM13 catalogs for different magnitude thresholds. Warmer and cooler colors indicate regions where ESHM20 and ESHM13 catalog feature more events, respectively.</figcaption></figure>

## Unified Earthquake Catalogue: 2020 update for Iceland

Jónasson et al. (2021) provide an update of the instrumental earthquake catalogue for Iceland. In this updated version the authors have revised the Mw conversion acknowledging that the equations used in the 2013 earthquake catalogue (Grünthal and Wahlström, 2012; Grünthal et al., 2013) resulted in conservative Mw estimates for early events. The differences of Mw between the two catalogues is time-dependent as reported in Jónasson et al (2021), i.e. the average difference is 0.41 before 1970, 0.37 between 1970 and 1980, and 0.27 after 1980; third quartiles are 0.59, 0.47, and 0.36, respectively and the locations are reasonably accurate. Following the recommendation of the Icelandic colleagues, the core team used the earthquake catalogue by Jónasson et al. (2021) for ESHM20.

## 2.1 Catalogue Declustering

One of the assumptions of the time independent probabilistic seismic hazard analysis is that earthquake occurrence times can be represented by a stationary Poisson process. However, catalogues contain dependent events, such as foreshocks aftershocks or triggered earthquakes, and exhibit strong spatial and temporal clustering of seismicity.

Thus, the catalogues of recorded earthquakes usually reveal noticeable departures from the Poisson distribution. To obtain catalogs that follow a stationary Poisson process, the reported catalogs need to be declustered.

Declustering algorithms categorize the earthquakes in a catalogue into groups of independent and triggered occurrences and fulfil two main purposes:

- Creating a catalogue of expected mainshocks using a stationary Poisson process.

- Removing the spatial bias in the seismicity rate across the entire catalogue produced by seismicity's intensive spatio-temporal clustering.

There are numerous declustering techniques available, which reflect the complexity of assigning such labels. Often, the methods are used without considering their effect on the final results.

ESHM13 employed a straightforward windowing scheme based on Grünthal (1985) window sizes. To account for foreshock activity, the original method was slightly modified: if a greater magnitude event is discovered inside the time- and space-window of a detected mainshock, the algorithm retains the foreshock cluster, extends it to the window size of the larger event, and merges all events into one cluster, leaving just the largest event in the declustered catalog. This adjustment is reflected in the implementation of the hazard modeler toolkit (https://github.com/gem/oq-engine; Weatherill et al 2012) by the "fs time prop" input parameter (see table below).
As a first step, we used the SHEECv3.3 catalog from ESHM13 to derive this input parameter and thereby define a reference method. We discovered that setting "fs time prop" to 0.17 results in a declustered catalog that matches the published mainshock catalog for ESHM13. This method will be referred to as our reference approach (winGT fs0.17). In the framework of ESHM20 we investigated the following three main declustering methods:

1. Windowing techniques with space-time windows as defined by Gruenthal (personal communication), (Gardner & Knopoff, 1974) and (Uhrhammer, 1986). The corresponding declustering methods are referred to as winGT, winGK and winUH, respectively. All those windowing techniques are based on the idea that for each earthquake in the catalogue, the following events are labelled as aftershocks if they occur within a specified time-span, and within a distance range. Foreshocks are similarly treated and if the largest earthquake occurs later in the sequence, the foreshock is labelled as the aftershocks. Thus, the temporal-spatial windows are reset according to the magnitude of the largest earthquake in a sequence.

2. Cluster method introduced by (Reasenberg, 1985). The method aims at identifying dependent earthquakes by linking them to clusters according to spatial and temporal interaction zones. Such clusters thus typically grow in size when processing more and more events. The spatial interaction zone is chosen according to the stress distribution near the mainshock, and the temporal extension of the interaction zone is based on Omori's law. A cluster is then composed of all linked events, and the largest earthquake is labelled as the mainshock whereas smaller earthquakes are classified as fore- and aftershocks. The procedure applied here does not account for the uncertainties on the earthquake location.

3. Declustering based on a correlation metric (Zaliapin et al., 2008). The technique is based on a space-time metric to correlate earthquakes with each other. (Zaliapin et al., 2008) showed that by comparing rescaled times to rescaled distances it is possible to identify two distinct populations. These two categories consist of the time-stationary Poissonian mainshock events and the dependent events characterized by much smaller time and space inter-event distances. The method does not provide binary labels of mainshock and aftershock, but gives, for each earthquake, the associated probabilities that it is an independent event.

### 2.1.1 Sensitivity Analysis

We explore the sensitivity of the number of declustered events to different choices of declustering methods and corresponding parameters. We first considered the crustal part of the unified earthquake catalogue for ESHM20 as input for seismicity declustering. We define crustal seismicity as the part of the catalog that features (1) depth values < 40km or (2) *NaN* depth values. This definition leads to a total of 54791 crustal events. We report a large variability in terms of number of mainshocks when comparing the different techniques. This variability is more dominant than the variability given by parameter changes for a single selected method (see Fig. 2.2). Thus, we selected three methods for further investigation with the default input parameters given in bold font in Table 2. The results of the sensitivity analysis are given in Table 3. The method that generates the largest number of events after declustering is the Reasenberg method with 43563 events identified as main events, followed by Zaliapin method with 30982 main events and default time windows technique used in ESHM with 22708 main events.

The results of the sensitivity analysis are given in Table 3. The method that generates the largest number of events after declustering is the Reasenberg method with 43563 events identified as main events, followed by Zaliapin method with 30982 main events and default time windows technique used in ESHM with 22708 main events. However, these results do not provide a straightforward answer in which
method to be chosen. Thus, at a glance one could use the default time-window declustering technique, which was used in ESHM13. This method is calibrated for the seismicity in central Europe, and it might provide an aggressive declustering as the number of events remaining is the lower among the considered declustering techniques. A spatial comparison of the total number of events identified as mains by the three declustering methods is given in Fig. 2.3. The comparison is done in grid cells of 50km and only events of magnitude greater than 4.5 are illustrated. It can be seen that when compared with the WT-Grünthal the Reasenberg method identifies more main events in stable continental regions (Central Europe, Iberia, whereas in the southern Europe this trend is reversed. On the same figure (bottom plot), the spatial comparison between the WT-Grünthal the Zaliapin method, suggests that the Zaliapin method might be suitable for very active regions (i.e. Southern Europe), as in the central Europe the number of main events identified are below the WT-Grünthal method. However, these are just some assumptions not supported by any other evidences, hence we progress with a data-driven procedure to aid the selection of the suitable declustering technique. The results of this data driven staticial analysis indicate that the Reasenberg technique is an alternative method to WT-Grünthal method. Given the computational demand of the model, it was decided not to include an alternative branch for the declustering algorithm.

<table><caption>**Table 2** Input parameters of the main declustering techniques under investigation for use in declustering the unified earthquake catalogue. The bolded parameters are considered as the reference input parameters.</caption><thead><tr><th>Parameter</th><th>Declustering Method</th><th>Description</th><th>Range/(ESHM20 default)</th></tr></thead><tbody><tr><td>fs_timeProp</td><td>Window</td><td>fraction of the time window used for foreshocks</td><td>[1.0, 0.5, 0.1, 0.17, 0.01, 0.001]</td></tr><tr><td>taumin</td><td rowspan="4">Reasenberg</td><td>Minimum look ahead time for clustered events</td><td>[1.0] day</td></tr><tr><td>taumax</td><td>maximum look ahead time for clustered events</td><td>[10.0, 50.0 100 300] days</td></tr><tr><td>xmeff</td><td>effective magnitude to define magnitude cutoff (with xk=0.5)</td><td>[3.5, 4.5, 5.5,6.5, 7.5]</td></tr><tr><td>rfact</td><td>factor for interaction radius for dependent events</td><td>[10.0, 20.0 30.0]</td></tr><tr><td>fractal_dim</td><td rowspan="3">Zaliapin</td><td>spatial weighting factor</td><td>[1.4, 1.6, 1.8]</td></tr><tr><td>b_value</td><td>magnitude weighting factor</td><td>[0.9, 1.0, 1.1]</td></tr><tr><td>theta</td><td>temporal weighting factor</td><td>[0.5, 1.0, 2.0]</td></tr></tbody></table>

<table><caption>**Table 3** Results of the sensitivity analysis on the declustering techniques.</caption><thead><tr><th>Unified catalogue</th><th>ESHM20-full</th><th>ESHM20-winGT</th><th>ESHM20-rbOrig</th><th>ESHM20-Zaliapin</th></tr></thead><tbody><tr><td># Total Number</td><td>54791</td><td>22708</td><td>43563</td><td>30982</td></tr><tr><td># magnitude ≥ 5.0</td><td>6726 (100%)</td><td>3972 (59%)</td><td>5583 (82%)</td><td>3641 (54%)</td></tr><tr><td>Instrumental part</td><td>ESHM20-full</td><td>ESHM20-winGT</td><td>ESHM20-rbOrig</td><td>ESHM20-Zaliapin</td></tr><tr><td># Total Number</td><td>49436</td><td>18784</td><td>38470</td><td>27526</td></tr><tr><td># magnitude ≥ 5.0</td><td>5245</td><td>2795</td><td>4106</td><td>2337</td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>**Fig. 2.2** Sensitivity of number of declustered events to the choice of declustering methods and corresponding parameters.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 2.3** Spatial difference of the number of events with magnitude M >4.5 in a 50km grid cell, comparing the time-window based method used in ESHM13 (winGT fs0.17) with rb=Reasenberg (1985) (left), and with Zaliapin=Zaliapin et al. (2008) (right).</figcaption></figure>

## 2.2 Catalogue Completeness

The magnitude of completeness (Mc) of earthquake catalogs is a significant metric that has critical consequences for probabilistic seismic hazard assessment, principally because the correct calculation of the parameters of the Gutenberg Richter (GR) equation is contingent on the reliability of Mc. The importance of correct estimation of Mc is illustrated in **Fig. 2.4**.
<figure><img src="image_17.png" /><figcaption>**Fig. 2.4** a) Spatial distribution and magnitude of all earthquakes reported in the ESHM20 earthquake catalog; (b) b-value as function of assumed magnitude of completeness ($M_c$); error bars indicate the 95% confidence interval of the estimate of b-value; (c) overall empirical magnitude distribution for the entire catalog; bestfit GR law obtained assuming $M_c = 4.5$ for the entire catalog; Maximum likelihood estimates of the a and b-value are indicated.</figcaption></figure>

Assuming a too small value of $M_c$ leads to an underestimation in the b-value. As the value of $M_c$ is increased, the bGR-value first increases and then saturates at a constant level. Increasing the value of $M_c$ beyond the saturation point, then only increases the error in the estimate of the b-value, as the number of earthquakes with magnitudes $\ge M_c$ decreases exponentially.

Using the stability in the estimates of b-values as a proxy for the completeness of the catalog (Cao & Gao, 2002), we can consider the shallow ESHM20 catalog to have a bulk? completeness for $M \ge 4.5$. Using $M_c = 4.5$, we find that the global estimates of (aGR, bGR) values for the ESHM20 catalog are (8.43, 0.91), respectively.

In PSHA, one frequently deals with earthquake catalogs very heterogeneous in space and time, with varying degrees of completeness and non-uniform histories. As a result, an average $M_c$ for the entire region, especially large ones, and period of the catalog may not be the most appropriate method, resulting in a skewed estimation of region-specific GR parameters.

This issue is frequently addressed in PSHA by first identifying completeness zones (SZs), which are assumed to have homogeneous completeness history. Then methods that essentially employ changes in earthquake reporting rates in different magnitude bands to identify time periods and corresponding magnitudes of completeness are applied.

To the best of our knowledge, the applications of these methods are semi-quantitative, in that an expert determines both the time-intervals and the corresponding magnitude of completeness by visually analyzing the slope changes in the graph of either the cumulative number of earthquakes with time or from "Stepp" plots (Nasir et al., 2013; Stepp, 1972).

We address the primary deficiency of these methods by making them completely quantitative, replacing expert judgment with appropriate statistical tests. We primarily focus on quantifying the temporal course of earthquake frequency "TCEF" method and combine it further with the Maximum curvature method to achieve the best of both methods.
### 2.2.1 Completeness zones (CSZs)

To simplify the problem of estimation of the spatio-temporal variation of $M_c$, as a first step, the study region is divided into completeness zones (CSZs).

These CSZs are delineated starting from the same regions as those of the EMEC catalogue harmonization (See Fig. 1.3), with additional criteria based on the homogeneity of earthquake reporting through time history, i.e., the earthquake reporting rates are thought to be spatially homogeneous within each of the CSZs.

**Fig. 2.5(a)** shows the polygons outlining the boundary of the 48 defined CSZs. The unique Id assigned to each of these zones is shown within the corresponding polygon in **Fig. 2.5(a)**. The color of the polygons shown in Figure 5a scales with the log of the number of earthquakes reported within each polygon. The number of earthquakes reported within the polygons varies approximately by a factor of 3, according to the underlying seismicity of the region and the duration and characteristics of the catalog.

As shown in **Fig. 2.5(b)**, the duration of the earthquake catalog also varies by a factor of 2, with minimum and maximum duration being 100 and 1000 years. The variation in the duration of the catalog already points to the necessity of this spatial zonation, so as to ensure that regions with different completeness histories are analyzed independently of each other.

**Fig. 2.5(b)** also indicates the range of magnitudes reported in each CSZs. While the upper magnitude range is both a factor of maximum possible magnitude within a region as well as the completeness history, the variation in the reported lower magnitude range points to the heterogeneity of the completeness history from one CSZ to another and to the lower threshold of the regional catalogue.
<figure><img src="image_18.png" /><figcaption>**Fig. 2.5** a) Boundaries of the completeness super-zones and their IDs; the color of each completeness zone scales with the log of the number of events reported within it in the catalog; (b) duration and magnitude range of the catalog reported within each completeness zone.</figcaption></figure>
### 2.2.2 Improved temporal course of earthquake frequency (TCEF) method

Because the temporal course of earthquake frequency, herein after TCEF approach, is primarily reliant on the premise that earthquake rates are Poissonian in time, its first step consists of declustering the data.

To do this, a declustering technique is employed to first preprocess the earthquake catalogs by eliminating clustered seismicity. In this work, we decluster the catalog before finding $M_c$ using the window-based declustering approach with (Gruenthal, 1985) scaling relations.

Following this procedure, the remaining declustered catalog is subdivided into magnitude intervals. A time series of the cumulative number of earthquakes is presented for each magnitude interval. If the catalog is adequately declustered, it should generally follow a Poissonian distribution; consequently, the time series of the total number of earthquakes should seem straight. However, changes in the reporting rate of earthquakes may cause changes in the slopes of these cumulative time series, which can be used as a proxy for changes in the catalog's completeness.

These change points are graphically identifiable in the time series of all magnitude intervals. The change points for distinct magnitude intervals are ordered using the assumption that $M_c$ diminishes with time as more and smaller earthquakes are recorded. This is done in a methodical manner, beginning with the greatest magnitude bin and progressing to the lowest magnitude bin. The constraint that the change point in a particular magnitude bin may only be greater than or equal to the change point in the previous magnitude bin is followed each time. The majority of these transition points are detected visually.

To avoid the qualitative aspect of this approach, we enhance it by including statistical checks to detect the change points. First, we begin by defining a multiple magnitude threshold that may be used to evenly split the reported magnitudes.

Secondly, all earthquakes with magnitudes greater than the highest magnitude criteria are selected. The timings of these earthquakes are then used to undertake a statistical test of uniformity. In this test, we calculate the difference between the observed empirical temporal distribution's cumulative density function (CDF) and the theoretical CDF predicted from a perfectly uniformly distributed sample.

Next, this difference is compared to the distribution of difference calculated for random data simulated from a uniform distribution. Finally, the Poissonian hypothesis may be rejected if the observed time series distance is much greater than the distances for the simulated time series.

This procedure is repeated by gradually extending the lengths of the periods chosen for the uniformity test. The moment when the test can no longer be rejected is regarded as the change point for that magnitude threshold.

The magnitude of completeness for this change-point and magnitude-threshold combination is assessed by choosing all earthquakes since this change point and over this magnitude threshold and then using the maximum curvature technique to detect $M_c$ for the chosen data.

This method is then repeated for each magnitude threshold, with the condition that the change point for a lower magnitude threshold may only be greater than or equal to the change point for a higher magnitude threshold. Furthermore, the $M_c$ corresponding to a lower magnitude threshold must be less than or equal to the $M_c$ corresponding to a higher magnitude threshold. A detailed step-by-step algorithm of the method is provided in the Nandan et al 2022 (in preparation).
Furthermore, we apply the Maximum-Curvature approach (Wiemer & Wyss, 2000) in the proposed algorithm to estimate the magnitude of completeness for each given completeness period. We chose Maximum-Curvature since it is non-parametric and rather easy to apply. However, more complex approaches for measuring the amount of completeness are available (Mignan & Woessner, 2012).

### 2.2.3 Time Series of Mc

Fig. 2.6 shows a time series of the magnitude of completeness obtained by applying the proposed approach to CSZ02 - Switzerland.

Panels (a) and (d) of each of these images depict the completeness zone and the data to which the suggested technique is applied. Panel (d) shows the location of the earthquakes, with the size of the dots according to the magnitude of the earthquakes. Panel (a) displays the magnitude-time plot of the data for the chosen SZs. Because the approach is only applied to the declustered catalog, a distinction is made between clustered (black) and background seismicity markers (orange). The completeness time series generated by applying the suggested technique to the zone is shown using a solid blue line. A dashed blue line is used to represent the completeness times for each zone proposed in previous studies or estimated for the ESHM13.

Note that these previous completeness estimates are based on different earthquake catalogs. The time intervals in Figure 6a, for example, were produced using the national catalogs utilized in the different national seismic hazard analyses for Switzerland (Wiemer et al 2015).

We discover a significant difference between the completeness time series produced in this investigation and those reported in previous studies. One probable explanation for this disparity is the difference in both the analyzed catalogs and the different method used to define the complete periods.

The time intervals of completion for all CSZs are presented in the Appendix A of this Report.

### 2.2.4 Objective quality assessment of time series of Mc

To determine whether the completeness from the various studies or the one determined using the algorithm proposed in this study is a better choice for the ESHM20 earthquake catalog, we measure the consistency between the theoretical and empirical annual frequency magnitude distribution, which can be derived by evaluating completeness time periods.

Using the duration of the complete period $T(m)$ of each magnitude bin, one can obtain the number of earthquakes in such a period in each magnitude bin $N(m)$, which can then be annualized by dividing by $T(m)$ to obtain $N_{annual}(m)$.

Finally, $N_{annual}(m)$ can be used to obtain the empirical annualized frequency magnitude distribution by converting $N_{annual}(m)$ to $N_{annual}(\ge m)$. The empirical annualized frequency magnitude distributions resulting from the two completeness assessments are shown in panel (c) and (d) of Figure 6 using orange circles. Using the data above completeness and the length of the complete periods, we can also obtain the parameters of the GR law, aGR and bGR values, which best describes the corresponding complete data. The maximum likelihood estimate of these two parameters is obtained for the two completeness assessments and are indicated at the top of panel (b) and (c) in Fig. 2.6. By comparing the panels (b) and (c) in Fig. 2.6, we find that the theoretical GR distribution obtained using the completeness assessment of this study seems (visually) more consistent with corresponding
empirical annualized frequency magnitude distributions than for the case of the regional completeness. We then quantify this goodness of fit in terms of the KS distance between the empirical and theoretical cumulative density functions.

Furthermore, **Fig. 2.7 (a)** shows the plot of KS distance obtained from the completeness time steps in this study versus those obtained from the regional completeness time periods for all the SZs. Out of 48 zones, for which both types of completeness estimates are available, in 38 cases the completeness times obtained in this study result in lower KS distance, indicating a higher degree of consistency between the empirical and theoretical cumulative density functions of magnitudes. The average KS distance corresponding to the regional completeness times is 0.24, while for the periods obtained in this study is only 0.14. Using the paired Student's t-test, we verify that the mean KS distances for the regional completeness are significantly larger than those obtained from the method proposed in this study. The null hypothesis that the two means are equal can be rejected at a significance level of 0.01 with a p-value of ~2.2× 10⁻⁶.

<figure><img src="image_19.png" /><figcaption>**Fig. 2.6** a) Comparison of the magnitude of completeness time steps obtained using the improved TCEF method (solid blue line) and expert-driven magnitude of completeness time steps (dashed blue line) for four SZs; Orange and black circles indicate the time and magnitude of the declustered and clustered earthquake events; Only declustered data is used to obtain the completeness time steps; (b) Consistency between the empirical annualized frequency magnitude distribution obtained using the expert-driven magnitude of completeness (orange circles) and best fit theoretical GR distribution; Maximum likelihood estimate of the a and b-value and their 95% confidence interval are indicated at the top of the panel; (c) Same as panel (b) but for the completeness time steps obtained using the improved TCEF method; (d) area of completeness super zone highlighted in blue; location of the declustered and clustered seismicity.</figcaption></figure>
In Fig. 2.7 (b), we also plot the *log N*<sub>comp</sub> or log of the number of complete events identified by the two types of completeness for each of the SZs. We find that the proposed method identifies fewer complete events (~1400) than the regional completeness. This tendency of being over-optimistic about completeness is also reflected in the estimates of b-value obtained (Figure 3b) being generally smaller for the case of regional completeness. Of the 48 SZs, the b-values estimated using the regional completeness are smaller for the 34 SZs. Using the Student's t-test, we again verify if the mean b-value for the regional completeness is significantly smaller than the mean bGR-value for the completeness obtained in this study. We can reject the null hypothesis that the two means are equal at a significance level of 0.01 with a p-value of 0.0021. The tendency of lower resulting b-value indicates that the regional completeness is likely more optimistic in its definition of complete events.

<figure><img src="image_20.png" /><figcaption>Fig. 2.7 (a) KS distance between the empirical annualized frequency magnitude distribution and the best fit theoretical GR law obtained using the expert driven completeness and completeness time steps inferred in this study. The red and blue markers indicate that KS distance for the expert driven completeness is larger (worse) than the completeness time steps inferred in this study and vice versa; (b) Number of events above completeness times steps for the two types of completeness time steps.</figcaption></figure>
# 3 Seismogenic Source Model

The seismogenic source model of ESHM20 is built upon the foundation of the ESHM13 (Wössner et al 2015), enhanced with newly developed local and regional seismogenic sources. The main seismogenic source model consists of four distinct source models aimed at capturing the spatial and temporal uncertainty of the earthquake rate forecast across the entire pan-European region:

- area sources model, which is assumed to be the pan-European consensus model, incorporating the national area sources provided by local experts and fully cross-border harmonization.

- active faults and background smoothed seismicity, a hybrid seismicity model that combines the updated active faults datasets (chapter 3, section 3.3.) with the background seismicity in regions where faults are identified. In regions without active faults – the smoothed seismicity model represents an alternative to the area sources model.

- subduction sources depicting both the subduction interface and in-slab of the Hellenic, Cyprian, Calabrian and Gibraltar Arcs.

- non-subducting deep seismicity sources describe the nested seismicity with depth in Vrancea, Romania as well as the cluster of deep seismicity in southern Iberia Peninsula.

In this Chapter an overview of these seismogenic source models, their components and properties, as well as their characterization in terms of earthquake rate forecast is provided.

## 3.1 Unified Area Sources Model

Area sources are the most common seismogenic source representation (e.g., Grünthal et al., 1999a, 1999b, 2018, Jimenez et al., 2003, Meletti et al. 2008; Vilanova et al. 2014, Danciu and Giardini 2015, Danciu et al 2018, Sesetyan et al 2018, Pagani et al 2018).

The development of the area source model is achieved by updating the area sources of ESHM13 with the contribution from the national hazard models. Since 2013, the release of the ESHM13, updates of the seismic hazard assessment at national level have been conducted across Europe.

The following countries have been providing updates of the seismogenic sources: Belgium, Bulgaria, Switzerland, Germany, Macedonia, Romania, Slovenia, Spain, Switzerland, Turkey; contribution for Northern Africa are taken from global mosaic of hazard models of Global Earthquake Model (Pagani et al 2018), the models for the Balkans are the same as in the BSHAP2 project (Salic et al 2016).

From every national model provided by the regional and national experts we have retained the area source model, or the branch model, that was recommended when more than one model was available (i.e. Model C for Germany, Grünthal et al 2018).

The area sources model describes regions of homogenous seismicity classified as: shallow crustal, volcanic, subduction in-slab and deep seismicity (i.e. Vrancea region, Romania). The delineation of the area sources was revisited and then modified to ensure cross-border harmonization. Mainly, the overlapping area sources at national borders were resolved while preserving the original information (i.e. depth distribution, style-of-faulting, maximum magnitude). As guidance for cross-border harmonization, we prioritize the national source models within the country boundaries.

The cross-border sources are guided by seismotectonic evidence, active faults, and major geologic/tectonic features if available, if not, then the seismicity patterns are used: historical earthquake locations or recent clusters of seismicity. If none is available, then area sources from the most recent seismogenic model are subjectively preferred. Furthermore, the curation of the area sources includes removing duplicates, simplifying too complex area source boundaries (i.e. area boundaries with inner spaces), and eliminating multiple segments or points.
After the necessary curation of the input files, for the final harmonization of the area, several questions had to be addressed by the core team (similarly with ESHM13 and other regional hazard models, Danciu et al 2018, Pagani et al 2018:

* In area sources with few events, how can we constrain the magnitude-frequency distributions (MFDs)? We observed that a significant number of seismogenic sources have no events of moment magnitude Mw > 3.5. In fact, about 35% of the area sources have less than 10 events, 7% of the area sources contain 10 to 15 events, 18% of the area sources have 30 to 60 events, and about 21% have greater than 60 events of Mw>3.5. This statistic considers the un-declustered earthquake catalogue without completeness constraints. Hence, to have robust earthquake statistics the use of super zones, or of intermediate layers, is required. Such super zones are constituted of a number of adjacent sources, assembled according to tectonic criteria. The well-defined magnitude-frequency distributions evaluated at the super zone level can help constrain the distributions in source zones with little data.

* Is there additional tectonic or geologic data that justify small sources with only few seismicity, or can a small source be merged with adjacent sources? A minimal number of earthquakes has to be captured within each area source, for a robust statistical estimate of the activity parameters. Area sources without known significant seismicity or known tectonic/geologic features should be considered as background sources i.e. sources with a minimum activity rate assigned as earthquake occurrence cannot be excluded.

Noteworthy, the area sources were further discussed with the local experts in bilateral meetings or workshops and their feedback and changes were iteratively integrated. Each area source model is characterized by various parameters including tectonic settings, style-of-faulting parameters, depth distribution, range of maximum magnitudes, lower and upper seismogenic depth, some of which are inherited from the ESHM13. The procedure to estimate the seismic activity parameters are given in the next section. **Fig. 3.1** illustrates the geographical distribution of the shallow crust area, as well as the deep seismicity and subduction inslab sources of Europe and Turkey.

### 1.1.1 The use of Super Zones

Super zones provide a spatial proxy to describe tectonic features, geological fault systems and seismicity patterns. Some parameters evaluated at the super zone level are then applied at the level of seismogenic sources. Super-zones cover large geographical regions, and we distinguish three types (see **Fig. 3.2**): completeness super zones (CSZ), tectonic super zones (TECTO) and the maximum magnitude zones (MAXMAG).

The delineation of the CSZ closely follows the regions used for the compilation and magnitude harmonization of the instrumental catalogue (EMEC, chapter 3 section 1.2). The primary use of the CSZ is to provide the geo-delineation of the regions for estimating the complete magnitude and time windows. Assessing the catalogue completeness is more feasible for super-zones than for individual seismogenic sources due to the increased number of earthquakes that improves the statistical analysis.

TECTO model consists of super zones extracted from the national hazard models and then updated for consistency with the unified earthquake catalogue, active faults and the tectonic regionalization given in Chapter 3 section 3.5. The main use of TECTO is to constrain regional seismicity by facilitating a robust estimation of the activity parameters given the larger number of earthquakes within each super-zone. The TECTO model also provides the basis for a smoothed seismicity model and delivers the regional b-values necessary for estimating the seismic productivity of active faults.

MAXMAG zonation is an update of the maximum magnitude super zones of ESHM13 (Meletti et al SHARE Deliverable D3.3). The main use of this geo-regionalization is the estimation of the maximum
magnitude (Mmax) based on the earthquake catalogue. It shall be noted that none of the super zone
models presented herein are considered as a direct input for seismic hazard calculation.

<figure><img src="image_21.png" /><figcaption>**Fig. 3.1** ESHM20: Area Sources Model to describe the shallow crust seismicity (left) and the deep seismicity and subduction inslab sources (right) across the pan European region</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 3.2** ESHM20 Completeness super zones (top left), tectonically delineated TECTO super zones (top right) and MAXMAG maximum magnitude super zones (bottom).</figcaption></figure>

1.1.2 Activity Rates Parameters and Magnitude Frequency Models

The activity rates of the individual area sources are estimated based on the assumptions that the regional seismicity follows a memoryless Poisson process characterized by a stationary mean rate of occurrence described by an exponential distribution i.e. Gutenberg- Richter model: Log₁₀ N = aGR – bGR*M. The N is the cumulative number of earthquakes per year equal to or greater than a magnitude M (in this model we use the moment magnitude Mw), and aGR and bGR are constants. The aGR-value
represents the total seismic productivity of a given source (the log of number of events with $M \ge 0$); the $b_{GR}$-value is the negative slope of the recurrence curve expressing the average ratio of exponentially distributed small and large magnitude earthquakes.

Earthquake recurrence parameters of each area source are derived from the declustered earthquake catalogue (see chapter 2, section 2.1) using the completeness intervals as defined for each CSZ-zonation (see chapter 2, Section 2.2). The activity parameters are estimated, accounting for the uncertainties of the aGR, bGR parameters, and their confidence intervals, for both data and prediction. Methodology is outlined in Bollinger et al (1989, 1993), Weichert (1980), Berril and Davis (1980) and Herrmann (1977).

The uncertainties of the $a_{GR}$, $b_{GR}$ parameters are defined by their covariance and the associated standard error. To estimate the range of $a_{GR}$, $b_{GR}$ values we used a random sampling technique. Firstly, a multivariate normal distribution is used for generating samples of correlated $a_{GR}$ and $b_{GR}$ values, based on their values, their covariance and the individual standard error.

Next, we bootstrap 1 million samples and the discrete approximation method of probability distributions (Miller and Rice 1983) is used to obtain the range of the $a_{GR}$, $b_{GR}$ values. The 16th, 50th, and 84th percentiles are used to represent the lower, median and the upper values of the $a_{GR}$ and $b_{GR}$, respectively. The activity parameters were estimated for the completeness regions (CSZ), for the large-scale tectonic zones (TECTO) and as well as for the area source model (ASM).

Examples of the uncertainty range of the magnitude frequency distribution for three sources in Spain (ESAS134), France (FRAS176) and Italy (ITAS317) are given in Figure 3. Note the different number of events used in each source ($n\_used$ is the number of complete events used for estimating the aGR, bGR values). There is a relevant correlation between the number of events used for statistical fitting to obtain the GRs parameters and the uncertainty range (green dashed curved in the same plot); the uncertainties decrease as the number of events increases and vice-versa.

The mfds of the upper and lower bGR-value are well constrained by observations for the source with fewer events (EAS134, i.e. 33 events) and decreases as the data increases (FRAS176, i.e. 78 events, ITAS317, i.e. 220 events).

Another critical factor affecting the statistical fitting of the seismicity parameters of the source models is the magnitude bin sizes (Stromeyer and Grunthal, 2015). The maximum likelihood method (Weichert 1980) is sensitive to the number of events in the most populated bins, usually the first magnitude bins, which often describe the frequency of modern instrumental seismicity.

Herein, the size of the magnitude bin is set up to 0.2 units of magnitude, and the statistical fitting will apply for a minimum of two consecutive bins with at least three events in each bin when available, otherwise a proxy on bGR from the overlaying tectonic super zones (i.e. TECTO) is used. Furthermore, for area sources with more than 30 complete events, the GR parameters are computed with the automatic maximum likelihood method using the earthquakes falling in the area source, whereas for the area sources with less than 30 events, $b_{GR}$ from TECTO is used (Danciu et al 2018, Grünthal et al 2018, Sesetyan et al 2018).
<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

**Fig. 3.3** Examples of MFDs (left panels) for individual area sources in Spain (top), France (mid) and Italy (bottom). The right panels illustrate the magnitude-time windows intervals for which the earthquakes are assumed to be complete.
In this automatic procedure, the estimation of the bGR-value is implemented so that first the statical
fitting is done with all the magnitude bins; if the resulting bGR is outside a predefined range of values (i.e.
0.70 to 1.15), then the first magnitude bin is removed and the bGR is recomputed and checked again
against the predefined bGR values.

The procedure allows for iteratively removing a maximum of five magnitude bins. If bGR does not converge to values within the validity range (0.7 to 1.15), then the bGR from the TECTO model is considered. The exceptions are the volcanic sources for which the bGR values greater than 1.15 are expected.

When a regional bGR from TECTO is used as a proxy for the area sources with few data, the seismic activity is estimated either from a direct estimation of the aGR or by re-scaling the occurrence rates as a function of the number of complete earthquakes.

Specifically, the scaling factor is the ratio between the number of complete events in the area source and the number of complete events in the corresponding TECTO super zone, i.e. the aGR value of an area with limited events is equal to log₁₀(N_tecto (M>0) * (number of complete events in the super zone / number of complete events in the area source)).

In area sources with zero complete events, the aGR is computed as a function of the regional (i.e. aGR of the super zone) seismic productivity normalized per units of area, aGR = (aGR<sub>TECTO</sub> - log10(tecto superzone-area in km2))*area source area in km2.

**Fig. 3.4** displays the MFDs of individual area sources in an example super zone in Greece. The sum of the area source MFDs is compared with the regional seismic productivity as given by the TECTO super zone MFD.

It shall be noted that there are cases when the area source contains enough events, yet the resulting bGR lies outside the validity range, and hence the bGR is set to the regional bGR from the TECTO model. Other exceptions worth mentioning imply long sequence of historical seismicity (i.e. Southern Italy and off-shore Portugal), and in this case a two-range GR model (Utsu 1999) is used to define such a distribution (see Fig. 3.5).

A double truncation Gutenberg-Richter magnitude-frequency distribution is used to characterize the decay of the annual occurrence rates with different magnitude bins (i.e. from lower to upper magnitude bounds).

Additionally, a tapered Pareto distribution (Kagan 1993, Utsu 1999) is considered as an alternative to the double truncated MFD. The main reason is to provide an alternative recurrence model with a faster decay of the rates towards the maximum magnitude. In many area sources without events above Mw5, the rates obtained in the magnitude interval 5 to 6 from the GR model could be inflated. In these cases, the tapered Pareto distribution provides alternative estimates for the rates in the moderate to large magnitude range.

A comparison between the two distributions, i.e. double truncated GR distribution and the Pareto tapered distribution model, are given in **Fig. 3.5** for an area source in Southern Spain. The differences of the two MFDs are obvious for the recurrence rates of Mw > 5.2, with the Pareto tapered distribution depicting a faster decay of the rates.

The tapered Pareto distribution requires the activity parameters aGR and bGR, and an additional parameter, the corner magnitude Mc, estimated here as a function of the observed maximum
magnitude. A default M corner =6.0 is used for area sources with an observed maximum magnitude lower than 5.5; otherwise the M corner is automatically estimated as a function of the maximum observed magnitude.

The evaluation of the maximum magnitude of area sources is based on the strategy defined in ESHM13 (Meletti et al 2010) using the updated tectonic super zones (i.e. see **Fig. 3.2** for MAXMAG zonation) and the unified earthquake catalogue.

Three Mmax values are used to represent the upper bounds of the expected magnitude: the lower bound of the maximum magnitude (Mmax01) is equal to the highest magnitude observed for the zone, increased of its standard deviation; the second value (Mmax02) is assumed equal to Mmax01 plus a magnitude increment of 0.3, and finally the third value is assumed to be conservatively equal to Mmax01 plus a magnitude increment of 0.6.

<figure><img src="image_26.png" /><figcaption>**Fig. 3.4** Comparison of the magnitude-frequency distribution of the TECTO super zone (TSZ020, in Northern Greece) with the sum of individual MFDs of area sources. The activity parameters of the TECTO zone are given in the header of the left panel; the completeness plot and the geometry of the area sources together with the activity parameters of the individual area sources are given in the right panel.</figcaption></figure>

The magnitude increment was statistically estimated from the entire unified catalogue; the value 0.3 presents the statistical uncertainty (i.e. standard deviation) of the moment magnitude in the unified earthquake catalogue. Wheeler (2009) states that two added magnitude increments would imply an extension of the historical record by a factor of 1.35 to 1.7, assuming a regional b-value of about 0.9 to 1.00. With the two factors, the maximum magnitude obtained would correspond to a historical record of about ~ 4000 and ~ 5000 years (Danciu et al 2018).

For the tapered Pareto distribution, the Mmax is controlled by the Mcorner, and its value might be different from those estimated in the super-zones, as it can be observed in **Fig. 3.6**. For a fixed annual cumulative rate (i.e. 10⁻⁴) the Mmax values of all distributions can be compared: the lower Mmax values for Pareto's MFD are 6.1, 6.3 and 6.5 whereas the range of the remaining Mmax values span from 6.5 to 7.1.

In summary, shallow crust area sources are characterized by specific earthquake recurrence (aGR- and bGR-values) parameters with associated uncertainties, and a distribution of Mmax values. Details of the model implementation are given in chapter 7.
<figure><img src="image_27.png" /></figure>

**Fig. 3.5** Illustration of the two-range magnitude frequency distributions (Utsu 1999) used for two area sources in southern Italy (upper plot) and off-Shore Portugal (lower plot).

<figure><img src="image_28.png" /></figure>

**Fig. 3.6** Example of an area source in Southern Spain, modeling of the recurrence based on two alternative magnitude-frequency distributions. The double truncated GR distribution (black line) and its uncertainties (dashed-blue lines) and the Pareto tapered distribution (green line) are given in the left panel.
## 3.2 Active faults and background seismicity model

### 3.2.1 Active Faults

The concept of active faults and background seismicity has become the preferred solution for regional seismogenic source models (see Danciu et al 2018 for EMME14, Sesetyan et al 2018 for Turkey, and global mosaic of hazard models of GEM, Pagani et al 2018).

The ESHM13 model for active faults and background area sources is abandoned in ESHM20. The main reasons are twofold: 1) background area sources were delineated to match the area source model, preventing the independency assumption of the alternative seismogenic source models, 2) uneven spatial completeness of the faults within the background area sources. In the updated version, we combine the active faults and background smoothed seismicity.

The fault-based source model is built upon the European Fault-Source Model 2020, an update of the European Database of Seismogenic Faults (EDSF) and several other recent active fault compilations in the Euro-Mediterranean regions (Basili et al 2020). The crustal faults are represented with a down-dip planar geometry that defines its three-dimensional plane in the Earth's crust. In this representation, the trace of the fault upper edge is extruded downward basing on the dip and depth values.

The minimum set of basic fault parameters required for constructing a seismogenic source model refer to Geometry (Location: Lat, Lon, Depth; Size: Length, Width; Orientation: Strike, Dip) and Behavior (Rake and Slip Rate). This set of attributes is associated with each geographic feature, and then a magnitude-frequency distribution (MFD) can be derived for each fault source.

To this end, the geologic information contained in the fault datasets can be used to estimate two key parameters: 1) the seismic productivity by converting the long-term geological or geodetic slip rates into activity rates, assuming moment conservation, and 2) a proxy for the moment magnitude upper bound of the MFD. It is assumed that the geological slip rates encompass several seismic cycles of large-magnitude earthquakes on a fault, offering advantages over historical seismicity when used to estimate earthquake frequencies. Furthermore, another important assumption is that the fault location provides a spatial proxy for occurrence of future large magnitude earthquakes (M>6.0).

In the ESHM20, we use the Model 2 (i.e. arbitrary area) of Anderson and Luco (1983) that is a truncated exponential magnitude-frequency distribution to characterize the seismic activity of each entry of the active faults. The functional form of the recurrence relationships is provided in the original paper of Anderson and Luco (1983), and not repeated here.

For all faults, the following values were considered as default: crust shear modulus (μ) value of 3.0 x 1e11 dyne/cm2; typical values (c=16.05 and d=1.50) for the magnitude-moment scaling coefficients as originally proposed by Kanamori and Anderson (1975).

A key assumption when using the active faults is that the seismic productivity of these faults will satisfy the GR relation. In this view, we use the bGR-values of the TECTO-zones that host the faults. The TECTO super zones are representative of regional seismicity. The active faults, although simplified herein, are not isolated, but form branching structures and networks.

The maximum magnitude of each fault source is obtained by use of magnitude - fault geometry scaling relationships, hereinafter referred to as fault scaling laws (FSL).
The generalized functional form between rupture dimensions (L, W, A, D) and moment magnitude (Mw) is $(X) = a + bM_w$, where X is the rupture dimension under consideration (e.g., L or W in km, A in km2, D in m) and the coefficients a and b are empirically determined. We here propose to adopt the FSL developed by Leonard (2014) for crustal ruptures and by Allen and Hayes (2017) for subduction interface ruptures.

A sanity-check revealed discrepancies between the fault capabilities to accommodate the assigned $M_{max}$. Specifically, there are small faults whose average area may be smaller than the rupture area associated with the higher $M_{max}$ values when considering all the variability of the fault parameters and the uncertainty in the FSL; as such, the rupture corresponding to this $M_{max}$ might extend beyond the fault geometry limits.

An analysis of the effects of various key parameters of seismic productivity indicates that two parameters are responsible for a large variability of the magnitude frequency distribution: the slip-rates used and the $M_{max}$ (see Fig. 3.7).

For a constant slip rate, increasing the $M_{max}$ will result in decrease in the recurrence rates of low-to-moderate magnitude events, this is because the large magnitude events constitute the majority of the total seismic moment rate; while increasing the $M_{max}$ will require a subtraction of many smaller earthquakes to preserve the same seismic moment budget (Youngs and Coppersmith 1985).

The fault productivity is also sensitive to the fault area, the larger the area the lower the seismic productivity and vice-versa. The effect of the fault area is correlated with the effect of the $M_{max}$, thus only the latter is retained.

To conclude, the main input parameters are the bGR-values, fault-slip rate and the maximum magnitude ($M_{max}$) and their uncertainty on the seismic productivity of each fault is captured as such: three values of the slip-rates, three values for the maximum magnitude, a single b-value (estimated from robust statistics in the TECTO zones) and an average area of each fault.

<figure><img src="image_29.png" /><figcaption>**Fig. 3.7** Effects of various key input parameters in seismic productivity of a single fault: maximum magnitude (i.e. two values of Mw=7.0 and 7.3) and variable slip-rates values (5mm/year and 10mm/year)</figcaption></figure>
3.2.2 Smoothed Seismicity

The smoothed seismicity model is based on the TECTO zonation, for which the activity parameters (aGR- and bGR-values) were estimated on the declustered and complete earthquake catalogue using the smoothing algorithm of Nandan et al (2021, *EFEHR Technical Report, in preparation*).

Several choices are required to make earthquake forecasts based on the smoothed seismicity approach. Critical decisions include choosing the (1) adaptiveness of the bandwidth of the kernel, (2) smoothing parameters, (3) declustering algorithm, and (4) declustering parameters.

Training and validation sets are typically required to determine the best possible combination of options. The training set is used to generate the smoothed seismicity model for a given combination of the options mentioned above, while the validation set is used to rank the options in terms of mutual information gain.

We first propose a penalized-maximum-likelihood system that allows us to find optimal smoothing parameters using only the training set without invoking the need for a validation set. We then set up a pseudo-prospective test bench to validate the different combinations of choices used to obtain the smoothed seismicity forecasts.

In this test bench, we use all the data through 2006 to train a series of smoothed seismicity models (SSMs) and then compare their forecasting ability using the earthquakes between the entirety of 2007 through 2015. Based on this validation, we recommend the hyper-parameter choices for obtaining an alternative SSM for ESHM20. Finally, we perform a long-term retrospective validation of SSM based on the reference declustering algorithm, which is used in ESHM20, to demonstrate its consistency with the observed number of earthquakes in different TECTO polygons.

<figure><img src="image_30.png" /><figcaption>**Fig. 3.8** Conceptual sketch of active faults and background buffers for removing the double counting of seismicity.</figcaption></figure>
### 3.2.3 Combining the Smoothed Seismicity with the Active Faults

In the vicinity of active faults, the background-smoothed seismicity is complementary to the seismic productivity together with the faults, whereas in regions with no faults, the smoothed seismicity is an alternative to the area source model (Northern Europe).

The main challenge when modeling the background seismicity with the active faults is the avoidance of double counting the seismicity in the fault's proximity. To avoid the double counting of the earthquake rate forecast due to active faults and background seismicity, a variable threshold magnitude was introduced. This concept requires a definition of fault proximity, which in this case is defined as a spatially symmetric buffer following the top-trace of the fault. The buffer zone results from combining the fault-width projection to the Earth surface and an arbitrary polygon that extends up to 5km along strike on the opposite side of the fault trace (see Fig. 3.8).

The magnitude threshold ($fMthr$) is fault dependent, and the main assumption is that events with $Mw > fMthr$ are associated with the faults and a minimum of two magnitude bins is imposed for each fault. The $fMthr$ for each fault is based on the Mmax ($fMaxMag01$) value and is estimated as following:

- faults with $fMaxMag01 \le 6.3, fMthr=5.9$

- faults with $6.3 \le fMaxMag01 \le 7.3, fMthr= fMaxMag01-dM$, where the $dM=0.6$ (the overall uncertainty of the $M_{max}$)

- faults with $fMaxMag01 > 7.3, fMthr=6.7$

Thus, within the fault buffer the occurrence rates of lower magnitude bins ($\le fMthr$) are represented by point sources, whereas the occurrence rates of the $M>fMthr$ are associated with the fault. Outside the buffer, the $Mmax$ of the smoothed seismicity is the $Mmax$ value attributed to the corresponding TECTO zone.

Examples of comparisons between the composite cumulative annual rates of the active faults, background sources and the rates observed are done at the TECTO level and given in **Fig. 3.9**. As expected, there are three cases: 1) total productivity of all active faults matches the observed earthquake rates (TECTO zone TSZ041, Italy); 2) total productivity of all active faults overestimates the observed earthquake rates (TECTO zone TSZ029, Southern Spain) and 3) total productivity of all active faults underestimates the observed earthquake rates (TECTO zone TSZ022, South West of Bulgaria).

There are various possible reasons for such discrepancies: over or under estimation of the fault slip-rates, completeness of the active faults within each region, subjectivity in the delineation of the TECTO zone, incompleteness of earthquake catalogue at large magnitudes, short observation period for the earthquake catalogue, etc.

This visual comparison aids the sanity check of the seismic productivity of a specific seismotectonic domain. We consider that the uncertainties of the future earthquake rates are depicted by both the area source model and the smoothed seismicity plus fault model.
TSZ029 : 43 faults ($a_{GR}=4.74$) and BG seismicity

Tecto TSZ029 1082<sub>decli</sub>:492<sub>Mc</sub>:492<sub>used</sub>

$a_{GR}=4.58$, $b_{GR}=1.01$

$M_{max}$ = [7.1 7.4 7.7]

<figure><img src="image_31.png" /></figure>

TSZ029 Completeness

<figure><img src="image_32.png" /></figure>

TSZ022 : 24 faults ($a_{GR}=2.91$) and BG seismicity

Tecto TSZ022 106<sub>decli</sub>:63<sub>Mc</sub>:0<sub>used</sub>

$a_{GR}=2.70$, $b_{GR}=0.74^{\circ}CSZ_{val}$

$M_{max}$ = [7.3 7.6 7.9]

<figure><img src="image_33.png" /></figure>

TSZ022 Completeness

<figure><img src="image_34.png" /></figure>

TSZ041 : 24 faults ($a_{GR}=4.73$) and BG seismicity

Tecto TSZ041 434<sub>decli</sub>:250<sub>Mc</sub>:236<sub>used</sub>

$a_{GR}=3.33$, $b_{GR}=0.80$

$M_{max}$ = [6.9 7.2 7.5]

<figure><img src="image_35.png" /></figure>

TSZ041 Completeness

<figure><img src="image_36.png" /></figure>

**Fig. 3.9** Comparison of total seismic productivity of all faults within a TECTO super zone, examples in Spain, Bulgaria and Italy. The left panel shows the summed MFDs (in blue) of all faults (in green) within the TECTO zone; the right panel indicates the complete events within each TECTO zone. The colored curves are the mfds for individual faults as described by their $a_{bGR}$ parameters in the middle panel, the green solid line are the summed mfd of faults and background seismicity, the background seismicity mfd is depicted in a dashed green line, the black line indicates the TECTO mfd, whereas the weighted mfd between the TECTO and faults plus smoothed seismicity is described by the red curve.
## 3.3 Subduction and Deep Seismicity Sources

### 3.3.1 Subduction Interface

Subduction sources cover the parts of the Europe-Africa plate boundary where there is still the presence of oceanic crust in the lower plate. They are the Gibraltar Arc, Calabrian Arc, Hellenic Arc, and Cyprus Arc (Figure 8). The zones of active continental collision are included in the dataset of crustal faults.

After the EDSF compilation from project SHARE (Basili et al., 2013), the most recent slab geometrical reconstructions are those for the Calabrian Arc by Maesano et al. (2017) and the global compilation by Hayes et al. (2018), as well as a number of studies that treated specific aspects or parts of these subduction zones. The included subduction sources have largely been inspired by these studies.

The subduction interfaces are modelled as complex faults whose geometry is sampled from the 3D models shown in **Fig. 1.7** considering estimates of the upper and lower boundary of the seismic interface. These same models, along with estimates of the slab thickness, are also used to derive the intraslab sources which are cut at a maximum of 300 km depth.

The recurrence parameters of the subduction zones (interface and intra-slab) are estimated based on the unified earthquake catalogue. Moreover, for the seismic interface sources we also used an alternative recurrence model based on tectonic convergence rates across the subduction zones and the double-truncated frequency-magnitude distribution model.

An estimate of the activity rates is obtained through an exploratory logic tree considering variable upper and lower seismogenic depths, magnitude upper bounds, depth-dependent rigidity, tectonic rates based on the geodetically observed convergence across subduction zones, seismic efficiency, and b values in an approach similar to Davies et al. (2019).

The average values of convergence rates are 1 mm/year for the Gibraltar Arc, 2.9 mm/year for the Calabrian Arc, 24 mm/year for the Hellenic Arc, and 12 mm/year for the Cyprus Arc (see Section 3.3.2). The magnitude upper bound of the magnitude-frequency distribution is constrained by the scaling relations of Allen and Hayes (2017), assuming that the entire seismic interface area may rupture in one single event.

The key parameters can be classified in two components: geometry related (the first three branching levels) and seismogenic parameters: the tectonic rates, seismic efficiency and regional bGR-values. The aim is to estimate the total variability for the number of events generated by the interface when cascading all key input parameters.

The conversion of the tectonic rates into seismic productivity was done with the model of Kagan (2002). The resulting distribution of activity rates consists of 2187 possible realizations of activity rates. However, such distribution is extremely computationally demanding for the hazard computation, thus the logic tree was trimmed down to 27 branches by combining the first two branching levels and a median crust rigidity model to describe the geometry and crust features; a single b-value is used, as well as three values of the tectonic rate and seismic coupling.

The values of the key input parameters are given in the table of **Fig. 3.10**. The maximum magnitude is coupled with the effective area of the subduction interface. The resulting magnitude frequency distributions are illustrated in **Fig. 3.11**.
<figure><img src="image_37.png" /></figure>

<table><caption>**Fig. 3.10 (Top)** Map view of the slab models, GiA - Gibraltar Arc, CaA - Calabrian Arc, HeA - Hellenic Arc, and CyA - Cyprus Arc. (Bottom) Key input parameters for the calculation of the magnitude-frequency distribution.</caption><thead><tr><th>Slab</th><th>Upper Seism. Depth (km)</th><th>Lower Seism. Depth (km)</th><th>Conv. Rate (mm/yr)</th><th>Seism. Eff.</th><th>b-value</th></tr></thead><tbody><tr><td>GiA</td><td>8-10-11</td><td>21-27-33</td><td>0.8-1-1.2</td><td>0.1-0.2-0.3</td><td>1-1-1</td></tr><tr><td>CaA</td><td>9-10-12</td><td>19-34-39</td><td>2.7-2.9-3</td><td>0.1-0.2-0.3</td><td>1-1-1</td></tr><tr><td>HeA</td><td>9-10-11</td><td>30-34-38</td><td>19-23.8-28.6</td><td>0.1-0.2-0.3</td><td>1.05-1.07-1.15</td></tr><tr><td>CyA</td><td>6-8-10</td><td>25-31-38</td><td>9.8-12.2-14.6</td><td>0.1-0.2-0.3</td><td>1.05-1.07-1.15</td></tr></tbody></table>
<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /><figcaption>**Fig. 3.11 (Top)** Alternative components used to derive a distribution for the activity rates of the subduction interface sources. (Bottom) Frequency-magnitude distributions of the four slab interfaces. The series of numbers in the distribution labels of each panel represents the selected component at each branch of the exploratory logic tree.</figcaption></figure>
3.3.2 Deep Seismicity: Vrancea, Romania

A particular case of subduction slabs is the deep seismicity observed beneath the Vrancea Region, Romania. We keep the same seismogenic model (Danciu and Wössner 2014) proposed in the ESHM13 as illustrated in **Fig. 3.12**. The most hazardous seismicity in the Vrancea region, situated beneath the southern Carpathian Arc, is located at a depth between 70-150km. As there is currently no well-defined complex fault structure (no well-defined 3D fault plane), we model the pattern with four zones located at depth and a different zonation for the upper crust.

The seismicity parameters (GR parameters) are derived from the declustered catalogue and the data-
driven derived completeness intervals applied to earthquakes occurring at a hypocentral depth greater
than 70km.

Maximum magnitude for the seismicity in the Vrancea region was assessed by the historical record, i.e. we added an uncertainty to the maximum observed event and then added further increments with decreasing weights in the computation. The Mmax-values for the seismicity at depths 70-90km depth are 7.5 and 7.8, the seismicity between 90-150km has Mmax values of 7.8 to 8.1.

The magnitude frequency distributions are the same as those used for the shallow crust regions, an exponential double truncation GR and tapered Pareto model, including the uncertainties.

ESHM20: Deep Seismicity Sources, Vrancea, Romania

<figure><img src="image_40.png" /></figure>

© ESHM20 - EFEHR

<caption>Fig. 3.12 Deep Seismicity Sources – Azimuth dependent in Vrancea Region, Romania</caption>
## 3.4 Seismogenic Source Model - Logic Tree

Development of a seismic hazard model requires decisions and wide interpretations at every step to overcome limitations on data, randomness associated with input data, incomplete knowledge as well as challenges due to large scale modeling and calculations. One of the key constraints on the development of the logic tree was to design a logic tree that is feasible for the hazard calculation at the entire scale of the pan European region including Turkey.

In the design of the logic tree, we acknowledge the circular role of various inputs and components which might be reflected in a correlated logic tree branch. The logic tree is made of two main branching levels, area source model and smoothed seismicity plus fault model (Fig. 3.13). Both branching levels are based on the unified earthquake catalogue, the completeness zones, and regional bGR proxies from the TECTO model. The active fault model and the area source model are treated as independent, albeit in the area source delineation fault locations were used in some parts of the model. The weight assigned to a logic tree branch should reflect the independence between datasets, the assumptions taken on the model building and ultimately it shall represent the likelihood of that branch being the true model. Herein, weights are assigned by the experts deriving the input models, based on their judgement of the likelihood of the alternative models.

The logic tree structure of the seismogenic source model is asymmetric. The first branching level depicts the two main regional source models: the area sources model and the active fault with background smoothed seismicity model. The second branching level addresses the individual area sources for which two alternative magnitude-frequency distributions are handled as individual branches: one for the exponential double truncation GR and one for the tapered Pareto distribution. At the level of individual area sources, the main branching levels are the a and b-values and maximum magnitudes for area sources, whereas for active faults the uncertainties on the slip-rates and on the maximum magnitudes are considered.

We assume a constant seismic rate (above a given magnitude) when applying the maximum magnitude uncertainties of the area sources, and a constant seismic moment rate when considering maximum magnitude uncertainties of active faults (see Figure 2 in Youngs and Coppersmith 1985).

The uncertainty of the activity parameters for individual area sources is obtained from a simulated distribution of the activity parameters using the standard deviation of the aGR, bGR and the covariance matrix of these two parameters.

The median, and the two quantiles, i.e. 5th and 95th percentiles, populate the three independent branches of individual area sources for the exponential double truncated GR distribution. The weights assigned to the lower, median, upper activity parameters are 0.2, 0.6 and 0.2 respectively, representing the 5-95th confidence interval. For the tapered Pareto distribution only the median aGR, bGR values are used.

Three branching levels for the Mmax are symmetrically applied to all models, however with different weights. The Mmax weights for the area sources are the same as in the ESHM13 in the case of exponential double truncated GR distribution and equal weights for the tapered Pareto distribution.

For the active faults, only one value is used for the bGR-value (from TECTO) because it was found that the impact of b-value uncertainties on the magnitude-frequency-distributions is lower in comparison with the impacts of the fault slip-rates and Mmax uncertainties.
Note that, each seismic source is characterized by alternative values of various source parameters (i.e. average depth and the predominant style-of-faulting) that were not treated as epistemic but rather aleatory.

<figure><img src="image_41.png" /><figcaption>**Fig. 3.13** ESHM20: Schematic representation of the logic tree of the seismogenic sources for shallow crust earthquakes.</figcaption></figure>

For deep seismicity of Vrancea, Romania and the subduction intraslab, the upper part of the logic tree is used with the same weights as indicated in Figure 13. To conclude, for the final calculation, the main seismogenic source model logic tree is supplemented with the subduction sources (both interface and in-slab) and Vrancea region in Romania for a complete model, and also combines ground motion logic tree model described in the next Chapter.
# 4 Ground Motion Model Characterization

## 4.1 Ground Motion Modelling Approach

The development of the ground motion model (GMM) logic tree for a region should aim to leverage as far as possible on the available strong motion data in order to characterize the expected ground motion from earthquakes, its aleatory variability, its epistemic uncertainty and its variability from region to region.

Previous models, including the 2013 European Seismic Hazard Model (ESHM13), have mostly attempted to do this by adopting a multi-model approach, in which a set of *best suited* GMMs are identified for each tectonic region of relevance (e.g., active shallow crustal seismicity, subduction interface, subduction in-slab, volcanic etc.) and weighted in a manner that reflects either their respective fits to the observed data (where available) or to expert views on suitability. Though practical, this approach has many theoretical limitations when it comes to their objectives of characterizing the epistemic uncertainty in the ground motions (Atkinson et al., 2014).

In the ESHM20 we adopt a different approach for characterising epistemic uncertainty, which is built around the concept of a *scaled backbone* ground motion model logic tree. In this approach a single ground motion model is calibrated (or selected from the literature) and to this model adjustment factors are applied that quantify the uncertainty in the expected ground motion as a result of the limited knowledge on the seismological properties in a region.

We not only quantify epistemic uncertainty in terms of a parameterised statistical distribution, we also set in place a framework through which new information and/or information at regional scale can adapt the backbone GMM logic tree and reduce the epistemic uncertainty over time. In the following we describe the scaled backbone ground motion model logic tree for the three main seismotectonic region types in Europe: shallow crustal seismicity (non-craton), seismicity in the stable craton region of north-eastern Europe, and subduction and deep seismicity (including the Hellenic, Calabrian, Cypriot and Gibralter arcs, as well as the Vrancea deep seismic zone). In addition, a number of *special cases* are considered in which the main approach is modified, or else decisions are made on the basis of insights and information from other data sets or studies.

A more comprehensive discussion of the theoretical and practical aspects of the ESHM20 GMM logic tree can be found in an accompanying EFEHR Report (Weatherill et al. 2022a) and in multiple journal publications that have arisen from this work (Kotha et al., 2020, 2022; Weatherill et al, 2020; Weatherill & Cotton, 2020; Weatherill et al.; 2022b).

## 4.2 Ground Motion Model Logic Tree for Shallow Seismicity

### 4.2.1 A regionally adjustable backbone Ground Motion Model (GMM) and logic tree

The backbone GMM logic tree for shallow crustal seismicity leverages heavily on the large volume of information in the Engineering Strong Motion (ESM) flatfile (Lanzano et al., 2019). With the available information a backbone ground motion model is constructed in which regional differences are quantified explicitly as random effects within the mixed effects regression process. Here, two prior regionalization models based on geology and tectonics are used to assign the earthquakes and the recording stations to their respective zones.

The two regionalisations (shown in Fig. 4.1) are i) a tectonic zonation based on an early version of the TECTO source model, and ii) a geological and tectonic zonation provided by Basili et al. (2019) and with small adaptations in the zone boundaries to improve the distribution of recording stations. Earthquake
sources are assigned to regions according to the TECTO model, while stations are assigned to the regions from the Basili et al. (2019) zonation.

<figure><img src="image_42.png" /><figcaption>Fig. 4.1 Regionalisation models used to calibrate the region-to-region random effects: i) TECTO model (left) and ii) adapted geological and seismotectonic model from Basili et al (2019) (right)</figcaption></figure>

For the development of the model, records are selected from the ESM database according to the following criteria: i) events must be classified as non-subduction (according to the classification scheme indicated in the subsequent section), ii) reported hypocentral depth must be shallower than 40 km, iii) events must have ≥ 3 records, iv) only records with high-pass filter frequency $f_{hp} \le 0.8 T$ are retained for each period in the regression. The final database yields 18,222 records from 927 events ($3.1 \le M_W \le 7.4$) recorded at 1829 stations ($0 \le R_{JB}(km) \le 545$). The predictor variables of the model are moment magnitude ($M_W$), taken from the updated EMEC catalogue (section 3), hypocentral depth ($z_h$) in km, Joyner-Boore distance (km) and 30 m averaged shearwave velocity ($V_{S30}$, m/s). The Joyner-Boore distance is taken as the distance to the surface projection of the rupture plane where available (for most events with $M_W \ge 5.5$) or treated equivalently to epicentral distance otherwise. More details of the selection and data processing are given in Kotha et al. (2020).

The functional form of the backbone model is (Kotha et al. 2020):

$$ \ln(\mu_{ijkl}) = e_1 + f_{R,g}(M_{W,i}, R_{JB,ij}) + f_{R,a}(R_{JB,ij}, k) + f_M(M_{W,i}) + \delta L 2 L_l + \delta B_{e,i}^0 + \delta S 2 S_j + \delta W_{e,ij} \quad (\text{Eq. 6.1}) $$

$$ f_M(M_{W,i}) = \begin{cases} b_1 \cdot (M_W - M_h) + b_2 \cdot (M_W - M_h)^2 & M_W \le M_h \\ b_3 \cdot (M_W - M_h) & M_W > M_h \end{cases} \quad (\text{Eq. 6.2}) $$

$$ f_{R,g}(M_{W,i}, R_{JB,ij}) = (c_1 + c_2 \cdot (M_W - M_{ref})) \cdot \ln \sqrt{\frac{(R_{JB}^2 + h_D^2)}{(R_{ref}^2 + h_D^2)}} \quad (\text{Eq. 6.3}) $$

$$ f_{R,a}(R_{JB,ij}, k) = \frac{c_3 + \delta c_{3,k}}{100} \cdot \left( \sqrt{(R_{JB}^2 + h_D^2)} - \sqrt{R_{ref}^2 + h_D^2} \right) \quad (\text{Eq. 6.4}) $$

Where $\mu_{ijkl}$ is the PGA or spectral acceleration at period T (Sa(T)) from event i located in event region l, recorded at station j located in station region k. The definition of the horizontal component of motion is RotD50 (Boore, 2010). $M_h$, $M_{ref}$ and $R_{ref}$ are period-independent constants, which take the values $M_W$ 6.2, $M_W$ 4.5 and $R_{JB}$ 30 km respectively in Kotha et al. (2020); however, the hinge-magnitude ($M_h$) was revised in the recent paper of Kotha et al. (2022) to $M_W$ 5.7 based on updated analysis of ground motion records from larger magnitude events in Europe and globally. $h_D$ is the near-source saturation plateau, a term that is dependent on hypocentral depth ($z_h$). This is fixed to a reference value depending whether the depth falls in one of three bins: $z_h \le 10$ km, $10 < z_h(km) \le 20$ or $z_h(km) > 20$, for $h_D$
takes the value 4 km, 8 km and 12 km respectively. Robust linear mixed effects regression is used, which down-weights outliers when fitting the fixed effects (e₁, b₁, b₂, b₃, c₁, c₂ and c₃).

In addition to the fixed effects, the model contains five random effects: $\delta L2L_l$, a source-region specific scaling parameter (depending on the location of the event), $\delta B_e^0$, the between-event residual filtered for source region-to-source region variability, $\delta S2S_j$ the site-to-site variability and $\delta c_{3,k}$ the between-region variability of the residual attenuation term $f_{R,a}$. Finally, $\delta W_{e,ij}$ describes the remaining region-, event- and site-corrected variability, representing therefore the apparent aleatory variability of the model.

These five random effects are described by Gaussian distributions such that: $\delta L2L_l = N(0, \tau_{L2L})$, $\delta B_e^0 = N(0, \tau_0)$, $\delta S2S_j = N(0, \phi_{S2S})$, $\delta c_{3,k} = N(0, \tau_{c3})$ and $\delta W_{e,ij} = N(0, \phi_0)$. We exploit the separation of region-specific random effects ($\delta L2L_l$ and $\delta c_{3,k}$) from the more traditional between-event, site-to-site and site-corrected within event random effects ($\delta B_e^0$, $\delta S2S_j$ and $\delta W_{e,ij}$) in order to be able to identify and apply regional calibrations of the GMM, and at the same time removing some of the ergodicity that would otherwise result in higher aleatory variability, and thus higher hazard.

The spatial distributions of the source-region to source-region random effects ($\delta L2L_l$) and the residual attenuation region random effects ($\delta c_{3,k}$) are shown in **Fig. 4.2** and their total variability in **Fig. 4.3**. Where $\delta L2L$ is higher we expect earthquakes in these regions to produce characteristically larger ground motions than the “average” across Europe, and where it is lower, we would therefore expect weaker ground motions. Similarly, $\delta c_{3,k}$ controls the rate of decay of ground motion over longer distances ($R_{JB} \ge 80$ km), but as $c_3$ is negative then a higher $\delta c_{3,k}$ results in a slower modelled attenuation, while a lower $\delta c_{3,k}$ indicates a faster attenuation.

<figure><img src="image_43.png" /><figcaption>**Fig. 4.2** Distribution of source-region specific residual $\delta L2L_l$ (left) and net residual anelastic attenuation term ($c_3 + \delta c_{3,k}$) (right) for Sa (0.2 s).</figcaption></figure>

The described by the terms $\tau_{L2L}$ and $\tau_{c3}$ effectively describe the full region-to-region variability in source and residual attenuation implied by the available strong motion data, conditional upon the assumed regionalization. As can be seen in **Fig. 4.2**, however, there are many regions for which $\delta L2L_l$ and $\delta c_{3,k}$ have not been determined owing to insufficient data. In this case we do not know the specific source-region or residual attenuation region term, but rather we assume that the probability distribution of these adjustments is described by the full region-to-region variability quantified by $N(0, \tau_{L2L})$ and $N(0, \tau_{c3})$. We use this information to create the *default* scaled backbone ground motion logic tree to be applied in regions where data are absent, but we believe the probability of a given $\delta L2L_l(T)$ and $\delta c_{3,k}(T)$ values follow the aforementioned distributions. The distributions themselves are mapped into branches of the logic tree using the N-branch discrete approximation to a Gaussian distribution proposed by Miller & Rice (1983). This approximation will be used throughout the logic tree, so the corresponding weights and $\epsilon$ values are for the N = 3, 5 and 7 branch cases as provided in Table 1.
<figure><img src="image_44.png" /><figcaption>**Fig. 4.3** Distribution of region-specific random effects for each spectral period: $\delta L_{2L_1}$ (left) and $\delta c_{3,k}$ (right). Markers are colour-scaled according to their resulting weight in the robust mixed effects regression, and the respective standard deviations shown in black line.</figcaption></figure>

The default backbone GMM logic tree has two branch sets, which correspond to the residual attenuation region and the source region variability respectively. This is shown in the top half of **Fig. 4.4**, in which the source-region distribution is mapped into 5 discrete branches, and the residual attenuation region distribution into 3 branches.

In the case of the source-region variability we take the larger of two values: $\tau_{L2L}$ and $\sigma_{\mu}$, the latter of which corresponds to the statistical uncertainty in the regression determined according to the approach of Al Atik & Youngs (2014).

By taking the larger of $\tau_{L2L}$ and $\sigma_{\mu}$ we ensure that over the range of $M_W$ and $R_{JB}$ scenarios that are well constrained in the regression (for which $\sigma_{\mu}$ is small) the epistemic uncertainty is not less than the source-region to source-region variability.

At the same time, for scenarios that are poorly constrained in the regression, typically large $M_W$ and small $R_{JB}$ the statistical uncertainty is then dominant. This is intended to minimize the degree of double counting of epistemic uncertainty where region-to-region variability may be present within the statistical uncertainty.

<table><caption>**Table 6.1:** Weightings and an $\varepsilon$ values for the 3-, 5- and 7-branch discrete approximations to a Gaussian distribution based on Miller & Rice (1983).</caption><thead><tr><th colspan="3">N = 3</th><th colspan="3">N = 5</th><th colspan="3">N = 7</th></tr><tr><th>Weight</th><th>ε</th><th>Weight</th><th>ε</th><th>Weight</th><th>ε</th><th>Weight</th><th>ε</th><th>Weight</th><th>ε</th></tr></thead><tbody><tr><td>0.167</td><td>-1.732051</td><td>0.011257</td><td>-2.85697</td><td>0.000548</td><td>-3.750440</td><td>0.666</td><td>0.000000</td><td>0.222076</td><td>-1.355626</td><td>0.030757</td><td>-2.366759</td></tr><tr><td>0.167</td><td>1.732051</td><td>0.533334</td><td>0.000000</td><td>0.240123</td><td>-1.154405</td><td>0.167</td><td>1.732051</td><td>0.222076</td><td>1.355626</td><td>0.457144</td><td>0.000000</td></tr><tr><td></td><td></td><td>0.011257</td><td>2.85697</td><td>0.240123</td><td>1.154405</td><td></td><td></td><td></td><td></td><td>0.030757</td><td>2.366759</td></tr><tr><td></td><td></td><td></td><td></td><td>0.000548</td><td>3.750440</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>
<figure><img src="image_45.png" /></figure>

<figure><img src="image_46.png" /><figcaption>**Fig. 4.4** Full scaled backbone ground motion model logic tree for shallow crustal seismicity in Europe, excluding the stable shield region.</figcaption></figure>

Explicit within the formulation of the *default* backbone GMM logic tree is the assumption that the epistemic uncertainty in the source-region and residual attenuation region is described by their respective full region-to-region variabilities.

For regions where data are available and region-specific values have been determined, however, the local adjustments to the backbone GMM and their respective uncertainties should be used in order to regionalize the GMM and reduce the epistemic uncertainty. In the case of $\delta L2L_l$ the regional variability has been explored and interpreted, and we find that not only are the specific spatial trends difficult to reconcile with previous observations of regional variation in ground motion in Europe, for many regions the values are constrained by very few events, and in most of these cases the events may come from a single seismic sequence.

The decision is taken *not* to adopt the local $\delta L2L_l$ distributions but rather maintain the full $\tau_{L2L}$ distribution. This should of course be revisited in future when more data are available and could still be re-considered at a local level should specific countries or users wish to adopt the local values.

In contrast with $\delta L2L_l$, our assessment of $\delta_{c3,k}$ led us to conclude that not only are these values better constrained by data, with more stations per calibrated region, they are also consistent with previous investigations into regional variability of strong motion attenuation within Europe (e.g., Kotha et al., 2016; Scasserra et al., 2009).

In this case a region-specific $\delta_{c3,k}$ term has been calculated for 42 zones (from 111 considered within the original regionalization), two of which are subsequently removed as they were found to be
constrained by only one station and their statistical uncertainty was almost equal to the total region-to-region variability.

These zones are Southern Betics (Spain) and Central Europe (covering Germany, Poland, Czechia, Slovenia and Denmark). From the 40 remaining zones, however, there are many regions showing similar period-dependent trends in $\delta_{c3,k}$, many of which are associated with moderate statistical uncertainty. The data show clear evidence of broader scale regional trends in attenuation, so rather than define a zone-specific $\delta_{c3,k}$ distribution we apply hierarchical cluster analysis to group these into five clusters (or cluster regions) that capture the main regional variations in attenuation.

Within each cluster region the period-dependent $\delta_{c3,k}$ values are grouped and using a Bayesian regression analysis (with the original distribution of $\tau_{c3}$ as a prior) we define a cluster region specific distribution of $\delta_{c3}$, which we define as $N(\delta_{c3,clp}, \tau_{c3,clp})$ where clp refers to cluster region $p$ for $p = 1, 2, ..., 5$. The spatial distribution of the clusters and the distributions of $c_3$ are shown in Fig. 4.5. Further details on this process can be found in Weatherill et al. (2020a).

<figure><img src="image_47.png" /><figcaption>**Fig. 4.5** Spatial distributions of the five residual attenuation cluster regions (left) and their corresponding distributions of $c_3$ (right). Figure from Weatherill et al. (2020)</figcaption></figure>

With the five cluster-specific distributions established and the regional trends in residual attenuation identified from the available data, the set of adjustments made to the scaled backbone ground motion model are then cluster specific.

This means that the N-branch logic tree illustrated in the top half of Fig. 4.4 is then repeated for each of the five cluster regions, with $c_3$ and $\varepsilon \cdot \tau_{c3}$ replaced by $c_{3,cl,p}$ and $\varepsilon \cdot \tau_{c3,clp}$ for $p = 1, 2, ..., 5$ respectively, while the default logic tree with full region-to-region variability is applied to those regions indicated by the semi-transparent white shaded regions in Fig. 4.5. The complete logic tree for shallow crustal seismicity excluding the stable shield region of north-eastern Europe, is shown in Figure 6.6.

### 4.2.2 Aleatory Uncertainty

As the region-to-region variability in source and attenuation characteristics retrieved from the mixed effects regression are used to construct the GMM logic tree, these terms are removed from the aleatory variability of the scaled backbone GMM. The resulting total aleatory uncertainty ($\sigma_T$) is therefore described by $\sigma_T = \sqrt{\tau_0^2 + \phi_0^2 + \phi_{S2S}^2}$. The corresponding between-event, site-to-site and site-corrected within event standard deviations provided by Kotha et al. (2020) are homoscedastic, i.e., independent of the scenario.
Many studies of ground motion variability in recent years, however, have identified heteroskedasticity in some (or all) of these residual terms (e.g., Rodriguez-Marek et al. 2013; Al Atik, 2015). Re-analysis of the uncertainty terms by Kotha et al. (2022) also showed that the assumption of homoskedasticity with respect to earthquake magnitude was not consistent with the data at large magnitudes.

The change in variability of the residual terms $\delta B_e^0(T)$ and $\delta W_{e,ij}(T)$ with respect to $M_W$ and $R_{JB}$ are explored by taking the observed standard deviations of the residuals within fixed bins of 0.5 $M_W$ units and a set of logarithmically spaced distance bins respectively (Weatherill et al. 2020). This analysis agrees with that of Kotha et al. (2022) in identifying reductions in $\tau_0$ and $\phi_0$ at higher magnitudes, while no compelling evidence was found for distance-dependent heteroskedasticity.

For the case of $\tau_0$ the magnitude-binned standard deviations are compared with a heteroskedastic model of $\tau$ proposed by Al Atik (2015), which is itself fit to a global database of ground motions from shallow crustal events and is therefore richer in observations from larger magnitudes than the ESM database. We find that for the small-to-moderate magnitudes the observed values of $\tau_0$ in the Kotha et al. (2020) homoscedastic model agrees well with the model of Al Atik (2015).

At larger magnitudes the observed $\tau_0$ is smaller than the heteroskedastic model of Al Atik (2015); however, with so few observations from large magnitude ($M \ge 6.0$) events in ESM we would consider the model of Al Atik (2015) better constrained. Therefore, we adopt the global $\tau$ model of Al Atik (2015) in place of the homoscedastic $\tau_0$ model of Kotha et al. (2020) without modification.

When considering the magnitude-dependence of observed $\phi_0$ we again compare this against a corresponding global model of $\phi_0$ proposed by Al Atik (2015). With the largest magnitude events in the ESM database being particularly well-recorded, we have greater confidence in the constraint of $\phi_0$ from the observed ESM data and thus attempt to calibrate a magnitude-dependent $\phi_0$ model.

Comparing the calibrated $\phi(M_W, T)$ values against the model of Al Atik (2015) we find good agreement at larger $M_W$ values, while for small-to-moderate magnitudes the observed $\phi_0$ are similar but slightly larger. In this case we replace the homoscedastic $\phi_0$ model of Kotha et al. (2020) with a new magnitude-dependent model calibrated upon the ESM data set. We do not adopt the Al Atik (2015) $\phi_0(M_W, T)$ model in the present case.

### 4.2.3 Site Amplification and Characterisation

The functional form of the Kotha et al. (2020) in equation 6.1 does not contain an explicit site amplification term, $f_{SITE}$. In their original paper, the authors present alternative models for $f_{SITE}$ dependent on either $V_{S30}$ or topographic slope. These are determined by regression of the $\delta S2S_S$ terms against various predictor variables, assuming a quadratic form such that $f_{SITE}(X) = g_0 + g_1 \cdot \ln\left(\frac{X}{X_{ref}}\right) + g_2 \cdot \ln\left(\frac{X}{X_{ref}}\right)^2$, where X refers to either topographic slope or measured $V_{S30}$. From their robust mixed effects regression $\delta S2S_S$ is available for more than 1,800 sites with 3 or more records, 644 of which would be considered "well constrained" by 5 or more records. From the full set of 1,800 site only 419 stations report a measured $V_{S30}$, the rest requiring a $V_{S30}$ that is inferred from slope using the method of Wald & Allen (2007).

In their estimation of $f_{SITE}(V_{S30})$ only the subset of measured observations was used, meaning that the resulting variability ($\phi_{S2S}$) represents only the site-to-site variability for measured $V_{S30}$ cases, not for
Europe as a whole. The spatial distributions of the measured and inferred V<sub>S30</sub> stations are shown in Fig. 4.6.

<figure><img src="image_48.png" /><figcaption>Fig. 4.6 Distribution of sites used in the construction of the Kotha et al. (2020) model for the case when a measured V<sub>S30</sub> available (left) and when V<sub>S30</sub> is inferred from topography (right)</figcaption></figure>

To capitalise on the model complete data set of δS2S<sub>S</sub>(T) including stations without a measured V<sub>S30</sub> case, a new f<sub>SITE</sub> regression is undertaken in which different coefficients and resulting φ<sub>S2S</sub>(T) are defined depending on whether the predictor V<sub>S30</sub> refers to a measured or an inferred site condition.

The trend in δS2S<sub>S</sub>(T) against V<sub>S30</sub> is shown in Fig. 4.7 for Sa (0.2 s) and Sa (1.0 s) considering both the measured and inferred cases. To define f<sub>SITE</sub> we use a two-segment piecewise linear function:

$$f_{\text{site}}(\theta, T) = \begin{cases} g_1(T) \ln\left(\frac{\theta}{\theta_{\text{ref}}(T)}\right) & \text{for } \theta \le \theta_c \\ g_1(T) \ln\left(\frac{\theta_c}{\theta_{\text{ref}}(T)}\right) & \text{for } \theta > \theta_c \end{cases} + \varepsilon_s \cdot \phi_{S2S}(T) \quad (\text{Eq. 6.5})$$

where θ is equal to V<sub>S30</sub> and g<sub>1</sub>(T), θ<sub>ref</sub>(T) and φ<sub>S2S</sub>(T) depend on whether the V<sub>S30</sub> is measured or inferred. θ<sub>c</sub> refers to the V<sub>S30</sub> above which f<sub>SITE</sub> is constant, which in the current implementation is fixed to 1100 m/s. The fitted model is also shown in Figure 6.9. The resulting function results in a smaller increase in amplification with decreasing V<sub>S30</sub> for the inferred V<sub>S30</sub> case than for the measured V<sub>S30</sub> case, but this is partly compensated for in the hazard calculations by the larger φ<sub>S2S</sub> term. For the measured V<sub>S30</sub> case, φ<sub>S2S</sub>(T) was found to be almost identical to that of the Kotha et al. (2020) V<sub>S30</sub>-based amplification model, and the degree of amplification comparable over the range of V<sub>S30</sub> constrained by data. This functional form does not consider nonlinearity in amplification for strong shaking on soft soil sites as the available database contains too few records of strong earthquakes on such sites to constrain a nonlinear amplification model empirically.

The dependence of the amplification model on whether the V<sub>S30</sub> corresponds to a measured or inferred site has implications for its application. The ESHM20 serves multiple purposes: i) calculation of seismic hazard on reference rock (Eurocode 8 class A, V<sub>S30</sub> 800 m/s) that may form an input into seismic resistant design codes, ii) calculation of ground shaking on the surface soil as input into seismic risk analysis including the 2020 European Seismic Risk Model (ESRM20).

In the case of the former, we define the reference rock as corresponding to a measured V<sub>S30</sub> case, and adopt the lower φ<sub>S2S</sub> value accordingly. This decision is made according to the context in which site amplification is applied in design codes.
As code amplification factors are generally calibrated with conservatism to account for uncertainty in site response, and in the current Eurocode 8 case higher amplification is required in the case in which the site has unknown $V_s$ properties, adoption of the inferred $V_{S30}$ $\phi_{S2S}$ would likely double count uncertainty already implicit in the code amplification factors.

<figure><img src="image_49.png" /><figcaption>**Fig. 4.7** Relation between $\delta S_2 S_S(T)$ and $V_{S30}$ for Sa (0.2 s) (left column) and Sa (1.0 s) (right column) for the measured (top row) and inferred $V_{S30}$ case (bottom row). Blue lines and grey shaded region indicate a non-parametric LOWESS regression fit, while red lines indicate the fit of the two-segment linear $f_{SITE}$ model.</figcaption></figure>

For the seismic risk case, including ESRM20, the corresponding measured or inferred $V_{S30}$ can be chosen depending on whether the sites in question take their properties from direct measurement, or detailed local scale microzonation, or whether they are inferred from regional scale proxies, such as topographic slope.

As the ESRM20 requires characterisation of the soil properties across all of Europe, the inferred $V_{S30}$ case is dominant. Though not described here, the ESRM20 actually adopts a form of equation 6.5 that is dependent on both slope *and* local geology, thus refining the degree of amplification depending on the geological environment. $\phi_{S2S}(T)$ values are found to be very similar regardless of whether topographically inferred $V_{S30}$ is used or slope directly (with or without local geology).

Further details on the construction, application and seismic hazard and risk implications of the site amplification model can be found in Weatherill et al. (2022, in preparation).

## 4.3 Scaled Backbone GMM Logic Tree for the Stable Craton

### 4.3.1 Defining and delineating the cratonic region of Europe

In the previous section it has been indicated that the *default* shallow crustal seismicity GMM logic tree is applied in regions for which we have insufficient data to calibrate $\delta L2L_l$ and/or $\delta c_{3,k}$. A fundamental assumption that underpins this, however, is that the seismological characteristics of the regions for which we lack data are within the range implied by the full region-to-region distributions determined using the ESM data (Lanzano et al 2019).
As data is lacking for much of northern Europe, we cannot necessarily say for certain to what extent this
is true. We therefore consulted other geophysical data sets whose properties may be relevant for
ground motion characterization, and from these sought to identify where regional differences exist.

In Fig. 4.8 we show the spatial patterns of 1-second Rayleigh wave attenuation quality factor ($Q_{LG}$) (Mitchell *et al.*, 2008), upper mantle shearwave velocity anomaly at 120 km depth ($\delta V_s$) (Mooney et al., 2012), heat flow (Lucazeau, 2019) and Moho depth (Grad et al. 2009; Szwillus et al. 2019).

From these data sets we observe a clear difference in geophysical properties when we contrast southern
and western Europe against northeastern Europe. Specifically, in northeastern Europe we observe
higher $Q_{LG}$ (implying slower attenuation), higher mantle shearwave anomaly, lower heatflow and deeper
Moho. The transition between the two environments is generally clear and corresponds to the Trans-
European Suture zone (TESZ).

<figure><img src="image_50.png" /><figcaption>Fig. 4.8 Geophysical data sets across Europe: $Q_{LG}$ (top left), $\delta V_s$ (top right), Heatflow (bottom left) and Moho depth (bottom right)</figcaption></figure>

As the observed strong motion data in the ESM database fall entirely on the southern and western size
of the TESZ, it is not possible to assume with certainty that the seismological characteristics of
northeastern Europe are represented within the region-to-region variability. Instead, we find that a
better analogue can be found in eastern North America.

This is not necessarily a novel insight, as previous seismic hazard models for northern Europe, including
the ESHM13, have adopted GMMs from the central and eastern United States (CEUS). We make use of
this analogue in order to leverage upon recent developments in ground motion modelling within the
CEUS as a means of potentially regionalizing our shallow crustal scaled backbone logic tree approach
further. We summarise this in the following, but full details describing the development of the GMM
model and scaled backbone GMM logic tree for northeastern Europe can be found in Weatherill &
Cotton (2020)
### 4.3.2 Calibrating the scaled backbone GMM for median ground motions

The Next Generation Attenuation East (NGA East) project represents a significant step forward in strong ground motion modelling for the central and eastern United States (Goulet et al. 2018). Among the wide range of outputs and publications from the project we find a suite of 20 new individual median GMMs for very hard rock ($V_s$ 3000 m/s) conditions, the characterization of a full GMM logic tree to represent epistemic uncertainty, a site amplification model (Stewart et al. 2020; Hashash et al. 2020) and an aleatory uncertainty model (Al Atik, 2015). In the absence of observed ground motion data for northeastern Europe, we utilize these products to undertake a pseudo-regionalisation of the scaled backbone ground motion model.

We begin with the 20 individual ground motion models for CEUS, which are developed using a mixture of simulated, empirical or hybrid ground motion modelling methods with different assumptions about the seismological properties of the crust.

We also add the model of Pezeshk et al. (2011), which precedes the NGA East project but adopts many of the same characteristics. Each model is assumed equally weighted. From the 21 GMMs we generate expected median ground motions (on very hard rock) for a broad range of scenario magnitudes ($4.0 \le M_W \le 8.0$) and distances ($1.0 \le R_{RUP} (km) \le 500$ km). The epistemic uncertainty on the median ground motions is then described by a log normal distribution:

$$P(\ln Y(M_W, R_{RUP}, T)) = \mathcal{N}(\mu(M_W, R_{RUP}, T), \sigma_{\mu}(M_W, R_{RUP}, T)) \quad (\text{Eq. 6.6})$$

The distribution is scenario-dependent, where $\mu$ is the median expected ground motion and $\sigma_{\mu}$ represents the epistemic uncertainty in the median.

Though this distribution could be used in this non-parametric form to capture the ground motion, we instead use this distribution of scenario ground motions to fit a parametric GMM adopting the same functional form as that of Kotha et al. (2020) (equations 6.1 – 6.4). Several modifications are made, however, the preferred distance metric of the NGA East GMMs is closest distance to rupture ($R_{RUP}$). As such $R_{ref}$ is now fixed to 1 km and $h_d$ to 5 km (no longer hypocentral depth dependent):

$$\ln Y_{\text{median}}(M_W, R_{RUP}, T) = e_1 + f_M(M_W, T) + f_{R,g}(M_W, R_{RUP}, T) + f_{R,a}(R_{RUP}, T) + \varepsilon_{\mu} \cdot \sigma_{\mu}(T) \quad (\text{Eq. 6.7})$$

This model is fit using robust least squares regression, with no mixed effects. The standard deviation ($\sigma_{\mu}$) therefore represents the model-to-model variability. As with the shallow crustal GMM logic tree, we use the N-branch discrete approximation to a Gaussian distribution by Miller & Rice (1983) to then map the epistemic uncertainty distribution ($\mathcal{N}(\mu, \sigma_{\mu})$) into a scaled backbone GMM logic tree.

The response spectra from the suite of NGA East models, the default (or general) shallow crustal seismicity logic tree, the proposed parametric ground motion model logic tree and a recent GMM developed for nuclear power plants in Finland (Fülöp et al. 2020) are compared for four scenarios in Fig. 4.9.

The parametric GMM and its epistemic uncertainty generally capture the range of median ground motions from the NGA East model well and are in good agreement with the GMM of Fülöp et al. (2020). In contrast with the default GMM logic tree, the ground motions are significantly higher for short periods ($T < 0.3$ s) though the centre and range of the two logic trees converge for $T \ge 1.0$ s.
<figure><img src="image_51.png" /><figcaption>**Fig. 4.9** Comparison of the GMM logic tree of the craton GMM model (N = 3 branches) (black) with the NGA East GMM suite (coral), the *default* crustal seismicity logic tree (grey) and the Fülöp et al. (2020) GMM (blue) for a reference rock site of V<sub>S30</sub> 800 m/s.</figcaption></figure>

4.3.3 Site Amplification

As the parametric form of the model is fit to the NGA East models, the ground motion values are initially
defined for a very hard rock (V<sub>S</sub> 3000 m/s) condition. To determine the ground motion on the reference
V<sub>S30</sub> 800 m/s rock condition, and to other site conditions, the CEUS site amplification models of Stewart
et al. (2020) and Hashash et al. (2020) are integrated into the model. These amplification models have
three components: a very hard rock to reference rock amplification factor (F<sub>760</sub>(T)), a linear
amplification factor (f<sub>lin</sub>(V<sub>S30</sub>, T)) and a nonlinear amplification factor (f<sub>nl</sub>(V<sub>S30</sub>, PGA<sub>760</sub>, T)), each of
which has an epistemic uncertainty that when combined forms an additional distribution that we define
here as σ<sub>μ,S</sub>. One critical difference between the present model and similar GMMs that have been
applied previously in Europe is that for sites on very stiff soil or rock the F<sub>760</sub> and f<sub>lin</sub> factors assume a
seismic velocity profile of the site with a strong impedance contrast as depth, typical of what might be
expected in regions of low tectonic deformation and glacially formed geomorphology. This impedance
contrast results in a stronger amplification for high frequency motion than that of previous models that
have assumed a more gradational velocity profile. The new GMM therefore predicts higher short-period
motion than the previous ESHM13 selection; however, we believe this is appropriate to the region in
question.

4.3.4 Aleatory variability

In keeping with the aleatory variability models assumed by both the NGA East ground motion project
*and* the default shallow crustal seismicity logic tree, we adopt the same heteroskedastic aleatory
uncertainty model from Al Atik (2015). In this case we use the Al Atik (2015) “global” models (fit to the
NGA West 2 database) with magnitude-dependent τ and magnitude-dependent φ<sub>0</sub>. Site-to-site
variability ($\phi_{S2S}$) is taken from the study of Stewart et al. (2019), which is used for the ergodic seismic hazard assessment in the 2018 United Stated National Seismic Hazard Model. No epistemic uncertainty on the aleatory variability term is considered in the present case.

### 4.3.5 Complete Logic Tree for GMM in the Stable Cratonic Region of NE Europe

To define the complete GMM logic tree for this region (shown in Figure 6.12) we return to our initial premise that the scaling and attenuation of ground motions in NE Europe are fundamentally different from those of southern and western Europe. Without strong motion data within the ESM to verify this assumption, we retain in the logic tree both the possibility that the ground motion characteristics are fundamentally different, and the possibility they are not. The logic tree if first split into two branches to accommodate both hypotheses. Initially we consider the possibility that some branches of the default crustal seismicity logic tree may be appropriate in NE Europe. This is assigned a weight of 0.2. We do not, however, consider the low stress region branches or the fast attenuation branches to be feasible and thus eliminate them from the tree. The branches are re-weighted to assign $W_{11} = W_{111} = W_{121} = 0.667$ weight to the high stress ($\epsilon_{11} = 1.732$) and to the slow attenuation branches with ($\epsilon_{111} = \epsilon_{121} = 1.732$), and just $W_{12} = 0.333$ weight to each of the “average” stress region and attenuation region branches ($W_{112} = W_{122} = 0.333$) respectively ($\epsilon_{12} = \epsilon_{112} = \epsilon_{122} = 0$).

<figure><img src="image_52.png" /><figcaption>**Fig. 4.10** Complete GMM logic tree for the stable cratonic region of northeastern Europe</figcaption></figure>

On the other side of the logic tree, we have the assumption that the ground motion characteristics of NE Europe are fundamentally different from those represented by any branch of the default crustal seismicity GMM logic tree, and therefore apply the parameteric craton GMM and its epistemic uncertainty. Two branching levels are considered to reflect both the epistemic uncertainty in the median ground motion on very hard rock ($\sigma_\mu$) and the epistemic uncertainty in the site amplification factor ($\sigma_{\mu,S}$). As the former is the dominant epistemic uncertainty, we apply a 5-branch discrete approximation to the Gaussian distribution, while for the latter we apply a 3-branch discrete approximation (see Table 6.1 for the $\epsilon$ values and weights).
## 4.4 Subduction and Deep Seismicity

The ESM database is rich in records from shallow crustal seismicity, and this has motivated the development of the regionalized scaled backbone ground motion model approach thus far. Some of most active seismic regions in Europe, however, and indeed some of the areas of highest seismic hazard come from seismicity that does not take place in the shallow crust but rather from Europe's subduction zones (Hellenic, Calabrian, Cypriot and Gibraltar) and from the Vrancea deep seismic zone in Romania. The ground motion scaling and attenuation characteristics from earthquakes in these regions can differ significantly from those of shallow crustal seismicity, prompting the use of ground motion models specifically calibrated for subduction earthquakes. Fortunately, the ESM database contains hundreds of records from earthquakes originating in the Hellenic and Calabrian arcs, and from the Vrancea deep seismic zone (DSZ).

### 4.4.1 Identifying “non-shallow crustal earthquakes”

In its initial form the ESM database does not specifically indicate the tectonic region of the earthquake, and therefore events are not classed as shallow, subduction, etc. The first step is to sort the records from subduction regions and the Vrancea DSZ out from the other shallow crustal events. In the case of the Vrancea DSZ this is done manually, first by identifying the subset of earthquakes in the Vrancea region whose hypocentral depths are between 60 and 250 km, then verifying these by inspection to check for misclassifications or geographical outliers from the observed deep seismicity. This process yields 938 records from 64 deep source events in the Vrancea region.

For the subduction regions the identification process is more complex as the earthquakes cannot be so easily discriminated by the simple depth criteria. Instead, the southern Hellenic region and Calabria experience earthquakes in the upper 50 km that may originate from shallow seismicity in the over-riding crust, subduction interface seismicity or even seismicity within the subducting slab or the outer rise.

Given the observed uncertainty in both the depth estimates of the earthquakes and the location of the crust and subducting slab, a fuzzy classification system was designed to classify earthquakes in terms of degree of membership of each tectonic region type given the event's depth, proximity to the subduction surface and, if available, the focal mechanism.

A comprehensive description of the fuzzy classification algorithm is given in the accompanying report of Weatherill et al. (2021) and has also been applied to similar problems of earthquake classification in subduction regions by Ghasemi et al. (2020). Setting a *degree of membership* threshold of 0.6 (where 0 is unambiguously not a classified according to a given tectonic region type and 1.0 is unambiguously classified as so) we apply the fuzzy classification system to categorise events as one of: active shallow crust, subduction interface, subduction in-slab, upper mantle, outer rise and non-subduction (i.e., deep events not seemingly associated with the subduction zone).

From this classification we identify 681 records from 164 subduction interface earthquakes and 704 records 121 from subduction in-slab events. The distribution of records for the three “non-shallow” seismicity types (subduction interface, subduction in-slab and Vrancea DSZ) are shown in Fig. 4.11.
<figure><img src="image_53.png" /><figcaption>**Fig. 4.11** Magnitude and distance composition of the ground motion databases for the subduction interface, subduction in-slab and Vrancea DSZ regions of Europe.</figcaption></figure>

### 4.4.2 Selecting and modifying a backbone deep/subduction seismicity GMM

In contrast to both the case of shallow seismicity (for which a substantial amount of data was available in ESM) and the cratonic region (for which no data were available), the subduction/deep records provide us with enough data to make some inferences about the ground motion characteristics of the region, but not enough to fit a full GMM. In particular, Figure 6.13 shows that while the number of records is quite high overall, the magnitude range is limited to $M \le 6.0$ to 6.5 for the subduction events while even the Vrancea database has only a limited number of records for $M \ge 6.0$ earthquakes.

As larger magnitudes in the range $6.5 \le M_w \le 8.5$ are the most relevant for seismic hazard we instead aim to construct a backbone model by first using the available data to identify the most appropriate candidate GMM from the suite of recent global subduction models and then attempting to calibrate this model to the attenuation characteristics of the southern Mediterranean and Vrancea regions.

An initial pre-selection process is applied considering the required range of spectral periods and earthquake scenarios for the ESHM20. At the time of undertaking this analysis the Next Generation Attenuation Subduction (NGA-Sub) ground motion models (Parker et al. 2020; Kuehn et al. 2020; Gülerçe & Abrahamson, 2020) were not available. Though the NGA-Sub models do introduce some calibrations for different subduction zones around the globe, no data from the eastern Mediterranean were included and thus no calibration is currently available.

The pre-selected models are compared against the observed data in each of the three region types quantitatively using multivariate log-likelihood scoring (Mak et al., 2017). From this pre-selection process we find that the Abrahamson et al. (2016) “BC Hydro” model is the most suitable candidate for the subduction regions, meeting the project requirements in terms of scenario and period range while also giving low log-likelihood scores across the full range of spectral periods in each of the three regions. For the Vrancea DSZ the in-slab version of the Abrahamson et al. (2016) model is also found to be the best suited.

With the candidate model identified, we then use the available strong motion data to explore the trends in the between-event and within-event residuals with respect to this model. A lack of data from large magnitude events, combined with limited a limited number of sites with measured $V_{S30}$ terms, means that we do not attempt to re-calibrate either magnitude or linear site scaling terms of the GMM.
However, with the large volume of data covering the 60 to 200 km distance range we use the observed trends in within-event residuals to calibrate the anelastic attenuation term of the model to better fit the local data.

We do this by first using the site amplification component of the GMM to scale the ground motions to reference rock. Then within-event residuals are calculated with respect to an adaptation of the model in which the anelastic term ($f_{R,A}(R_{RUP}) = \theta_6 \cdot R_{RUP}$) is set to 0. The locally calibrated $\theta_6$ is retrieved for each period from the slope of a linear fit to the within-event residuals.

The observed period-dependent trend in $\theta_6$ for both the interface and in-slab models is shown in Fig. 4.12, where it is compared against similar anelastic attenuation scaling coefficients from other relevant subduction GMMs as well as from regional estimates of this same parameter from Abrahamson et al. (2018), a preliminary update of the 2016 model using NGA Sub data.

<figure><img src="image_54.png" /><figcaption>**Fig. 4.12** Variation in anelastic attenuation term ($\theta_6$) with period for the fitted cases (red and black dots for in-slab and interface events respectively) and selected subduction GMMs and local regionalisations, where “Global”, “Cascadia”, “Central America”, “Japan”, “New Zealand”, “South America” and “Taiwan” correspond to the regions for which $\theta_6$ terms have been calibrated by Abrahamson et al. (2018)</figcaption></figure>

### 4.4.3 Quantifying the epistemic uncertainty and constructing the logic tree

Fig. 4.12, demonstrates that there is a significant diverge in anelastic attenuation coefficients from different regions of the globe, particularly at short periods where they diverge more significantly between two groups. Recognising this variability, and considering the potential uncertainties in our present estimates of $\theta_6$ for the regions of interest, we define an average model (solid red and black lines in Figure 6.14) and an upper and lower bound (dashed red and black lines) based around the global variability in this term, which we define as a period-independent ±0.0015.

For consistency with the scaled backbone approaches adopted elsewhere, we consider these upper and lower bounds as notionally equivalent to the 95th and 5th percentile of a Gaussian distribution, and hence we define a three-branch distribution in the logic tree of $\theta'_6 - 0.0015, \theta'_6$ and $\theta'_6 + 0.0015$ with weights of 0.167, 0.666 and 0.167 accordingly for the subduction zones, and weights of 0.2, 0.6 and 0.2 for Vrancea where we believe slightly higher weights on the upper and lower branches may be appropriate.
The calibration of $\theta_6$ may account for epistemic uncertainty in the anelastic attenuation, but its influence is mostly apparent only at longer distances. This does not address epistemic uncertainties due to large magnitude scaling and/or regional differences in the stress parameter term of the model.

For this we take advantage of an epistemic uncertainty distribution for median ground motions ($\sigma_\mu$) proposed by Abrahamson et al. (2018), who looked at both the region-to-region variability in stress parameter ($\theta_1$). Their model presents a $\sigma_\mu$ term that is constant with period for interface models and increases slightly at shorter period for in-slab events.

As their analysis is based on the stress-parameter, this term does not account for region-to-region variability in attenuation, meaning that it can be applied as an additional set of branches on top of the anelastic attenuation uncertainty model with a minimal risk of double counting.

An example of the resulting ground motions for a three-branch approximation of the distribution of $\sigma_\mu$ is shown in Figure 6.15 for a M 7.5 event at 25 km depth (for the interface) and 60 km depth (in-slab) assuming a site condition of 800 m/s and a forearc path. The distribution of $\sigma_\mu$ is shown in **Fig. 4.13**.

<figure><img src="image_55.png" /><figcaption>**Fig. 4.13** Ground motion attenuation for the ESHM20 subduction GMMs for a three-branch distribution of $\sigma_\mu$, where $\varepsilon$ corresponds to $\pm 1.732$ (according to the 3-point distribution shown in Table 6.1). Mw = 7.5, $V_{S30}$ = 800 m/s and hypocentral depth is 60 km for the in-slab model</figcaption></figure>

Combining the epistemic uncertainties on the stress-parameter (characterized here by $\sigma_\mu$) and the attenuation, we define a complete ground motion model logic tree for application to the European subduction interface and in-slab regions and Vrancea DSZ. The Abrahamson et al. (2016) in-slab model is used as the backbone for both the in-slab and Vrancea DSZ regions.

As the epistemic uncertainty with the largest influence on hazard, we map $\sigma_\mu$ into a 5-branch discrete approximation of the Gaussian distribution according to Miller & Rice (1983). The scaling factors and corresponding weights for the anelastic attenuation adjustments are stated previously. The complete logic trees for each of the three environments (subduction interface, subduction in-slab and Vrancea DSZ) are shown in **Fig. 4.14**.
#### 4.4.4 Modifying the forearc/backarc scaling term

The Abrahamson et al. (2016) GMM incorporates a forearc/backarc scaling term, one in which attenuation with distance is faster for sites in the backarc of the subduction zone.

Trends in within-event residuals were explored to identify if this forearc/backarc scaling is present in the observed data, and although data in the backarc are sparse there was sufficient evidence to suggest such a distinction is present. However, with so few data in the backarc we could not re-calibrate the coefficients of the scaling term within the GMM.

To define the boundary between the forearc and backarc for the Hellenic and Calabrian subduction zones we first consider the Slab 2.0 representation of the zone geometry. Then using the locations of predominant volcanic activity as the defining boundary we draw the forearc/backarc front based on the relative positions of the volcanic front with respect to the interface.

For Vrancea we follow the example of Vacareanu et al. (2015) in defining the forearc/backarc margin along the watershed of the Carpathian mountain belt. For the Gibraltar arc we cannot identify a clear front, and in the present case all sites are assumed to represent forearc sites. The lines of the fronts can be seen in Fig. 4.15.

<figure><img src="image_56.png" /><figcaption>**Fig. 4.14** Complete scaled backbone logic tree formulations for the subduction interface, in-slab and Vrancea DSZ regions.</figcaption></figure>

Though the analysis of the ground motion residuals does support the use of a forearc/backarc scaling term, in the Abrahamson et al. (2016) GMM this is treated as a binary term, $F_{FABA}$, taking 0 or 1 depending on whether the site is located in the forearc or backarc respectively. As the impact on the attenuation is significant, with lower ground motions in the backarc, this can result in a “cliff-edge” drop in seismic hazard across the front.

This was particularly pronounced in Romania where the Vrancea DSZ is the dominant seismic source over several hundreds of kilometres. To modulate this effect, we modify the forearc/backarc definition to describe $F_{FABA}$ as a continuous variable in the range $0 \le F_{FABA} \le 1$, dependent on the distance to the forearc/backarc front, thus effecting a smoother transition from one attenuation regime to the other across the front. To accommodate this a new site term is introduced, which defines the distance of the
site with respect to the forearc/backarc front (or volcanic front), $x_{vf}$, and this parameter is calculated for all three major subduction and deep seismicity regions. $F_{FABA}$ is related to $x_{vf}$ through a simple function that describes the smooth transition from “pure” forearc to “pure” backarc attenuation over a set distance range. Several functional forms were tested but the preferred relation is an S-function such that:

$$F_{FABA}(x_{vf}, a, b) = - \begin{cases} 0 & \text{for } x_{vf} \le a \\ 2 \cdot \left(\frac{x_{vf}-a}{b-a}\right)^2 & \text{for } a < x_{vf} \le 0.5 \cdot (a+b) \\ 1 - 2 \cdot \left(\frac{x_{vf}-a}{b-a}\right)^2 & \text{for } 0.5 \cdot (a+b) < x_{vf} < b \\ 1 & \text{for } x_{vf} \ge b \end{cases} \quad (\text{Eq. 6.8})$$

where a and b are -150 km and 150 km respectively. The spatial pattern of $x_{vf}$ is shown in Figure 6.17, along with an illustration of the influence of this term on the attenuation of ground motion from a scenario Mw 7.7 earthquake (with the volcanic front placed at 250 km from the source).

<figure><img src="image_57.png" /><figcaption>Fig. 4.15 Location of forearc/backarc front and $x_{vf}$ distance in eastern and southern Europe (left), and impact of smoothed functions for $F_{FABA}$ on attenuation for a scenario event of Mw 7.7 on rock with the volcanic front 250 km from the source</figcaption></figure>

A small region of conflicting overlap can be seen at the North Macedonia/Bulgaria border region, which falls into the backarc region with respect to the Hellenic arc but also the forearc region with respect to the Vrancea DSZ. In this narrow region the sites are classified as forearc in order to ensure that the contribution to the seismic hazard from Vrancea is dominant over the very distant Hellenic sources.

No epistemic uncertainty on the functional form and distance range of the $F_{FABA}$ is modelled here, though this is accommodated in the model implementation and thus may be considered in future projects.

## 4.5 Special cases and adjustments

The three main tectonic region types that have been considered so far, i.e., shallow crustal seismicity, stable craton and subduction/deep seismicity, cover most of the seismic sources in Europe.

With the adoption of the regionalization features in the shallow crustal seismicity model GMM we are also able to capture a further level of detail than has been assumed in the previous classifications adopted by ESHM13 and some existing national models.
As this approach has been largely data driven, however, the calibrations themselves can be readily quantified and adopted in regions with sufficient data *within the available ESM dataset*. There are areas, however, where ground motion data are available in a limited capacity, but were not included within the ESM data set and therefore were not part of the main regression process.

As the ground motion models and its epistemic uncertainties *should* represent the best state-of-knowledge at the time of construction we consider two *special cases* of shallow seismicity where we adapt the logic tree on the not on the basis of the data used for the regression but rather from other data, information and interpretations as well as feedback from third parties.

### 4.5.1 North-western Europe

The general strategy adopted within the shallow default crustal logic tree is to assign regions the full region-to-region variability in stress region scaling and residual attenuation where there are insufficient data to constrain a local calibration. This would be the case much of France (excluding the Pyrenees, Alpine and Rhine regions), Norway, the United Kingdom, Spain and Portugal. Subsequent analysis by Kotha (2020) explored in detail a database of strong and weak motion records from the French seismological and Geodetic Network (RESIF), which contains a large number of stations and events that were not available within the ESM database. In the majority of cases these events from RESIF were smaller magnitude ($2.5 \le M \le 4.5$) than those found in ESM. For these reasons, and accounting for differences in processing that may introduce potential biases between the two data sets, we do not include the RESIF data into the initial regressions.

The analysis by Kotha (2020) was able to define $\delta c_{3,k}$ terms for additional regions in France and parts of Spain, while also providing a general comparison point for overlapping regions in southeastern France and Italy. The regional variation in $\delta c_{3,k}$ for a period close to Sa (0.1 s) is shown in **Fig. 4.16**.

From this analysis the general spatial trend of slower residual attenuation to the north and west of the Alps contrasted against progressively faster attenuation in Italy is persistent. What the analysis on the RESIF data does suggest, however, is that the slow residual attenuation seen in the Pyrenees is potentially persistent throughout much of France west of the Alps.

However, it is important to note that the distribution of stations for in the RESIF data set for this region is skewed toward southwestern France and the Pyrenees, and to the northwest in Brittany and Pays de la Loire.

For central and northern France few stations exist in the RESIF database, which leads us to conclude that we would retain the default general crustal logic tree here, conceding that this region would most likely be a transition zone between the faster residual attenuation in the southeast and the particularly slow residual attenuation in the north and southwest.