# A regionally-adaptable ground-motion model for shallow crustal earthquakes in Europe

Sreeram Reddy Kotha, Graeme Weatherill, Dino Bindi, Fabrice Cotton

► To cite this version:

Sreeram Reddy Kotha, Graeme Weatherill, Dino Bindi, Fabrice Cotton. A regionally-adaptable ground-motion model for shallow crustal earthquakes in Europe. Bulletin of Earthquake Engineering, 2020, 18 (9), pp.4091-4125. 10.1007/s10518-020-00869-1. hal-04502061

HAL Id: hal-04502061

https://hal.science/hal-04502061v1

Submitted on 14 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A regionally-adaptable ground-motion model for shallow crustal earthquakes in Europe

Sreeram Reddy Kotha<fnref n="1" />,<fnref n="2" /> ID · Graeme Weatherill<fnref n="1" /> · Dino Bindi<fnref n="1" /> · Fabrice Cotton<fnref n="1" />,<fnref n="3" />

Received: 30 January 2020 / Accepted: 12 May 2020 / Published online: 21 May 2020
© The Author(s) 2020

## Abstract

To complement the new European Strong-Motion dataset and the ongoing efforts to update the seismic hazard and risk assessment of Europe and Mediterranean regions, we propose a new regionally adaptable ground-motion model (GMM). We present here the GMM capable of predicting the 5% damped RotD50 of PGA, PGV, and SA(T = 0.01 – 8 s) from shallow crustal earthquakes of $3 \le M_W \le 7.4$ occurring $0 < R_{JB} \le 545$ km away from sites with $90 \le V_{s30} \le 3000 \text{ m s}^{-1}$ or $0.001 \le \text{slope} \le 1 \text{ m m}^{-1}$. The extended applicability derived from thousands of new recordings, however, comes with an apparent increase in the aleatory variability ($\sigma$). Firstly, anticipating contaminations and peculiarities in the dataset, we employed robust mixed-effect regressions to down weigh only, and not eliminate entirely, the influence of outliers on the GMM median and $\sigma$. Secondly, we regionalised the attenuating path and localised the earthquake sources using the most recent models, to quantify region-specific anelastic attenuation and locality-specific earthquake characteristics as random-effects, respectively. Thirdly, using the mixed-effect variance-covariance structure, the GMM can be adapted to new regions, localities, and sites with specific datasets. Consequently, the $\sigma$ is curtailed to a 7% increase at $T<0.3$ s, and a substantial 15% decrease at $T \ge 0.3$ s, compared to the RESORCE based partially non-ergodic GMM. We provide the 46 attenuating region-, 56 earthquake localities-, and 1829 site-specific adjustments, demonstrate their usage, and present their robustness through a 10-fold cross-validation exercise.

**Keywords** Ground-motion model · Response spectra · Robust mixed-effects regression · Regionally adaptable · Seismic hazard and risk · Europe

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

<footnote>Sreeram Reddy Kotha
sreeram-reddy.kotha@univ-grenoble-alpes.fr</footnote>

<footnote n="1">¹ Helmholtz Centre Potsdam, GFZ German Research Centre for Geosciences, 14467 Potsdam, Germany</footnote>

<footnote n="2">² Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, IFSTTAR, ISTerre, 38000 Grenoble, France</footnote>

<footnote n="3">³ Institute of Geosciences, University of Potsdam, Potsdam, Germany</footnote>
# 1 Introduction

Ground-Motion Models (GMMs) characterise the random distribution of ground-motions at site for a combination of earthquake source, wave travel-path, and the affected site's geological properties. Typically, GMMs are regressed over a compendium of strong ground-motion recordings collected from several earthquakes recorded at multiple sites scattered across a variety of geographical regions. The necessity of compiling such large datasets is to expand the range of magnitude and distance, and diversity of site-types, in order to derive a GMM capable of predicting realistic ground-motions for prospective seismic risk scenarios, e.g. large magnitudes at short distances from a reference rock site. NGA-West2 (Ancheta et al. 2014) is one such dataset compiled of ground-motion observations recorded around the globe—primarily from Western US, and in smaller fractions from Alaska, China, Italy, Japan, Taiwan, Turkey, etc. Several GMMs have been derived from this dataset for application in probabilistic seismic hazard (PSHA) and risk assessments. Given the clear tectonic and geological diversity of the data, possible regional differences in observed ground-motions needed to be quantified (Douglas 2004). Four of the NGA-West2 GMMs (Abrahamson et al. 2014; Boore et al. 2014; Campbell and Bozorgnia 2014; Chiou and Youngs 2014) accounted for regional differences in ground-motions through region-specific regression coefficients. Through region-specific adjustments, these GMMs were able to capture and predict, for example, the faster attenuation of short-period ground-motions with distance in Japan compared to Western US, and the relatively slower attenuation in China. Earlier GMMs were not capable of such predictions simply owing to the lack of sufficient data from individual regions to quantify the differences.

In case of Europe and the Middle-East, RESORCE (Akkar et al. 2014b) is one such dataset compiled of data from Italy, Turkey, Greece, and other active seismic regions in pan-Europe. Using mixed-effects regression algorithms (Abrahamson and Youngs 1992), a few GMMs (Akkar et al. 2014a; Bindi et al. 2014) were derived and used in regional (Grünthal et al. 2018) and continental PSHA (Giardini et al. 2018; Woessner et al. 2015). However, these GMMs were not regionalised despite known geological differences between Italy and Turkey. In fact, even the regionalised GMMs (Bora et al. 2017; Kale et al. 2015; Kotha et al. 2016; Kuehn and Scherbaum 2016) limited the distinction to geopolitical boundaries, yet the geological diversity within these regions is far more complex. In essence, quantification of regional differences requires, first, a regionalisation scheme, and then, sufficient data from each region. With the arrival of the new European Strong-Motion (ESM) dataset (Bindi et al. 2018b; Lanzano et al. 2019a) and regionalisation models, and ongoing efforts to update the pan-European PSHA, a revision of the regionalised pan-European GMMs is proposed.

In this study, we present an upgrade of the RESORCE dataset based region- and site-specific GMM (Kotha et al. 2016, 2017) for shallow crustal earthquakes using the new ESM dataset. We often see GMMs evolve with progressively larger datasets, and supersede their older versions in terms of applicability (Bommer et al. 2010). However, with increasing data and complex parametrization, a reduction in the apparent aleatory variability ($\sigma$) of GMMs could not be achieved (Douglas and Edwards 2016; Strasser et al. 2009). Of course, increasing amount of data comes with increasing the spatiotemporal diversity of ground-motion observations, and thus an increasing $\sigma$. One approach would be to introduce new predictor variables into the GMM, but then, uncertainty or unavailability of predictor values then becomes an issue during application (Bindi et al. 2017; Kuehn and Abrahamson 2017). With this in mind, this revision of the Kotha et al. (2016) GMM will

<figure><img src="image_3.png" /></figure>
attempt to regionalise and refine the aleatory variability, while maintaining its original par-
ametrization. In addition, assuming a possible contamination of data and deviation from
the assumption of log-normality, instead of the usual ordinary least-square estimates of
GMM median and variances we compute their robust counterparts; while flagging outlier
events, stations, and records in the dataset. Table 1 summarizes the additional features we
introduced in this new GMM (abbreviated K20) with respect to the recent pan-European
GMMs: AK14 (Akkar et al. 2014a), B14 (Bindi et al. 2014), K16 (Kotha et al. 2016), L19
(Lanzano et al. 2019b).

## 2 Ground-motion data and selection criteria

Figure 1 compares the data distribution between the RESORCE and ESM datasets. The increase in amount of data between 2014 and 2018 for GMM development is dramatic. While the K16 GMM was regressed over 1251 records, the proposed revision (K20) derives from 18,222 records. One striking feature of the ESM dataset is the number of stations with $\ge 3$ ground-motion recordings. RESORCE dataset had about 150 stations with $\ge 3$ records, while ESM has 1077. This increase is highly sought in empirical site-specific GMM, PSHA, and seismic risk applications (Faccioli et al. 2015; Kohrangi et al. 2020; Kotha et al. 2017; Rodriguez-Marek et al. 2013).

Another feature is the remarkable increase in the number of small earthquakes with $M_W \le 4.5$. This could imply an increase in $\sigma$ of the revised GMM over K16, similar to the reported increase for NGA-West2 GMMs with respect to their superseded NGA-West counterparts. Nevertheless, ground-motions from frequent small-moderate sized earthquakes drive the hazard in low-moderate seismicity regions of Europe, such as France (Drouet et al. 2020) and Germany (Grünthal et al. 2018). It is necessary that the GMM is well behaved in these small-moderate magnitude ranges. Moreover, if the site-specific terms $\delta S2S_s$ were to be estimated with low uncertainty, data from several small-moderate sized earthquakes is indispensable.

The distance (Joyner-Boore metric, $R_{JB}$) range and density of data is also superior to RESORCE. The region-specific anelastic attenuation terms of K16, NGA-West2 and other GMMs (Sedaghati and Pezeshk 2017) were estimated from records at $R_{JB} \ge 80$ km. This revised GMM aims to refine the regionalisation of K16, and therefore, such increase in far-source data is quite useful. In addition, advanced studies on spatial and temporal variability of attenuation are now possible (Bindi et al. 2018a; Dawood and Rodriguez-Marek 2013; Kotha et al. 2019; Landwehr et al. 2016; Piña-Valdés et al. 2018; Sahakian et al. 2019).

The data visualized in Fig. 1 are from shallow crustal earthquakes in the ESM dataset. The full dataset contains ground-motions from other tectonic regimes as well, such as subduction interface, subduction in-slab, Vrancea, etc. To filter out these and other records not suitable for a shallow crustal GMM development, we adopt the following selection criteria:

1. To keep only the shallow crustal earthquakes, we select only those events classified as non-subduction events by Weatherill et al. (2020d). The selection removes in-slab, interface, outer-rise, and upper-mantle events from the regression dataset. The resulting 927 events with event depth $0 < D \le 39$ km are located in regions with $14 \le \text{Moho depth} \le 49$ km, as per the Moho map of Grad et al. (2009)

2. Only those events with $\ge 3$ records in the dataset are used in regression. Also, wherever available, the harmonised $M_W$ estimates from the European Mediterranean Earthquake
<table><caption>**Table 1** Comparison of features between the recent pan-European GMMs: AK14 (Akkar et al. 2014a), B14 (Bindi et al. 2014), K16 (Kotha et al. 2016), L19 (Lanzano et al. 2019b), and K20 is the new model</caption><thead><tr><th>GMM feature</th><th>AK14</th><th>B14</th><th>K16</th><th>L19</th><th>K20</th></tr></thead><tbody><tr><td>Event-to-event variability</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Event locality-to-locality variability</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Site-to-site variability</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>V<sub>530</sub> based linear site-response</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>V<sub>530</sub> based nonlinear site-response</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Topographic slope based site-response</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Depth dependent attenuation</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Regionalised anelastic attenuation</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Scalar style-of-faulting adjustments</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Anisotropic radiation pattern</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Median epistemic uncertainty</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Heteroscedastic variances</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Robust variance estimates</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Outlier detection</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

Shaded cell if black indicates a feature included, white for a feature not included, and grey for a feature tested but not included in the final GMM
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Data distribution comparing RESORCE (yellow) and ESM (black) datasets</figcaption></figure>

Catalogue [EMEC] (Grünthal and Wahlström 2012) are preferred over the ESM default values

3. We keep all sites in the dataset irrespective of whether their $V_{S30}$, measured from geo-technical investigations, is provided or not in ESM. This is to estimate the site-specific terms ($\delta S2S_s$) at as many sites as possible, and then explore various site-response proxies to characterize them (Kotha et al. 2018; Weatherill et al. 2020b)

4. Choice of distance metric is $R_{JB}$ where available, otherwise the epicentral distance $R_{epi}$— but only for events with $M_W \le 6.2$. The distance range is not truncated and extends up to $R_{JB} = 545$ km

5. Only those records with high-pass filter frequency $f_{hp} \le 0.8/T$, where T is the period, for both horizontal components are kept in the regression of spectral accelerations $SA(T)$. This is to ensure that the filter does not significantly affect the response spectral values. As a result the dataset varies with the *T* (Abrahamson and Silva 1997). Further detail on waveform processing can be found in Lanzano et al. (2019a).

Following the above criteria, the number of records available for GMM regression is 18,222, from 927 events ($3.1 \le M_w \le 7.4$) recorded at 1829 stations ($0 \le R_{JB} < 545$ km). These numbers decrease to 9698, 491, and 1341 respectively due to the high-pass filter criteria (#5) when $T=8$ s.

## 3 Regionalisation datasets

In order to attempt GMM regionalisation, we first need to define regions, classify the data into regions, and quantify the regional variabilities through regression. The mixed-effects regression estimates first the *group* random-variances, and then the random-effects for individual *levels* in the group (Bates et al. 2014). For instance, in estimating the well-known between-event term ($\delta B_e$), the preliminary step is to quantify observed event-to-event ground-motion variability as the between-event random-variance ($\tau^2$). Following this,
the between-event random-effects for each event are estimated from the total-residuals. The procedure is similar for the between-site terms ($\delta S_2 S_s$), and other random-effects in a mixed-effects regression. Essentially, within the event *group*, individual events are the *levels*. It is important to note that random-effect values for the levels (in a group) can be used in three ways: (1) in level-specific predictions accounting the epistemic uncertainty of level-specific adjustment, (2) investigated for any physical phenomena and proxy parameters or (3) ignored, while treating the group random-variance as an aleatory variability or an epistemic uncertainty in a GMM logic-tree (Douglas 2018a, b; Weatherill et al. 2020c).

We emphasize that an *exploratory* analysis on random-effects (of levels in a group) is to evaluate if they concur with known geophysical properties. If yes, the random-effect group remains in the regression; otherwise the grouping needs revision. In this study, we aim to capture variabilities among events, event sources (e.g. fault systems), attenuating regions, and recording sites. Anticipating significant random-variances, we first describe the various regionalisation datasets used to formulate the random-effect groups.

## 3.1 Attenuation regionalisation

Recent GMMs have demonstrated a strong between-region variability of anelastic attenuation, which is partially attributed to spatial variability in crustal characteristics, e.g. the 1 Hz Lg-coda Q values (Cong and Mitchell 1998), shear-wave velocity (Lu et al. 2018), etc. The first generation of regionalised GMMs, however, relied on national administrative boundaries for regionalisation of anelastic attenuation (e.g. Italy, Turkey, Japan, California, etc.)—which could yield incongruous estimates of random-variances. For example, spatial variability of attenuation imaged by 1 Hz coda Q maps of regions within Italy, Turkey (Cong and Mitchell 1998), France (Mayor et al. 2018), and Europe in general (Pilz et al. 2019) is rather high. Such spatial variability cannot be contained within political boundaries. Therefore, in this study, we adopt a more geological-geophysical regionalisation developed by Basili et al. (2019) under the purview of the TSUMAPS-NEAM project.

Figure 2 shows the regionalisation and the number of records within each region polygon, as decided by the recording site location. Assuming that the anelastic attenuation is

<figure><img src="image_5.png" /><figcaption>**Fig. 2** Number of records in each of the 46 attenuation regionalisation polygons</figcaption></figure>

<figure><img src="image_6.png" /></figure>
a phenomenon dominant at far-source distances and in near-surface crustal layers, we let the recording site location decide to which region a particular record is to be assigned. Alternatively, regionalising records based on event location would cause some inconsistencies concerning event depth. Seismic waves from shallow and deep events sample different depths of the crust and upper mantle. TSUMAPS-NEAM regions used in this study are surficial only, and cannot accommodate depth dependence of attenuation. Therefore, records are allotted to the regions based on the recording surface site location.

There are 46 regions in this map (Fig. 2), which is substantially larger than K16 with Italy, Turkey, and rest of Europe as the three attenuating regions. This regionalisation scheme splits Turkey into at least eight regions, East and West Anatolia being the largest. Italy is highly fragmented as well. Figure 2 indicates that there are regions with a few thousand records (e.g. Central and Northern Apennines), and a few with less than a hundred records (e.g. only 13 records from Rhine Graben). Thus, random-variance of the *attenuation regions* group will be a quantification of the observed region-to-region variability in anelastic attenuation across these regions (each region is a level). A finer or coarser regionalisation may yield a different estimate of between-region variance. Customarily, only if the random-variance for this group is significant it would mean the regionalisation is appropriate, and that there is a significant regional variability in anelastic attenuation. Subsequently, we will inspect whether the region-specific random-effect values resemble any physical phenomenon.

## 3.2 Event localisation

Traditionally, earthquake-to-earthquake variability is captured by the between-event random-variance ($\tau^2$), and the event-specific random-effect values ($\delta B_i$) are estimated for each event in this group. Based on earlier works on fault-maturity (Bohnhoff et al. 2016; Manighetti et al. 2007; Radiguet et al. 2009), we hypothesized that events associated to particular earthquake source (or a fault system) show systematic differences in their ground-motions, and that event-to-event ground-motion variability is similar across various fault-systems.

In this study, we introduced an additional random-effect to quantify earthquake locality-to-locality variability, similar to the location-to-location variability ($\delta L_{2L_i}$) defined in Al Atik et al. (2010). We assign the ESM shallow crustal events to seismotectonic zones defined in the European Seismic Hazard Model 2020 (ESHM20) towards development of seismic source models. The seismotectonic source zonation (referred to as ‘TECTO’) is designed to be large-scale and does not attempt to resolve smaller scale seismogenic features. As such it is intended to reflect the local tectonics influencing the seismic source generation but not to resolve highly localised features, which makes it a good candidate for our purpose—and more appropriate for source localisation than the TSUMAPS-NEAM model used for attenuation regionalisation.

The random-effects group will be referred to as between-locality ($\delta L_{2L_i}$) from hereon, with each level (locality $l$) being a tectonic locality with allotted shallow crustal events. Figure 3 shows the distribution of recordings across the tectonic localities with at least one shallow crustal event associated to them. Once again, this group is introduced to quantify the earthquake locality-to-locality variability of ground-motion in the dataset through the between-locality variance ($\tau^2_{L_{2L}}$), which if close to zero indicates no regional variability of earthquake characteristics.
<figure><img src="image_7.png" /><figcaption>Fig. 3 Distribution of records from shallow crustal events in each of the 56 event localisation polygons (tectonic localities)</figcaption></figure>

## 4 Functional form

A mixed-effects GMM is composed of fixed-effects and random-effects. Fixed-effects part of the GMM predicts generic ground-motions as a continuous function of predictor variables, which in this case are event magnitude $M_W$ and distance metric $R_{JB}$ (in km). Random-effects can serve as adjustments to the generic fixed-effects regression coefficients in order to produce ground-motion predictions specific to an event, tectonic locality, region, and site. Equation (1) illustrates the decomposition of observed ground-motion values into mixed-effects. We note that Eq. (1) is not an exact mixed-effects formula, but is used here only as an aid to explain the approach.

$$ \ln(GM_{e,l,s,r}) \sim N \left( \ln(\mu) + \left| \begin{array}{c} \delta B_{e,l}^0 \\ \delta L2L_l \\ \delta S2S_s \\ \delta c_{3,r} \end{array} \right|, \phi^2 \right) \quad (1) $$

In the LHS of Eq. (1), $\ln(GM_{e,l,s,r})$ is the natural-log of ground-motion ($GM_{e,l,s,r}$ which are rotD50 measures of SA(T) here) produced by an event (e) originating in a tectonic locality (l in Fig. 3), and recorded by a site (s) located in region (r in Fig. 2). These observed ground-motions are decomposed into mixed-effects. Since events are exclusive to tectonic localities, and sites to regions, the formulation contains nested mixed-effect grouping. Alongside, since events in a locality can be recorded by sites in multiple surrounding regions, the formulation contains crossed mixed-effect grouping as well. Therefore, this GMM is developed using a crossed and nested mixed-effects approach (Stafford 2014).

In the RHS of Eq. (1), $\ln(\mu)$ is the generic fixed-effects formula without any specificities related to event, locality, site, and region. The $\delta B_{e,l}^0$, $\delta L2L_l$, $\delta S2S_s$, and $\delta c_{3,r}$ are random-effect adjustments to the coefficients of the fixed-effect $\ln(\mu)$; and are aimed to be specific to the event (e), locality (l), site (s), and region (r), respectively. $\phi$ is the standard-deviation of the left-over residuals $\epsilon = N(0, \phi^2)$, which contain phenomenon
that could not be contained in the crossed and nested mixed-effects. Subsequent sections will explain these terms in more detail.

## 4.1 Fixed-effects

Equations (2) through (5) present the fixed-effects formula of the GMM. In Eq. (2), ln(μ) is the generic median prediction as a function of M<sub>W</sub> and R<sub>JB</sub>; where, e<sub>1</sub> is the generic offset, f<sub>R,g</sub>(M<sub>W</sub>, R<sub>JB</sub>) represents the geometric spreading, f<sub>R,a</sub>(R<sub>JB</sub>) represents the apparent anelastic attenuation, and f<sub>M</sub>(M<sub>W</sub>) represents the magnitude scaling of ground-motion prediction. The ln(μ) is the predicted natural-log of RotD50 (Boore 2010) measures of PGA (in gal) and PGV (in cm s<sup>-1</sup>), and the 5% damped elastic response spectral ordinates in acceleration (SA, in gal) at 34 periods ranging from 0.01 to 8 s. Therefore, the fixed-effects coefficients e<sub>1</sub>, b<sub>1</sub>, b<sub>2</sub>, b<sub>3</sub>, c<sub>1</sub>, c<sub>2</sub>, c<sub>3</sub> in Eqs. (2–5) change with the period, PGA and PGV, and are generic to all events, tectonic localities, attenuating regions, and sites.

$$ \ln(\mu) = e_1 + f_{R,g}(M_W, R_{JB}) + f_{R,a}(R_{JB}) + f_M(M_W) \quad (2) $$

$$ f_{R,g} = (c_1 + c_2 \cdot (M_W - M_{ref})) \cdot \ln \sqrt[2]{(R_{JB}^2 + h_D^2) / (R_{ref}^2 + h_D^2)} \quad (3) $$

$$ f_{R,a} = \frac{c_3}{100} \cdot \left( \sqrt[2]{R_{JB}^2 + h_D^2} - \sqrt[2]{R_{ref}^2 + h_D^2} \right) \quad (4) $$

$$ f_M = \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_h < M_W \end{cases} \quad (5) $$

The fixed-effects component of this GMM remains similar to that of K16, but with a few minor changes based on non-parametric analyses. In Fig. 4, the top row of the plot shows non-parametric scaling of SA(T = 0.1s) with R<sub>JB</sub>; wherein, for clarity, the data is split into magnitude bins and plotted in separate panels. The bottom row of the plot shows non-parametric scaling of SA(T = 0.1 s) with M<sub>W</sub>; where, the data is split into distance bins and plotted in separate panels. The smooth (coloured) curves are the loess fits (Jacoby 2000) illustrating non-parametric scaling of SA(T = 0.1 s) with magnitude and distance. Given the variety of hypocentral depths in the dataset, and that we are using the depth-insensitive R<sub>JB</sub>, we categorized the events into three depth bins (D < 10 km, 10 km ≤ D < 20 km, 20 km ≤ D), and plotted unique non-parametric attenuation curves for each bin. We note that, although the non-parametric plots presented here are for SA(T = 0.1 s) and the shape of loess fits varies gradually with periods, the differences were not large enough to make changes to the GMM functional form. Based on such plots we introduced some changes into the GMM with respect to K16:

1. In the SA(T = 0.1 s) versus R<sub>JB</sub> panels, it is evident that deeper events (20 km ≤ D) have a longer near-source saturation plateau—a feature controlled by the so-called h-parameter in GMMs (h<sub>D</sub> in Eqs. 3 and 4). Shallower events have a shorter near-source saturation plateau, and a steeper decay with distance. However, these differences are only prominent at R<sub>JB</sub> ≤ 30 km.
<figure><img src="image_8.png" /><figcaption>**Fig. 4** Non-parametric plots of SA(T = 0.1s) versus R<sub>JB</sub> (top panels) for various magnitude ranges, and SA(T = 0.1s) versus M<sub>W</sub> (bottom panels) for various distance ranges. The smoothed curves are loess fits to the data, color-coded to distinguish scaling for events in three hypocentral depth bins</figcaption></figure>

The marginally larger SA(T = 0.1 s) values of 20 km ≤ D event at R<sub>JB</sub> > 30 km are also reasonable (Derras et al. 2012). Deeper events are closer to Moho; producing weaker surface ground-motions in the epicentral zone due to longer travel paths of direct arrivals, and relatively stronger ground-motions at far-source distances from the more efficient propagation in deeper crustal layers.

Based on these observations, we allowed the h-parameter to vary with the depth bins, as indicated by the subscript D in $h_D^2$ of Eqs. (3) and (4). In K16 GMM, this parameter was depth-independent, and had no subscript. Moreover, instead of making this parameter a regression coefficient, we keep the regression linear by assigning a priori values (based on preliminary non-parametric analyses and nonlinear regression trials), which are independent of magnitude and period: $h_D$ = 12 km for deep 20 km ≤ D events, $h_D$ = 8 km for events of intermediate depth 10 km ≤ D < 20 km, and for shallow events with D < 10 km we assign $h_D$ = 4 km

2. In addition to the above, since ground-motions appear reasonably depth independent at R<sub>JB</sub> ≥ 30 km, we set R<sub>ref</sub> = 30 km instead of the 1 km in K16. M<sub>ref</sub> in Eq. (3) is the reference magnitude, and remains the same as in K16 i.e., M<sub>ref</sub> = 4.5

3. In the SA(T = 0.1s) versus M<sub>W</sub> panels, although the depth dependence of near-source attenuation is not evident, we do observe saturation towards large magnitudes. In the panel showing non-parametric ground-motion scaling with M<sub>W</sub>, the evidence suggests that saturation begins at M<sub>W</sub> ≥ 6.2. However, this is a feature most noticeable at R<sub>JB</sub> ≤ 30 km and at short periods; here SA(T = 0.1 s). Towards longer periods and at longer distances, the saturation is less pronounced. Therefore, unlike in K16 where the hinge-magnitude was set as M<sub>h</sub> = 6.75, in this GMM we set M<sub>h</sub> = 6.2. M<sub>h</sub> is period independent; saturation or otherwise beyond M<sub>h</sub> is captured by b<sub>3</sub> in Eq. (5).
## 4.2 Random-effects

In the mixed-effects formulation of Eq. (1), $ln(μ)$ is the fixed-effects component. Equations (2) through (5) describe its constituents, wherein $e_1, b_1, b_2, b_3, c_1, c_2, c_3$ are the fixed-effects regression coefficients. These generic fixed-effect coefficients can be adjusted using random-effect estimates to predict event, tectonic locality, attenuation region, and site specific ground-motions. Regarding the random-effects components:

1. $τ_{c3}$ quantifies the between-region variability of anelastic attenuations across the attenuation region group described earlier (in Fig. 2). This means that, along with a generic $c_3$ and random-variance $τ_{c3}^2$, region-specific adjustments $δc_{3,r}$ are estimated as random-effects. These random-effects follow a Gaussian distribution $Δc_3 = N(0, τ_{c3}^2)$. The generic fixed-effect coefficient $c_3$ in Eq. (4) can be adjusted with a region's $δc_{3,r}$ to achieve region-specific coefficient for that region $r$, as in $c_{3,r} = c_3 + δc_{3,r}$. Unlike the generic $c_3$, $c_{3,r}$ is the apparent anelastic attenuation term that varies with region (subscript $r$ identifies regions in Fig. 2). In K16, $r$ identified the regions Italy, Turkey, and the rest of pan-Europe as Other. This grouping ensured that each of the regions had sufficient data to estimate statistically reliable $δc_{3,r}$. Accordingly, $τ_{c3}$ quantified the regional variability of anelastic attenuation when RESORCE dataset is grouped into Italy, Turkey, and Other. ESM contains much more data from Italy and Turkey, and also several other nations. Since the number of regions is 15 times (46 regions) that of K16 (three regions), the quantified regional variability in anelastic attenuation, in terms of $τ_{c3}$, is also larger than that of K16. Regions with sufficient data also benefit from a well-constrained region-specific adjustment $δc_{3,r}$.

2. Between-locality variability of observed ground-motions are captured by the random-effect $ΔL2L = N(0, τ_{L2L}^2)$, where the mixed-effects regression quantifies the variability as $τ_{L2L}$, at each $T$, and earthquake locality-specific terms as $δL2L_l$ (subscript $l$ identifies regions in Fig. 3). This random-effect can be used to adjust the fixed-effect coefficient $e_1$ in Eq. (2), as in $e_{1,l} = e_1 + δL2L_l$, to achieve ground-motion predictions specific to the tectonic locality of the events' origin

3. Event-to-event variability in this GMM is the *traditional* between-event random-effects $ΔB_e = N(0, τ^2)$ filtered for between-locality variability $ΔL2L = N(0, τ_{L2L}^2)$, and is now captured by the $ΔB_e^0 = N(0, τ_0^2)$; where, for an event $e$ located in tectonic locality $l$, the event-specific term can be seen as $δB_{e,l}^0 ≈ δB_e - δL2L_l$. $τ_0$ is the generic event-to-event variability corrected for locality-to-locality variability, and does not vary with tectonic locality $l$

4. Site-to-site response variability is captured by the site-specific random-effects $ΔS2S = N(0, φ_{S2S}^2)$. The potential of $δS2S_s$ in site-specific GMMs is well-known, and are useful in studying regional differences in site-response scaling with $V_{S30}$ (time-averaged shear-wave velocity in 30 m top-soil) as in K16 or other site-response proxies (Kotha et al. 2018; Weatherill et al. 2020b)

5. The left-over residuals $ε = N(0, φ^2)$ contain the unexplained natural variability of ground-motion observations, and thus represent the *apparent* aleatory variability of the model. These residuals can be investigated for less dominant phenomenon, such as the anisotropic shear-wave radiation pattern (Kotha et al. 2019)

In all, this GMM has four random-effect groups, i.e. one degree-of-freedom more than K16, to explain more than 15 times the data. Those common with K16 are refined with a
more physical regionalisation scheme, and greater geographical coverage of shallow crus-
tal events and recording sites. With this configuration of mixed-effects GMM, we run a
robust linear mixed-effects regression (Koller 2016) independently for the 36 RotD50 IMs
of 5% damped SA for T = 0.01 – 8 s, PGA and PGV. Along with the regression coeffi-
cients, we estimate and provide also the fixed-effects variance–covariance matrices needed
to estimate the GMM epistemic uncertainty (Atik and Youngs 2014; Bindi et al. 2017) and
to update the GMM coefficients in a Bayesian framework (Kowsari et al. 2019, 2020; Kühn
and Scherbaum 2015).

**5 Regression method adapted to large datasets**

Unlike K16, where the GMM is regressed using an ordinary least-squares mixed-effects
regression algorithm (Bates et al. 2014), in this study we employ a robust mixed-effects
regression algorithm (Koller 2016). As any other real-life data, ground-motion data may
contain outliers and other contaminations. Even minor contamination may drive the clas-
sical ordinary least-square mixed-effect estimates away from those without contamination.
Robust linear mixed-effects (rlmm) regression is then quite useful in limiting the influence
of outlier events, sites, and records on the GMM median and variances. It is important to
note that, these outlier data points are not entirely removed from the GMM regression but
are simply down-weighted.

A feature of rlmm relevant to this study is that the random-effects (attenuation regions, events, event localities, and sites) and residuals (records) with values beyond ±1.345 standard-deviations of their respective normal random-distributions are assigned progressively lower weights (<1); whereas, in ordinary least-squares all data is assigned unit weight. Any event, site, and record with non-unit weight are considered a possible outlier, and needs to be examined for its peculiarity. Although it is possible to tune the regression parameters, e.g. those of the Huber loss function (Huber 1992) in robust regressions, we chose to remain with the suggested default values (Koller 2016) optimized for efficiency and robustness in detecting outliers.

**6 Results and discussion**

The regression results comprise of fixed-effect coefficients and covariance matrices, ran-
dom-effect values including weights and standard-errors, residuals, and variances. It is cus-
tomary to check the behaviour of GMM fixed-effects component, and its epistemic uncer-
tainty in various magnitude and distance ranges. Random-effects and residuals are checked
for any noticeable biases or trends with predictor variables. We discuss them separately.

**6.1 Fixed-effects**

Figure 5 shows the GMM’s median SA(T = 0.1 s) and SA(T = 1 s) (in gal) predictions over magnitude and distance ranges. Along with the median prediction (lines), its epistemic uncertainty in terms of asymptotic standard-deviation (±σµ) is shown as well (coloured ribbons).
<figure><img src="image_9.png" /><figcaption>**Fig. 5** GMM median predictions of SA(T = 0.1 s) and SA(T = 1 s). The left panel shows predicted scaling of SA(T = 0.1 s). (top row) and SA(T = 1 s) (bottom row) with distance, for shallow (*D* < 10 km), intermediate (10 km ≤ *D* < 20 km), and deep (20 km ≤ *D*) events in regions with *average*, *slower*, and *faster* than median anelastic decay at long distances. The right panel shows scaling of ground-motion with magnitude for the same combinations. In both panels, the epistemic uncertainty on GMM median is shown by the orange ribbon around prediction for *intermediate* depth events in *average* attenuating regions. Also shown are ground-motion observations corrected for δL2Lᵢ, δB₀ᵢ, and δS2Sₛ random-effects as grey markers</figcaption></figure>

### 6.1.1 Distance scaling

First, we discuss the predicted SA(T = 0.1 s) and SA(T = 1 s) scaling with $R_{JB}$ in the left panel of Fig. 5. Looking at the curves for M4, we notice the impact of depth-dependent $h_D$ in rendering a shorter near-source saturation (plateau ~0 to 3 km) for shallower events, compared to the intermediate depth events (plateau ~0 to 5 km) and deeper events (plateau ~0 to 10 km). The three curves merge at about 30 km, which is our $R_{ref}$. Beyond $R_{ref} = 30$ km, the depth-dependence of distance scaling is non-existent. The predictions show good agreement with the non-parametric trends in Fig. 4.

Beyond the $R_{ref}$, we see the region-dependent anelastic adjustments coming into play, which show the impact of adjusting $c_3$ with $\delta c_{3,r} = 0, -\tau_{c3}, +\tau_{c3}$. A region with faster than average attenuation will have a $c_{3,r}$ more negative (with $\delta c_{3,r} = -\tau_{c3}$) than generic average, and vice versa for slower attenuation (with $\delta c_{3,r} = +\tau_{c3}$). Which of the 46 regions in pan-European apparently attenuate faster or slower than pan-European average (with $\delta c_{3,r} = 0$) will be discussed in the following sections, and elaborated in a follow-up study. Effect of
$\delta c_{3,r}$ at $R_{ref} < 30$ km is negligible, as it should be. The 27 curves appear as three on either side of $R_{ref} = 30$ km, because near-source ($h_D$) and far-source ($\delta c_{3,r}$) adjustments have their exclusive domains of influence.

In Fig. 5, we also show the epistemic uncertainty on median predictions. The orange ribbon is almost too thin to be noticeable for M4 and M5.5 predictions. Only for M7 and larger events, the ribbon is visibly wide because of the limited data from large magnitude events at near-source distances in the ESM dataset, and the uncertainty on the coefficient $b_3$, which describes the magnitude scaling beyond $M_h=6.2$.

### 6.1.2 Magnitude scaling

In the right panels of Fig. 5 we show the $SA(T = 0.1 \text{ s})$ and $SA(T = 1 \text{ s})$ scaling with $M_W$. Here as well we show 27 curves, but this time for $R_{JB} = 10, 50, 150$ km instead of the three $M_W$ values. The features in scaling with distance discussed in reference to the left panels also prevail here; $\delta c_{3,r}$ is effective at $R_{JB} = 150$ km, while $h_D$ is effective at $R_{JB} = 5$ km, and neither are at $R_{JB} = 50$ km. More important in this context, is the difference in scaling with magnitude at $R_{JB} = 10$ km compared to $R_{JB} = 50, 150$ km. Evidently, the scaling of $SA(T = 0.1 \text{ s})$ is more gradual (less steep) at near-source distance for $M_W < M_h$ and oversaturates for $M_W \ge M_h$. This is a known physical behaviour wherein, near-source groundmotions, especially the short period SAs, are less sensitive to $M_W$ (Campbell 1981; Schmedes and Archuleta 2008). A few previous GMMs observed the same with various datasets. Boore et al. (2014) and K16 allowed oversaturation at large magnitudes ($M_W \ge 6.75$), but whether this is realistic or not needs to be verified with specially compiled near-source ground-motion datasets, e.g. the one recently made available by Pacor et al. (2018). The oversaturation effect, however, disappears at longer periods, as shown here for $SA(T = 1 \text{ s})$.

Figure 6 shows the predicted response spectra for the same scenarios shown in Fig. 5. The effect of event depth and regionalized anelastic attenuation are the same as in Fig. 5. An interesting feature is that the peak of the response spectra shifts towards longer periods for larger magnitude events, reflecting the decreasing corner-frequency with increasing magnitude of a Brune (1970) source spectra.

<figure><img src="image_10.png" /><figcaption>**Fig. 6** Predicted response spectra of the GMM for various scenarios, differentiated by colour for shallow (D < 10 km), intermediate (10 km ≤ D < 20 km), and deep (20 km ≤ D) events in the left panel. In the right panel are the response spectra differentiated by regional anelastic attenuation (Slower, Average, Faster)</figcaption></figure>
An oddity in the response spectra is the *kink* appearing between $T = 3.5$ s and $T = 4.5$ s, which is most likely from the sudden loss of approximately 3000 usable records (in regression) due to the high-pass filtering frequency criterion #5 in Ground-Motion Data and Selection Criteria. Similar kinks can be seen in a few other subsequent plots as well. Traditionally, GMM fixed-effect regression coefficients are *smoothened* in a post-processing step, in order to smoothen the predicted response spectra and to constrain the GMM behaviour in sparsely sampled magnitude-distance ranges, e.g. Abrahamson et al. (2014). However, in this study we chose not to alter the fixed-effect coefficients, and to maintain full consistency with their variance–covariance matrices; which will be necessary in testing and updating of the GMM.

## 6.2 Random-effects and residuals

Figure 7 shows the random-effect and residual standard-deviations of the GMM. The total ergodic standard-deviation of the GMM is $\sigma = \sqrt{\tau_{L2L}^2 + \tau_0^2 + \phi_{S2S}^2 + \phi^2}$, when considering all source and site variabilities as aleatory. The solid lines in this plot correspond to the variance estimates of the GMM from ESM dataset (this study), while the dashed lines indicate those of the K16 GMM from RESORCE dataset. Note that the K16 GMM does not regionalise the events, and thus between-locality variability component $\tau_{L2L}$ is absent. For the sake of comparison between the two models, we combine the two source related random-variances of this study into the traditional between-event variability $\tau = \sqrt{\tau_{L2L}^2 + \tau_0^2}$.

We show that the total standard-deviation $\sigma$ of the new GMM is considerably larger than that of K16 GMM. The largest increase in $\sigma$ is at short-moderate periods, and is from the increase in between-event $\tau$ and between-site $\phi_{S2S}$ variabilities. The increase in between-event variability can be attributed to the increased regional diversity of earthquakes, the increase of number of $4 \le M_W \le 5.5$ events from 164 in RESORCE to 778 in ESM, and the additional 76 events with $3 \le M_W < 4$ in ESM. Between-site variability is clearly from the increase in number of recording sites from 385 in RESORCE to 1829 in ESM. The residual variability $\phi$ of the new GMM, however, is equal or smaller than that of K16, despite the increase in regional diversity, the sample size, and the recording distance range

<figure><img src="image_11.png" /><figcaption>**Fig. 7** Residual and random-effect variance estimates of the GMM for $T = 0.01 – 8$ s (solid lines). The colours identify the various random-effect and residual components. For comparison, variance estimates for the K16 GMM from RESORCE dataset (Kotha et al. 2016) are overlain in dashed lines</figcaption></figure>
from 300 to 545 km. In our view, this plot emphasizes the need to move from ergodic to partially non-ergodic region- and site-specific ground-motion predictions to inhibit the effect of increased σ on PSHA (Anderson and Brune 1999; Bommer and Abrahamson 2006; Kotha et al. 2017; Rodriguez-Marek et al. 2013).

### 6.2.1 Anelastic attenuation variability

Figure 8 shows the region-specific $\delta c_{3,r}$ adjustments (coloured lines and ribbons) of the 45 regions for $T = 0.01 – 8$ s. Most of these curves lay within the $\pm \tau_{c3}$ (red lines) bounds. As indicated by $\tau_{c3}$, the region-to-region variability of apparent anelastic attenuation is the highest at short periods, and decreases gradually towards longer periods. High-frequency ground-motions are more susceptible to strong anelastic decay in the crust, which could be related to the crustal properties. Long period SAs are not affected as much by anelastic attenuation, therefore regional differences are relatively smaller, and $\tau_{c3}$ is smaller by ~50% at $T \ge 1$ s.

In RESORCE dataset, K16 observed that short period SAs attenuate faster in Italy than in Turkey, which was observed earlier in NGA-West2 dataset by Boore et al. (2014), and confirmed later in Fourier domain by Bora et al. (2017). However, these observations were based on distinguishing regions by administrative boundaries and not geological or geophysical features. Since $\tau_{c3}$ estimated from the new regionalisation is not zero, we can assert that regional variability of anelastic attenuation exists. Of course, regions with fewer data have larger epistemic uncertainty (standard-error) on their $\delta c_{3,r}$, but the largest epistemic uncertainty is always lower than the aleatory $\tau_{c3}$, and decreases with increasing data.

Figure 9 shows the spatial variability of apparent anelastic attenuation as captured by the regionalisation model of Fig. 2. In this figure, red polygons identify regions with slower than pan-European average of anelastic attenuation $c_3$ i.e. $\delta c_{3,r} > 0$, and $\delta c_{3,r} < 0$ for the faster attenuating blue regions. Regions with insufficient data, and thereby large epistemic uncertainty on their $\delta c_{3,r}$, are white in color, and given $\delta c_{3,r} = 0$. There are few interesting features in these maps:

1. Regions with similar attenuation characteristics are spatially clustered, although this is period dependent. In general, regions characterized with high seismic activity (e.g. Italy

<figure><img src="image_12.png" /><figcaption>**Fig. 8** $\delta c_{3,r}$ for $T = 0.01 – 8$ s. Each line corresponds to one of the 46 attenuation regions, with colours indicating their weight in robust regression. Overlain red curves mark the $\pm \tau_{c3}$ values. Regions with $\delta c_{3,r}(T)$ beyond $\pm 1.345\tau_{c3}(T)$ are given lower than a unit weight</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 9** Variation of $\delta c_{3,r}$ across the 46 attenuating regions ($r$) for $T = 0.1$ s (top panel) and $T = 1$ s (bottom panel). Blue polygons identify regions with anelastic attenuation faster than the pan-European average, red polygons locate regions with slower attenuation, and white polygons are regions close to the average. Regions with fewer ground-motion observations, thereby larger epistemic error on $\delta c_{3,r}$, are more transparent and appear white. Note that the color scale is limited to vary between $\pm 3 * \tau_{c3}(T)$ for each period $T$.</figcaption></figure>

and Greece) show strong attenuation compared to those with lower seismic activity (e.g. central Europe).

A Bayesian clustering of regions with similar attenuation characteristics is proposed in the companion paper by Weatherill et al. (2020c); wherein, similar regions ($r$) are clustered ($cl$) and assigned cluster-specific $\delta c_{3,cl}$ with lower standard-error and smaller cluster-specific $\tau_{c3,cl}$. The epistemic uncertainty (standard-error) on cluster-specific $\delta c_{3,cl}$ is smaller than that of $\delta c_{3,r}$ due to the accumulation of data from multiple regions; while the cluster-specific $\tau_{c3,cl}$ is smaller than the overall $\tau_{c3}$ due to a smaller within-cluster variability of region-specific $\delta c_{3,r}$.

2. The best-sampled regions are in central Italy with 5703 records from Northern and central Apennines W (West), 3438 records from Northern and central Apennines E (East). While attenuating faster than the pan-European average, there appears to be a strong contrast between these adjacent regions.
3. The fastest attenuation in the Aegean Sea is observed in the Gulf of Corinth, where the sites record highly attenuated ground-motions traversing across the Aegean volcanic arc

A contrast in short period attenuation is apparent around the Alps regions. In addition, some differences are noticeable between west, north, and central Anatolia. Although not conclusive, rapidly changing crustal thickness (Grad et al. 2009) and associated crustal properties may partially explain the rapid change in attenuation properties in these regions. Regional variability of anelastic attenuation may in fact be a combination of variability of crustal shear-wave velocity, crustal quality factor (e.g. coda Q), mantle temperature influencing the rigidity of the crust, and other parameters. Further elaboration is preferred in Fourier domain, rather than with response spectra here, in a follow-up study.

6.2.2 Source variability

Earthquake source variability is split in two random-effect groups: between-event $\Delta B_e^0 = N(0, \tau_0^2)$ and between-locality $\Delta L2L = N(0, \tau_{L2L}^2)$. Between-event terms can be estimated for the recorded earthquakes, but are difficult to predict for prospective earthquakes because of their spatiotemporal variability. Even when correlated with stress-drop or a source parameter that can explain the relative differences in ground-motions, prediction of stress-drop for the next event is not yet possible. Traditionally therefore, between-event variability is considered purely aleatory. Between-locality random-effect is intended to quantify a portion of the between-event spatial variability into tectonic localities.

Figure 10 shows the $\delta L2L_t(T = 0.01 - 8 \text{ s})$ for 56 tectonic localities. As always, the standard-error (or epistemic uncertainty) of $\delta L2L_t$ values is smaller than the $\tau_{L2L} \cdot \tau_{L2L}$ values are non-negligible, and the GMM fit improves (based on analysis of variance) with introduction of between-locality random-effect group. Therefore, we consider it an effective random-effect grouping. Although it is tempting to relate $\delta L2L_t$ to some source parameter, it is preferably done in the Fourier domain in a follow-up study.

Figure 11 maps the various tectonic localities (indexed *l*) color coded to their $\delta L2L_t(T = 0.1 \text{ s})$ and $\delta L2L_t(T = 1 \text{ s})$ values. In the top panel, corresponding to $\delta L2L_t(T = 0.1 \text{ s})$, a clear difference between central Italy and north-western Anatolia can be seen. Apparently, earthquakes located in central Apennine region produce substantially lower short period ground-motions than those in north-western Anatolia.

<figure><img src="image_14.png" /><figcaption>**Fig. 10** $\delta L2L_t$ for $T = 0.01 - 8 \text{ s}$. Each line corresponds to one of the 56 tectonic localities, with colours indicating their weight in robust regression. Overlain red curves mark the $\pm \tau_{L2L}$ values. Tectonic localities with $\delta L2L_t(T) \text{ beyond } \pm 1.345\tau_{L2L}(T)$ are given a lower weight than one</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. 11** Variation of $\delta L_{2L_i}$ across the 56 tectonic locality ($l$) for $T = 0.1$ s (top panel) and $T = 1$ s (bottom panel). Blue polygons identify tectonic localities that produced earthquakes less energetic than the pan-European average, red polygons identify localities with more energetic earthquakes, and white polygons are localities close to the average. Tectonic localities with fewer ground-motion observations, thereby larger epistemic error on $\delta L_{2L_i}$, are more transparent and appear white. Note that the color scale is limited to vary between $\pm 3 * \tau_{L_{2L}}(T)$ for each period $T$</figcaption></figure>

Similarly, there is an apparent distinction between the central Apennines and Po-plain earthquakes. The central Apennines tectonic locality contains the recent M6.5 Norcia earthquake (2016) and associated shocks, and several historic earthquakes in that region. The Po-plain locality contains data from the substantially stronger M6.45 Friuli earthquake (1976) and a few recent earthquakes. At a first glance, it may appear as if the spatial patterns are due to the predominant focal mechanisms in the regions, but the diversity of focal mechanisms within each region (especially among smaller events) disputes this hypothesis. Spectral decomposition of the ESM dataset by Bindi and Kotha (2020) revealed a much larger stress-drop of the M6.45 Friuli earthquake compared to M6.5 Norcia earthquake—the $\delta L_{2L_i}$ of regions containing these earthquakes also contrast similarly (Bindi et al. 2019). Once again, correlating the random-effects obtained in response spectra domain (as here) to the earthquake parameters (such as stress-drop)
may not hold well at short periods. We reserve further elaboration of random-effects
with the Fourier version of the ESM dataset.

Between-event variability, now partially corrected for spatial variability through
$\Delta L2L = \mathcal{N}(0, \tau_{L2L}^2)$, is quantified in the distribution $\Delta B_e^0 = \mathcal{N}(0, \tau_0^2)$. Customary checks for
$\delta B_{e,l}^0$ include dependencies on magnitude, depth, and style-of-faulting. A few points to note
here:

1. EMEC estimates of $M_W$ are used in the GMM regression. However, uncertainties in $M_W$ are ignored despite their known impact on the $\tau_0$ estimates (Kuehn and Abrahamson 2017)

2. Depth in our case is hypocentral depth of the event. Buried ruptures are likely to produce stronger ground-motions than the exposed ruptures. This phenomenon is introduced in some of the NGA-West2 GMMs e.g., Abrahamson et al. (2014) modelled the events with depth-to-top-of-rupture $z_{tor} \ge 20$ km to have 0.75–3 times larger SAs than exposed ruptures, depending on the period

3. In this study, we use the typical style-of-faulting classification provided in the ESM dataset: Normal (NF), Thrust (TF), Strike-slip (SS), and Unknown (U). However, instead of introducing a style-of-faulting random-effect in the GMM regression, we queried the $\Delta B_e^0$ and $\Delta L2L$ distributions and found no systematic differences. This result is consistent with the fact that style-of-faulting factors cannot be quantified using isotropic factors with no azimuthal variations (Kotha et al. 2019). Centroid Moment Tensor (CMT) solutions are however available for too few events to treat style-of-faulting with azimuth dependent factors as in Kotha et al. (2019). Even if available, the diversity of crustal structure across pan-Europe makes it difficult, without substantially more metadata (e.g. take-off angles) in the ESM dataset

The left column of Fig. 12 shows $\delta B_{e,l}^0 (T = 0s, 0.1s, 1s)$ trends with $M_W$. No significant offsets or trends imply the magnitude scaling of the GMM sufficiently captures the magnitude dependence of SAs. The error bars show the median absolute deviance (MAD) of $\delta B_{e,l}^0$ within bins of size $M_W = 0.5$. MAD is a robust estimate of variance when normality of distribution is not necessarily satisfied within each bin. The MAD estimates appear to be magnitude dependent, indicating heteroscedasticity of $\tau_0$. However, we do not provide a heteroscedastic variance model now without investigating first its significance; given the small number of large magnitude events compared to small-moderate magnitude event sample (see Fig. 1 for data distribution).

The right column of Fig. 12 shows the $\delta B_{e,l}^0 (T = 0.01s, 0.1s, 1s)$ trends with depth. We see no significant trend with depth. Although deeper events are likely to produce stronger ground-motions, there exists several shallower events with similarly large $\delta B_{e,l}^0$, and correcting $\delta B_{e,l}^0$ for depth shows no remarkable reduction in $\tau_0$—as indicated by the bin-wise $\tau_0$ values (error bars in right column of Fig. 12).

6.2.3 Site-response variability

As in K16 GMM, the KiK-net (Dawood et al. 2016) based Kotha et al. (2018) GMM, and Brooks et al. (2020) study on North-sea ground-motion data, we did not introduce a site-response scaling parameter in the fixed-effects for two reasons: (1) Only 419 of the 1829 sites have measured (and not inferred) $V_{S30}$ available in the ESM, and using data from only the sites with measured $V_{S30}$ leads to a large reduction in data for GMM regression; (2)
<figure><img src="image_16.png" /><figcaption>**Fig. 12** $\delta B_{e,1}^0$ versus $M_W$ (left column) and *Depth* (right column) for $T = 0 \text{ s}$ $0.1 \text{ s}$, $1 \text{ s}$. The mean (red line) and median absolute deviance (red error bars) of $\delta B_{e,1}^0$ within the bins of size $0.5 M_W$ and $5 \text{ km}$. The two horizontal dashed lines mark the $\pm \tau_0$ for the three periods</figcaption></figure>

investigation of site-to-site variability can be performed in a subsequent step. For instance, depending on the application, one can regress a relation between $\delta S2S_s \sim V_{S30}$ measured or inferred from topographic slope and geology (Thompson et al. 2014; Vilanova et al. 2018; Wald and Allen 2007) or relate $\delta S2S_s$ directly to topographic slope and geology at a regional scale (Crowley et al. 2019; Weatherill et al. 2020a, b).

Given the importance of site-response in hazard and risk assessments, the complexities in finding a compromise between a site-response proxy, its availability at different regional scales, and the propagation of uncertainties from GMMs to risk assessments, we intend to perform a separate investigation. Nevertheless, a database of $\delta S2S_s(T = 0.01 - 8s)$ for the 1829 ESM sites along with their robust regression weights and standard-errors is provided.

Note that we always use the measured $V_{S30}$ in the analysis shown in this study, and never the inferred values. Early analyses of $\delta S2S_s(T = 0.0, 0.1, 1.0 \text{ s})$ trends with measured $V_{S30}$ and topographic slope are presented in Fig. 13. The left column of Fig. 13 plots the mean and MAD (robust standard-deviation estimate) of $\delta S2S_s(T = 0.0, 0.1, 1.0 \text{ s})$ within ranges of $V_{S30}$ coinciding with Eurocode8 site classes A, B, C, and D, and *slope* divided into 9 bins of equal width between 0.001 and 1.000 m m$^{-1}$. Trend with $V_{S30}$ (not measured but)
<figure><img src="image_17.png" /><figcaption>**Fig. 13** Plot showing the trend of $\delta S2S_s(T = 0, 0.1, 1s)$ (from top to bottom rows) with $V_{S30}$ (left column) and topographic slope at site location (right column). Note that, only 419 sites have measured $V_{S30}$, compared to 1829 sites with slope. At short-moderate periods ($T=0, 0.1s$), the binned mean and MAD (red error bars) show a poor non-parametric trend, which improves at longer periods ($T=1s$). The blue line represents the parametric forms (Eqs. 6 and 7) fit to the data</figcaption></figure>

inferred from a $V_{S30} \sim$ slope correlation model of Wald and Allen (2007) is comparable to the trend with slope from which it is inferred in the right column of Fig. 13.

The correlation $\delta S2S_s \sim V_{S30}$ (measured only) at short-periods $T = 0.0, 0.1s$ is rather poor, as indicated by similar mean and MAD (the robust standard-deviation) for $V_{S30} < 800 m/s$ in the left column of Fig. 13. It appears that short-period site-responses of soft rock, stiff, and soft soils (EC8 class B, C, D) in this dataset are not adequately distinguishable. However, it is interesting to see that the mean of EC8 class A ‘rock’ sites with $V_{S30} \ge 800 m/s$ is much lower than the rest, along with a considerably larger variability. Short-period linear soil-response of rock sites is known to be highly variable compared to softer soils, whose nonlinear soil-response may suppress the otherwise high variability from linear-only amplification (Bazzurro and Cornell 2004). On the contrary, the long-period site-response of rock sites is less variable than that of softer soils (e.g. $180 m/s \le V_{S30} < 800 m/s$). Similar inferences can be drawn from the $\delta S2S_s \sim$ slope plots in the right column of Fig. 13, where higher slopes are (usually) indicative of rock sites on steep hillsides, and lower slopes at softer sites located on flatter sediments.
For completeness, along with a database of δS2S_s (T = 0.0–8.0 s) for the 1829 ESM sites, we derive a continuous empirical models for both δS2S_s ~ V_{S30}(measured) and δS2S_s ~ slope correlations seen in Fig. 13. We chose a quadratic functional form instead of the traditional piecewise linear function; as shown in Eqs. (6) and (7). Here, the measured V_{S30} is in m s⁻¹ and slope in m m⁻¹, the regression coefficients g₀, g₁, g₂ are different for SR^{V_{S30}} and SR^{slope}, and change with period. The residuals from δS2S_s ~ V_{S30} (measured) and δS2S_s ~ slope regressions are δS2S_s^{V_{S30}} and δS2S_s^{slope}, with robust standard-deviations φ_{S2S}^{V_{S30}} and φ_{S2S}^{slope}, respectively.

$$SR^{V_{S30}} = g_0 + g_1 \cdot \ln\left(\frac{V_{S30}}{800}\right) + g_2 \cdot \left(\ln\left(\frac{V_{S30}}{800}\right)\right)^2 \quad (6)$$

$$SR^{slope} = g_0 + g_1 \cdot \ln\left(\frac{slope}{0.1}\right) + g_2 \cdot \left(\ln\left(\frac{slope}{0.1}\right)\right)^2 \quad (7)$$

Robust linear fits using an M estimator (Venables and Ripley 2002), at each of 34 periods between T = 0.01s – 8 s, PGA(T = 0 s) and PGV, are derived for δS2S_s ~ V_{S30} correlation of 419 sites with measured V_{S30} available, and δS2S_s ~ slope of the 1829 sites with slope derived from digital elevation models provided by Shuttle Radar Topography Mission.

Although heteroscedastic models of φ_{S2S}^{V_{S30}} and φ_{S2S}^{slope} appear reasonable in Fig. 13, we chose not to propose one without testing its significance; given the uneven distribution of sites in different bins. Figure 13 shows the reduction in between-site variance from using V_{S30} and slope as site-response proxies (Eqs. 6 and 7) in the GMM. For comparison, the variances of K16 GMM are also shown in this plot. Note that the K16 GMM comes in two variants: one without a site-response component, and another with measured V_{S30} as a proxy for linear site-response. Note that, the K16 standard-deviations shown in Fig. 7 are those when not using V_{S30} as a parameter, while those in Fig. 14 are when using V_{S30} as site-response proxy—hence, lower in the latter.

<figure><img src="image_18.png" /><figcaption>**Fig. 14** Between-site and total variance estimates (φ_{S2S}, σ) of the GMM for T=0.01–8 s (solid lines) compared with those of K16 GMM from RESORCE dataset (dashed lines). Reduction of (φ_{S2S}, σ) to (φ_{S30}^{V_{S30}}, σ^{V_{S30}}) using V_{S30}, and to (φ_{S2S}^{slope}, σ^{slope}) using Slope as site-response proxies is shown here. Note that the (φ_{S2S}, σ) of K16 GMM are those using V_{S30} as site-response proxy, and are smaller than those shown in Fig. 7, which are from the K16 GMM version without a site-response parameter</figcaption></figure>
A significant reduction in between-site standard-deviation can be achieved using an efficient site-response proxy or a combination of proxies. Since only a few sites are provided with measured $V_{S30}$ values, $\phi_{S2S}^{V_{S30}}$ is substantially smaller than $\phi_{S2S}^{slope}$. For a new site with $V_{S30}$ or slope available, Eq. (6) or (7) can be appended to the $\ln(\mu)$ in Eq. (2), while replacing $\phi_{S2S}$ (when estimating $\sigma$) with $\phi_{S2S}^{V_{S30}}$ or $\phi_{S2S}^{slope}$, respectively. For site with neither site-response proxy available, but with sufficient strong ground-motion recordings, the site-specific $\delta S2S_s$ term can be estimated and used for site-specific ground-motion predictions.

### 6.2.4 Aleatory variability

The last component of the GMM is the apparent aleatory variability, quantifying the natural randomness of the ground-motion data—that which is not captured by the mixed-effects. The aleatory residuals, $\epsilon = N(0, \phi^2)$, are tested for event depth and recording site distance dependencies in Fig. 15. For all periods ($T = 0.0, 0.1, 1.0$ s), we observe no significant trends in binned means and MAD, which implies the $f_{R,g}(M_W, R_{JB})$ and $f_{R,a}(R_{JB})$ of the GMM (Eqs. 2 and 3) explain the distance scaling of ground-motion sufficiently well. Heteroscedasticity of $\phi$ is not evident either in these plots.

The moderate-long period residuals of this dataset show a clear evidence of anisotropic shear-wave radiation pattern in the near and intermediate distance ranges of $R_{JB} \le 80$ km, similar to that reported by Kotha et al. (2019) with the KiK-net (Dawood et al. 2016; Kotha et al. 2018) and NGA-W2 (Ancheta et al. 2014; Boore et al. 2014) datasets and their GMMs. In addition, there is a strong evidence of SmS phases from Moho reflection in the $60 \text{ km} < R_{JB} \le 200 \text{ km}$ distance ranges, e.g. Bindi et al. (2006), especially from events deeper than 10 km and close to Moho boundary. However, we chose not to discuss these features in the response spectral domain in which the GMM residuals are estimated, but instead with those in the Fourier domain as a follow-up.

<figure><img src="image_19.png" /><figcaption>**Fig. 15** Aleatory residual $\epsilon$ trends with distance for different event depths ($D < 10$ km, $10 \text{ km} \le D < 20$ km, $20 \text{ km} \le D$) for $T = 0.0, 0.1, 1.0$ s</figcaption></figure>
# 7 Application

The GMM presented in this study has no new explanatory parameters in its functional form compared to previous pan-European GMMs. The median predictions rely only on the two generic parameters $M_W$ and $R_{JB}$, which constitute the fixed-effects. Of course, the standard-deviation estimates of the new model are significantly larger than those of K16, but this is to be expected given the 15-fold increase in data: from a greater variety of sites, tectonic localities, $M_W \le 5.5$ events, etc. To explain the variability, without introducing new parameters, we instead resolved the apparent aleatory variability into various possible contributions. Therefore, the model can be used ignoring the region-to-region, locality-to-locality, and site-to-site variabilities, but at the cost of increased aleatory variability. We provide three application possibilities:

## 7.1 Ergodic application

The first approach is by ignoring all repeatable effects, i.e. the region, locality, and site-specific adjustments. The between-locality, between-event, between-site, residual standard-deviations can be combined into an ergodic, total standard-deviation $\sigma = \sqrt{\tau_{L2L}^2 + \tau_0^2 + \phi_{S2S}^2 + \phi^2}$, as shown in Fig. 7. In this case, the regional differences in anelastic attenuation, quantified by $\tau_{c3}$, will be treated as an epistemic uncertainty on far-source distance scaling. The epistemic uncertainty on the regionalised anelastic attenuation coefficient $c_3$ in Eq. (4) is $\tau_{c3}$. This uncertainty can be handled with a GMM logic tree consisting of a slower ($c_{3,r} = c_3 + \eta \cdot \tau_{c3}$), average ($c_{3,r} = c_3$), and faster ($c_{3,r} = c_3 - \eta \cdot \tau_{c3}$) attenuating branches. Following the Miller III and Rice (1983) three point approximation of a Gaussian distribution, with $\eta \approx 1.732$ one can use logic tree branch weights of 0.167, 0.666, and 0.167, respectively for slower, average, and faster branches. Consequently, the ground-motion prediction is a weighted mixture of three Gaussian distributions $N(\ln(\mu), \sigma^2)$, where $\ln(\mu)$ is estimated from Eq. (2) for three values of $c_{3,r} \in (c_3, c_3 + \eta \cdot \tau_{c3}, c_3 - \eta \cdot \tau_{c3})$.

Within the context of an ergodic application, if a site has either the measured $V_{S30}$ or slope information available, but no site-specific ground-motion recordings, then the $\phi_{S2S}$ in the estimation of $\sigma = \sqrt{\tau_{L2L}^2 + \tau_0^2 + \phi_{S2S}^2 + \phi^2}$ can be replaced with $\phi_{S2S}^{V_{S30}}$ and $\phi_{S2S}^{slope}$, whilst appending the fixed-effects of Eq. (2) with Eqs. (6) and (7), respectively. Figure 14 shows the consequent reduction of $\sigma$ to $\sigma^{V_{S30}}$ and $\sigma^{slope}$ values, when using site-response proxies $V_{S30}$ and slope, respectively. Consequently, the ground-motion predictions follow the mixed Gaussian distribution $N(\ln(\mu) + SR^{V_{S30}}, \sigma^{V_{S30}}^2)$ or $N(\ln(\mu) + SR^{slope}, \sigma^{slope}^2)$, where $\ln(\mu)$ is estimated from Eq. (2) for three values of $c_{3,r} \in [c_3, c_3 + \eta \cdot \tau_{c3}, c_3 - \eta \cdot \tau_{c3}]$, and $SR^{V_{S30}}$ or $SR^{slope}$ are estimated from Eq. (6) or (7). The reduced aleatory variability from using a site-response proxy is beneficial until when enough ground-motion data can be collected at a site, and then $\delta S_{2S_s}$ for the new site can be estimated using equations provided in Kotha et al. (2017); Rodriguez-Marek et al. (2013); Sahakian et al. (2018); Villani and Abrahamson (2015), etc.

## 7.2 Region-specific application

For region-specific applications, the predictions can be upgraded with the region-specific anelastic attenuation and tectonic locality specific adjustments. Anelastic attenuation is regionalised by adjusting the generic coefficient $c_3$ in Eq. (4) with a region-specific value
$\delta c_{3,r}$, as in $c_{3,r} = c_3 + \delta c_{3,r}$, where region $r$ is decided by the location of the site. Since $\delta c_{3,r}$ are estimated from a smaller region-specific sample of ground-motion recordings, it should be treated as epistemically uncertain. For this purpose, standard-error on $\delta c_{3,r}$ are provided as well, and these are always smaller than $\tau_{c3}$. Treating the $SE(\delta c_{3,r})$ as uncertainty on mean of a normally distributed sample, the 95% confidence interval of $\delta c_{3,r}$ would be $\delta c_{3,r} \pm 1.6SE(\delta c_{3,r})$.

Similarly, depending on the tectonic locality of the earthquake, the GMM predictions can be further regionalised by adjusting $e_1$ in Eq. (2) to $e_{1,l} = e_1 + \delta L2L_l$, where $l$ identifies locality tectonic locality of the earthquake. Since $\delta L2L_l$ is estimated from a smaller locality-specific ground-motion sample, the 95% confidence interval is bounded by $\delta L2L_l \pm 1.6SE(\delta L2L_l)$. A reduction of up to 10% in $\sigma$ is achieved by dropping the $\tau_{L2L}$ from aleatoric variance, resulting in a smaller $\sigma_r = \sqrt{\tau_0^2 + \phi_{S2S}^2 + \phi^2}$, as shown in Fig. 16. Further reductions to $\sigma_r^{\text{slope}}$ are $\sigma_r^{V_{S30}}$ possible if the site-response scaling with slope and $V_{S30}$ (Eqs. 6 and 7) are used in predictions, respectively.

It is interesting to compare the total standard-deviation $\sigma$ (dashed line in Fig. 16) of the regionalised, $V_{S30}$ based K16 GMM with $\sigma_r^{V_{S30}}$ of this study. At $T \le 0.3$ s, the $\sigma_r^{V_{S30}}$ is about 6% larger than its K16 counterpart; while at $T > 0.3$ s, $\sigma_r^{V_{S30}}$ is as much as 14% smaller. Despite the enormous increase in data, it appears that regionalisation has helped in substantially curtailing the increase in apparent aleatory variability.

## 7.3 Region- and site-specific application

Partially non-ergodic region- and site-specific ground-motion predictions are possible for those sites with $\delta S2S_s$ provided with this GMM or for new sites with sufficient ground-motion data. $\delta S2S_s$ for the 1829 sites in the ESM dataset are provided, along with the $\delta c_{3,r}$

<figure><img src="image_20.png" /><figcaption>**Fig. 16** Reduction in total standard deviation estimates ($\sigma$) of the GMM for $T=0.01–8$ s (solid lines) from ergodic with and without site-response proxy, region-specific with and without site-response proxy, and region- and site-specific values. Standard deviation estimates when using site-response proxies are indicated by the annotations with corresponding superscript, i.e. ($\sigma^{V_{S30}}$) with $V_{S30}$, and to ($\sigma^{\text{slope}}$) with slope. Annotations with subscript $r$ correspond to the variances for regionalized predictions, i.e. from discounting region-to-region variability $\tau_{L2L}$ from $\sigma$. Values annotated with subscript $r$, $s$ correspond to region- and site-specific predictions discounting also the $\phi_{S2S}$. The black dashed line is the total standard-deviation (same as in Fig. 14) of the regionalized and $V_{S30}$ based Kotha et al. (2016) GMM, and compares with the $\sigma_r^{V_{S30}}$ of this study</figcaption></figure>
of their region, and $\delta L2L_l$ of nearby tectonic localities (earthquake sources). Since between-site and between-locality variabilities do not apply for site-specific predictions, the reduction in apparent aleatory variability is enormous, i.e. $\sigma_{r,s} = \sqrt{\tau_0^2 + \phi^2}$ is about 40% smaller than $\sigma = \sqrt{\tau_{L2}^2 + \tau_0^2 + \phi_{S2S}^2 + \phi^2}$, as shown in Fig. 16. However, the reduction in aleatory variability will be accompanied by additional epistemic uncertainty. In addition to those in region-specific predictions, uncertainty on (the mean) $\delta S2S_s$ should be accounted with $\pm 1.6 \cdot SE(\delta S2S_s)$. Region- and site-specific predictions therefore are a mixture of 27 Gaussian distributions $N(\ln(\mu_{r,s}), \sigma_{r,s}^2)$, where $\ln(\mu_{r,s})$ is estimated from Eq. (1) for three values of $c_{3,r} \in [c_3 + \delta c_{3,r}, c_3 + \delta c_{3,r} + 1.6SE(\delta c_{3,r}), c_3 + \delta c_{3,r} - 1.6SE(\delta c_{3,r})]$, three values of $e_{1,l,s} \in [e_{1,l} + \delta S2S_s, e_{1,l} + \delta S2S_s + 1.6SE(\delta S2S_s), e_{1,l} + \delta S2S_s - 1.6SE(\delta S2S_s)]$, wherein $e_{1,l} \in [e_1 + \delta L2L_l, e_1 + \delta L2L_l + 1.6SE(\delta L2L_l), e_1 + \delta L2L_l - 1.6SE(\delta L2L_l)]$.

In complement, a more practical application is presented in the companion study by Weatherill et al. (2020c), where a more thorough comparison of this GMM with contemporary models is provided. Weatherill et al. (2020c) also provides a more complete discussion on how the various epistemic uncertainties in this GMM can be handled in the logic tree framework of a PSHA.

## 7.4 Towards non-ergodic ground-motion predictions

For tectonic localities, attenuating regions, and sites with sufficient amount of recordings the epistemic uncertainty on the random-effect adjustments are negligible with respect to the random-effect and standard-deviations. Collecting more ground-motion recordings is principal in moving towards non-ergodic predictions. The benefits in resolving the ergodic assumption and progressing towards region- and site-specific in ground-motion prediction is demonstrated Fig. 17. In this plot, predictions for the M6.5 Norcia event of the central Italy sequence, which occurred on 30th October, 2016, are compared to the response spectra recorded at three sites covered by the Italian strong motion network (Gorini et al. 2010). These sites are identified by the network code IT in the ESM dataset: (1) permanent, free-field station LSS (Leonessa) with $V_{S30} = 1091$ m/s located 25 km from the event epicentre, (2) permanent, free-field station MVB (Marsciano Monte Vibiano) with $V_{S30} = 1046$ m/s located 65 km from the event epicentre and, (3) permanent, free-field station PSC (Pescasseroli) with $V_{S30} = 1000$ m/s located 110 km from the event epicentre. The three columns in Fig. 17 correspond to the three stations.

These event and stations are selected to demonstrate progressively (in Fig. 17) the impact of moving from ergodic prediction relying on $V_{S30}$ as site-response proxy (top row), through region-specific predictions (middle row) considering regional (Northern and central Apennines West) anelastic attenuation ($c_{3,r} = c_3 + \delta c_{3,r}$) and adjustment specific to the tectonic locality ($e_{1,l} = e_1 + \delta L2L_l$) containing the event (locality ID: "PTTC007"), to region- and site-specific predictions (bottom row) from an additional site-specific adjustment ($e_{1,l,s} = e_1 + \delta L2L_l + \delta S2S_s$). Both the median prediction and standard-deviation change in process, which is reflected by the width of the coloured ribbon in Fig. 17. The $\delta S2S_s(T = 0.01 - 8s)$ of these sites are estimated from 29, 15, and 20 records from predominantly small-moderate earthquakes (details in the Fig. 17 panels). A few comments on this figure:

1. The ergodic median predictions (central line) and one $\sigma^{V_{S30}}$ interval (ribbon) are systematically above the observed response spectra at the three rock sites, located at near
<figure><img src="image_21.png" /><figcaption>**Fig. 17** Comparison of ergodic (top row), region-specific (middle row), and sites-specific (bottom row) ground-motion median (solid colored line) and variability (colored ribbon) predictions for the M6.5 Norcia earthquake with the observed response spectra (solid red line) at three sites in Italy (column wise)</figcaption></figure>

(25 km), intermediate (65 km), and far (110 km) source distances. This is likely because the M6.5 Norcia event produced relatively weaker ground-motions compared to other large magnitude events recorded in Greece and Turkey—as quantified by their respective $\delta B_{e,l}^0$ values. Since the ergodic predictions consider all event, region, site, and record variabilities as aleatory, the $\sigma^{V_{S30}}$ (Fig. 16) is large yet not large enough to contain the M6.5 event observations within $\pm\sigma^{V_{S30}}$ boundaries.

2. Region-specific ground-motion predictions for these sites are achieved by adjusting the GMM with the $\delta c_{3,r}$ of the Northern and central Apennines (West), in which the sites are located; and the $\delta L2L_l$ of the tectonic locality PTTC007, in which the M6.5 event occurred (along with a few other prominent events and aftershocks). The epistemic uncertainties of these adjustments are relatively very small given the large number of recordings. In the middle row of Fig. 17, we notice the observed response spectra are closer to region-specific predictions than to the ergodic predictions. $\delta c_{3,r}$ and $\delta L2L_l$ of

<figure><img src="image_22.png" /></figure>
these regions are both lower than the pan-European average (which is zero), meaning the region attenuates short-period ground-motions faster and the events on average produce weaker ground-motions than elsewhere in pan-European region.

It is interesting to note that, at short distances (site IT_LSS) the $δc_{3,r}$ has no effect on region-specific predictions, and the shift is mostly from $δL2L_f$—and so is at intermediate distance (site IT_MVB). At far-source distances (site IT_PSC), the combined effect of $δc_{3,r}^V$ and $δL2L_f$ worked well to capture the observed response spectra within the narrower $±σ_r^{V30}$ (Fig. 16) range about the region-specific median.

3. Region- and site-specific predictions (site-specific in short) for the three sites are shown in the bottom row of Fig. 17. Along with the curves, details on the number of recordings, magnitude and distance ranges (1st and 3rd quantile) of the recordings are provided. The additional adjustment to the preceding region-specific predictions is through $e_{1,r,s} = e_1 + δL2L_f + δS2S_s$. While most of the data in estimating the $δS2S_s$ of these sites is from small-moderate sized earthquakes, the site-specific predictions fit quite well with the observations for the large M6.5 event. Since $δS2S_s$ are used to adjust the region-specific ground-motions, the $V_{S30}$ becomes irrelevant, and $φ_{S2S}^{V30}$ is dropped from $σ_r^{V30}$ resulting in a smaller $σ_{r,s}$.

In the above example, we demonstrated that applying region-specific adjustments noticeably improved the match between observations and predictions. The best agreement was clearly from using region- and site-specific adjustments. To substantiate this claim, we performed a 10-fold cross-validation exercise to verify if the introduction of the various random-effects into GMM functional form indeed improves its prediction capabilities. In doing so, we re-run the regression with three functional forms:

* Ergodic model with no regionalisation of anelastic attenuation, no localisation of event terms, and no site-specific adjustments, i.e. no random-effects $δc_{3,r}$, $δL2L_f$, $δS2S_s$

* Regional model with regionalisation of anelastic attenuation and tectonic localisation of events, leaving out the random-effect to capture site-to-site variability, i.e. only the random-effects $δc_{3,r}$ and $δL2L_f$, but no $δS2S_s$

* Site-specific model identical to the GMM presented here i.e., regionalisation of anelastic attenuation, localisation event terms, and the site-to-site random-effect

To perform the cross-validation, the dataset is split into 10 parts with non-overlapping events. Meaning, earthquakes (and records) are exclusive to their subsets and do not feature in any other subset. We perform regression of the three models on any nine subsets combined, and test the predictions on the tenth subset. Root-mean-squared-errors are estimated for each trial and then averaged over the ten trials. This exercise is repeated for IMs:PGA, PGV, SA($T$ = 0.1, 0.2, 0.5, 1, 2, 4 s).

Figure 18 shows a histogram of RMSE for each IM, from the three models. The reduction in RMSE from ergodic to region-specific GMMs is clear and most prominent at short periods. Towards longer periods, the improvement is less pronounced but is still substantial. This is because the regionalisation random-effects are focused on the capturing variabilities in short-period ground-motions. In Fig. 7, we notice the overall variability $σ$ peaks at short periods. In Figs. 8 and 10, the largest regional variability of anelastic attenuation $τ_{c3}$ and tectonic localities $τ_{L2L}$ is also at short periods $T ≤ 0.3$ s. Without regionalisation, along with the increased the short period ground-motion
<figure><img src="image_23.png" /><figcaption>**Fig. 18** 10-fold cross validation statistics comparing the predictive capabilities of ergodic, regionalized, and site-specific GMM versions at various periods</figcaption></figure>

prediction variability and reduced precision of an ergodic median, the predictive capability (measured as RMSE) of the GMM is reduced as well.

Across all the periods, i.e. the entire response spectra, the best predictive capabilities are those of region- and site-specific GMM. In Fig. 16, the largest reduction in $\sigma$ is achieved not from using $V_{S30}$ or slope as the site-response proxy, but from using $\delta S2S_s$ as the site-specific adjustment. However, it is unlikely that every site has sufficient ground-motion data to estimate its site-specific $\delta S2S_s$. In that case, alternative site-response proxies are sought to predict the $\delta S2S_s$, as in Kotha et al. (2018); Weatherill et al. (2020b). However, even in these studies, while the long period site-response could be partially explained using some geotechnical parameters, short-period site-response is much more variable—even among the so-called reference rock sites (Bard et al. 2019; Pilz et al. 2020).

# 8 Summary

In this study, we present the update of the Kotha et al. (2016) region- and site-specific GMM, using the recently published European Strong Motion dataset (ESM) by Lanzano et al. (2019a) and Bindi et al. (2018b). The update derived from a dataset 15 times larger than the older RESORCE dataset (Akkar et al. 2014b), with only minor changes to fixed-effects formula, and one additional random-effect grouping. Leveraging on the recently developed regionalisation datasets, TSUMAPS-NEAM by Basili et al. (2019) and the tectonic zonation defined in the purview of European Seismic Hazard Model 2020 project, we have regionalised the ground-motion data to capture spatial variation of anelastic attenuation and earthquake characteristics, respectively.

Due to the increased spatial and temporal diversity of the data, the random-effect and residual variances are larger than the older pan-European GMMs. This increased variance motivates the development of region-dependent ground-motion models since a larger ergodic $\sigma$ will severely affect the probabilistic seismic hazard and risk assessments. We propose making use of the attenuation region, tectonic locality, and site-specific terms (random-effects), along with their epistemic uncertainty, to develop partially non-ergodic ground-motion predictions with a reduced $\sigma$. The improved predictive power of the region- and site-specific version of the GMM is substantiated with a 10-fold cross-validation. All the random-effect values will be provided on request, along with the fixed-effects
coefficients, and the variance-covariance matrices needed to update the model in a Bayesian framework.

Only the customary random-effect and residual analyses are presented here, leaving out their physical evaluation for a follow-up study. A similar exercise is actually carried out on the Fourier amplitude version of the ESM dataset. Random-effects and residuals of the Fourier GMM are easier to associate to geophysical and geological parameters that are not involved in the GMM regression. If indeed the spatial variability of source, path, and site random-effects can be attributed to a globally available parameter, e.g. a crustal velocity map, we can attempt migrating the GMM developed from seismically active regions to the less active regions with sparse ground-motion recordings. The GMM is developed with the intention of perpetual updates, when and wherever new datasets become available. Further elaboration on random-effect analyses, GMM update and application procedures are ongoing, and will complement this study. Meanwhile, this new GMM has the widest applicability yet for the pan-European region; as it is derived from a dataset stretching from Pyrenees in the west to Iran in the east, Rhine-Graben in the north to Hellenic arc in the south, and constituted of manually processed ground-motion data shallow crustal earthquakes of $3 \le M_W \le 7.4$ recorded at $0 \le R_{JB} \le 545\text{km}$.

## 9 Data and resources

The European Strong Motion flatfile is available at https://esm.mi.ingv.it/flatfile-2018/ with persistent identifier PID: 11099/ESM_flatfile_2018. The analyses in this study have been performed in R software (Team 2013). In particular, we used the libraries *rlmm* (Koller 2016), *dplyr* (Wickham et al. 2019b), *ggplot2* (Wickham et al. 2019a), *ggmap* (Kahle et al. 2019), *viridis* (Garnier 2019), etc. The electronic supplement provides all the data used, derived, and presented in this study. We further encourage readers to contact the authors for more details, GMM implementation in OpenQuake™ (Pagani et al. 2014), regression scripts, etc.

**Acknowledgements** Open Access funding provided by Projekt DEAL. We are grateful to the swift and insightful reviews by Prof. Sinan Akkar and an anonymous reviewer, and Prof. John Douglas in his capacity as the handling editor of this manuscript at Bulletin of Earthquake Engineering. The contributions of the Sreeram Reddy Kotha (corresponding author) in this research are funded by the SIGMA2 consortium (EDF, CEA, PG&E, SwissNuclear, Orano, CEZ, CRIEPI) under Grant—2017–2021. The model development has benefitted immensely from feedbacks provided by Dr. Paolo Traversa, SIGMA2 scientific committee, and the collaborators in Horizon 2020 (Grant No. 730900) “Seismology and Earthquake Engineering Research Infrastructure Alliance for Europe (SERA)” project.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
References

Abrahamson N, Silva WJ (1997) Empirical response spectral attenuation relations for shallow crustal earthquakes. Seismol Res Lett 68:94–127

Abrahamson N, Youngs R (1992) A stable algorithm for regression analyses using the random effects model. Bull Seismol Soc Am 82:505–510

Abrahamson NA, Silva WJ, Kamai R (2014) Summary of the ASK14 ground motion relation for active crustal regions. Earthq Spectra 30:1025–1055. https://doi.org/10.1193/070913eqs198m

Akkar S, Sandikkaya M, Bommer J (2014a) Empirical ground-motion models for point-and extended-source crustal earthquake scenarios in Europe and the Middle East. Bull Earthq Eng 12:359–387

Akkar S et al (2014b) Reference database for seismic ground-motion in Europe (RESOURCE). Bull Earthq Eng 12:311–339. https://doi.org/10.1007/s10518-013-9506-8

Al Atik L, Abrahamson N, Bommer JJ, Scherbaum F, Cotton F, Kuehn N (2010) The variability of ground-motion prediction models and its components. Seismol Res Lett 81:794–801

Ancheta TD et al (2014) NGA-West2 database. Earthq Spectra 30:989–1005

Anderson JG, Brune JN (1999) Probabilistic seismic hazard analysis without the ergodic assumption. Seismol Res Lett 70:19–28

Atik LA, Youngs RR (2014) Epistemic uncertainty for NGA-West2 models. Earthq Spectra 30:1301–1318

Bard P, Bora SS, Hollender F et al (2019) Are the standard VS-Kappa host-to-target adjustments the only way to get consistent hard-rock ground motion prediction? Pure Appl Geophys. https://doi.org/10.1007/s00024-019-02173-9

Basili R et al (2019) NEAMTHM18 documentation: the making of the TSUMAPS-NEAM Tsunami Hazard Model 2018

Bates D, Mächler M, Bolker B, Walker S (2014) Fitting linear mixed-effects models using lme4. arXiv pre-print arXiv:14065823

Bazzurro P, Cornell CA (2004) Nonlinear soil-site effects in probabilistic seismic-hazard analysis. Bull Seismol Soc Am 94:2110–2123. https://doi.org/10.1785/0120030216

Bindi D, Kotha SR (2020) Spectral decomposition of the Engineering Strong Motion (ESM) flat file: regional attenuation, source scaling and Arias stress drop. Bull Earthq Eng 18:2581–2606. https://doi.org/10.1007/s10518-020-00796-1

Bindi D, Parolai S, Grosser H, Milkereit C, Karakisa S (2006) Crustal attenuation characteristics in north-western Turkey in the range from 1 to 10 Hz. Bull Seismol Soc Am 96:200–214

Bindi D, Massa M, Luzi L, Ameri G, Pacor F, Puglia R, Augliera P (2014) Pan-European ground-motion prediction equations for the average horizontal component of PGA, PGV, and 5%-damped PSA at spectral periods up to 3.0 s using the RESOURCE dataset. Bull Earthq Eng 12:391–430

Bindi D, Cotton F, Kotha SR, Bosse C, Stromeyer D, Grünthal G (2017) Application-driven ground motion prediction equation for seismic hazard assessments in non-cratonic moderate-seismicity areas. J Seismol 21:1201–1218. https://doi.org/10.1007/s10950-017-9661-5

Bindi D, Cotton F, Spallarossa D, Picozzi M, Rivalta E (2018a) Temporal variability of ground shaking and stress drop in Central Italy: a hint for fault healing? Bull Seismol Soc Am 108:1853–1863

Bindi D, Kotha S, Weatherill G et al (2018b) The pan-European engineering strong motion (ESM) flatfile: consistency check via residual analysis. Bull Earthq Eng 17:583–602. https://doi.org/10.1007/s10518-018-0466-x

Bindi D, Picozzi M, Spallarossa D, Cotton F, Kotha SR (2019) Impact of magnitude selection on aleatory variability associated with ground-motion prediction equations: part II—analysis of the between-event distribution in Central Italy. Bull Seismol Soc Am 109:251–262

Bohnhoff M, Martínez-Garzón P, Bulut F, Stierle E, Ben-Zion Y (2016) Maximum earthquake magnitudes along different sections of the North Anatolian fault zone. Tectonophysics 674:147–165

Bommer JJ, Abrahamson NA (2006) Why do modern probabilistic seismic-hazard analyses often lead to increased hazard estimates? Bull Seismol Soc Am 96:1967–1977

Bommer JJ, Douglas J, Scherbaum F, Cotton F, Bungum H, Fäh D (2010) On the selection of ground-motion prediction equations for seismic hazard analysis. Seismol Res Lett 81:783–793

Boore DM (2010) Orientation-independent, nongeometric-mean measures of seismic intensity from two horizontal components of motion. Bull Seismol Soc Am 100:1830–1835

Boore DM, Stewart JP, Seyhan E, Atkinson GM (2014) NGA-West2 equations for predicting PGA, PGV, and 5% damped PSA for shallow crustal earthquakes. Earthq Spectra 30:1057–1085

Bora SS, Cotton F, Scherbaum F, Edwards B, Traversa P (2017) Stochastic source, path and site attenuation parameters and associated variabilities for shallow crustal European earthquakes. Bull Earthq Eng 15:4531–4561. https://doi.org/10.1007/s10518-017-0167-x
Brooks C, Douglas J, Shipton Z (2020) Improving earthquake ground-motion predictions for the North Sea. *J Seismol.* https://doi.org/10.1007/s10950-020-09910-x

Brune JN (1970) Tectonic stress and the spectra of seismic shear waves from earthquakes. *J Geophys Res* 75:4997–5009

Campbell KW (1981) Near-source attenuation of peak horizontal acceleration. *Bull Seismol Soc Am* 71:2039–2070

Campbell KW, Bozorgnia Y (2014) NGA-West2 ground motion model for the average horizontal components of PGA, PGV, and 5% damped linear acceleration response spectra. *Earthq Spectra* 30:1087–1115. https://doi.org/10.1193/062913eqs175m

Chiou BS-J, Youngs RR (2014) Update of the Chiou and Youngs NGA model for the average horizontal component of peak ground motion and response spectra. *Earthq Spectra* 30:1117–1153

Cong L, Mitchell B (1998) Lg coda Q and its relation to the geology and tectonics of the Middle East. In: Mitchell BJ, Romanowicz B (eds) Q of the Earth: global, regional, and laboratory studies. *Pageoph topical volumes, Birkhäuser, Basel*. https://doi.org/10.1007/978-3-0348-8711-3_15

Crowley H et al (2019) Methods for estimating site effects in risk assessments, vol Deliverable 26.4, Final edn. *Seismology and Earthquake Engineering Research Infrastructure Alliance for Europe (SERA)*

Dawood HM, Rodriguez-Marek A (2013) A method for including path effects in ground-motion prediction equations: an example using the Mw 9.0 Tohoku earthquake aftershocks. *Bull Seismol Soc Am* 103:1360–1372

Dawood HM, Rodriguez-Marek A, Bayless J, Goulet C, Thompson E (2016) A flatfile for the KiK-net database processed using an automated protocol. *Earthq Spectra* 32:1281–1302

Derras B, Bard PY, Cotton F, Bekkouche A (2012) Adapting the neural network approach to PGA prediction: an example based on the KiK-net data. *Bull Seismol Soc Am* 102:1446–1461

Douglas J (2004) An investigation of analysis of variance as a tool for exploring regional differences in strong ground motions. *J Seismol* 8:485–496. https://doi.org/10.1007/s10950-004-3094-7

Douglas J (2018a) Calibrating the backbone approach for the development of earthquake ground motion models. Best practice in physics-based fault rupture models for seismic hazard assessment of nuclear installations: issues and challenges towards full seismic risk analysis

Douglas J (2018b) Capturing geographically-varying uncertainty in earthquake ground motion models or what we think we know may change. In: Pitilakis K (ed) Recent advances in earthquake engineering in Europe: 16th European conference on earthquake engineering—Thessaloniki. Springer, Cham, pp 153–181. https://doi.org/10.1007/978-3-319-75741-4_6

Douglas J, Edwards B (2016) Recent and future developments in earthquake ground motion estimation. *Earth Sci Rev* 160:203–219

Drouet S, Ameri G, Le Dortz K et al (2020) A probabilistic seismic hazard map for the metropolitan. France. *Bull Earthq Eng* 18:1865–1898. https://doi.org/10.1007/s10518-020-00790-7

Faccioli E, Paolucci R, Vanini M (2015) Evaluation of probabilistic site-specific seismic-hazard methods and associated uncertainties, with applications in the Po Plain, northern Italy. *Bull Seismol Soc Am* 105:2787–2807

Garnier S (2019) Viridis: default color maps from “ matplotlib”. 2018. https://github.com/sjmgarnier/ viridis. R package version 03 4:27

Giardini D et al (2018) Seismic hazard map of the Middle East. *Bull Earthq Eng* 16:3567–3570

Gorini A et al (2010) The Italian strong motion network. *Bull Earthq Eng* 8:1075–1090

Grad M, Tiira T, Group EW (2009) The Moho depth map of the European Plate. *Geophys J Int* 176:279–292

Grünthal G, Wahlström R (2012) The European-Mediterranean earthquake catalogue (EMEC) for the last millennium. *J Seismol* 16:535–570

Grünthal G, Stromeyer D, Bosse C, Cotton F, Bindi D (2018) The probabilistic seismic hazard assessment of Germany—version 2016, considering the range of epistemic uncertainties and aleatory variability. *Bull Earthq Eng* 16:4339–4395

Huber PJ (1992) Robust estimation of a location parameter. In: Kotz S, Johnson NL (eds) Breakthroughs in statistics. Springer series in statistics (perspectives in statistics). Springer, New York, NY. https://doi.org/10.1007/978-1-4612-4380-9_35

Jacoby WG (2000) Loess: a nonparametric, graphical tool for depicting relationships between variables. *Electoral Stud* 19:577–613

Kahle D, Wickham H, Kahle MD (2019) Package ‘ggmap’

Kale Ö, Akkar S, Ansari A, Hamzehloo H (2015) A ground-motion predictive model for Iran and Turkey for horizontal PGA, PGV, and 5% damped response spectrum: investigation of possible regional effects. *Bull Seismol Soc Am* 105:963–980
Kohrangi M, Kotha SR, Bazzurro P (2020) Impact of partially non-ergodic site-specific probabilistic seismic hazard on risk assessment of single buildings. *Earthq Spectra* (in-review)

Koller M (2016) robustlmm: an R package for robust estimation of linear mixed-effects models. *J Stat Softw* 75:1–24

Kotha SR, Bindi D, Cotton F (2016) Partially non-ergodic region specific GMPE for Europe and Middle-East. *Bull Earthq Eng* 14:1245–1263

Kotha SR, Bindi D, Cotton F (2017) From ergodic to region- and site-specific probabilistic seismic hazard assessment: method development and application at European and Middle Eastern sites. *Earthq Spectra* 33:1433–1453. https://doi.org/10.1193/081016EQS130M

Kotha SR, Cotton F, Bindi D (2018) A new approach to site classification: mixed-effects ground motion prediction equation with spectral clustering of site amplification functions. *Soil Dyn Earthq Eng.* https://doi.org/10.1016/j.soildyn.2018.01.051

Kotha SR, Cotton F, Bindi D (2019) Empirical models of shear-wave radiation pattern derived from large datasets of ground-shaking observations. *Sci Rep* 9:1–11. https://doi.org/10.1038/s4159 8-018-37524-4

Kowsari M, Halldorsson B, Hrafnkelsson B, Snæbjörnsson Jb, Jónsson S (2019) Calibration of ground motion models to Icelandic peak ground acceleration data using Bayesian Markov chain Monte Carlo simulation. *Bull Earthq Eng* 17:2841–2870

Kowsari M, Sonnemann T, Halldorsson B, Hrafnkelsson B, Snæbjörnsson Jb, Jónsson S (2020) Bayesian inference of empirical ground motion models to pseudo-spectral accelerations of south Iceland seismic zone earthquakes based on informative priors. *Soil Dyn Earthq Eng* 132:106075

Kuehn NM, Abrahamson NA (2017) The effect of uncertainty in predictor variables on the estimation of ground-motion prediction equations. *Bull Seismol Soc Am* 108:358–370. https://doi.org/10.1785/0120170166

Kuehn NM, Scherbaum F (2016) A partially non-ergodic ground-motion prediction equation for Europe and the Middle East. *Bull Earthq Eng* 14:2629–2642

Kühn NM, Scherbaum F (2015) Ground-motion prediction model building: a multilevel approach. *Bull Earthq Eng* 13:2481–2491

Landwehr N, Kuehn NM, Scheffer T, Abrahamson N (2016) A nonergodic ground-motion model for California with spatially varying coefficients. *Bull Seismol Soc Am* 106:2574–2583

Lanzano G, Sgobba S, Luzi L et al (2019a) The pan-European engineering strong motion (ESM) flatfile: compilation criteria and data statistics. *Bull Earthq Eng* 17:561–582. https://doi.org/10.1007/s1051 8-018-0480-z

Lanzano G, Luzi L, Pacor F, Felicetta C, Puglia R, Sgobba S, D'Amico M (2019b) A revised ground-motion prediction model for shallow crustal earthquakes in Italy. *Bull Seismol Soc Am* 109:525–540

Lu Y, Stehly L, Paul A, Group AW (2018) High-resolution surface wave tomography of the European crust and uppermost mantle from ambient seismic noise. *Geophys J Int* 214:1136–1150

Manighetti I, Campillo M, Bouley S, Cotton F (2007) Earthquake scaling, fault segmentation, and structural maturity. *Earth Planet Sci Lett* 253:429–438

Mayor J, Traversa P, Calvet M, Margerin L (2018) Tomography of crustal seismic attenuation in Metropolitan France: implications for seismicity analysis. *Bull Earthq Eng* 16:2195–2210

Miller AC III, Rice TR (1983) Discrete approximations of probability distributions. *Manag Sci* 29:352–362
Pacor F et al (2018) NESS1: a worldwide collection of strong-motion data to investigate near-source effects. *Seismol Res Lett* 89:2299–2313

Pagani M et al (2014) OpenQuake engine: an open hazard (and risk) software for the global earthquake model. *Seismol Res Lett* 85:692–702. https://doi.org/10.1785/0220130087

Pilz M, Cotton F, Zaccarelli R, Bindi D (2019) Capturing regional variations of hard-rock attenuation in Europe. *Bull Seismol Soc Am* 109:1401–1418. https://doi.org/10.1785/0120190023

Pilz M, Cotton F, Kotha SR (2020) Data-driven and machine learning identification of seismic reference stations in Europe. *Geophys J Int*. https://doi.org/10.1093/gji/gga199

Piña-Valdés J, Socquet A, Cotton F, Specht S (2018) Spatiotemporal variations of ground motion in Northern Chile before and after the 2014 M w 8.1 iquique megathrust event. *Bull Seismol Soc Am* 108:801–814

Radiguet M, Cotton F, Manighetti I, Campillo M, Douglas J (2009) Dependency of near-field ground motions on the structural maturity of the ruptured faults. *Bull Seismol Soc Am* 99:2572–2581

Rodríguez-Marek A et al (2013) A model for single-station standard deviation using data from various tectonic regions. *Bull Seismol Soc Am* 103:3149–3163

Sahakian V, Baltay A, Hanks T, Buehler J, Vernon F, Kilb D, Abrahamson N (2018) Decomposing leftovers: event, path, and site residuals for a small-magnitude Anza region GMPE. *Bull Seismol Soc Am* 108:2478–2492
Sahakian V, Baltay AS, Hanks TC, Buehler J, Vernon FL, Kilb D, Abrahamson NA (2019) Ground-motion residuals, path effects, and crustal properties: a pilot study in Southern California. *J Geophys Res* 124:5738–5753

Schmedes J, Archuleta RJ (2008) Near-source ground motion along strike-slip faults: insights into magnitude saturation of PGV and PGA. *Bull Seismol Soc Am* 98:2278–2290

Sedaghati F, Pezeshk S (2017) Partially nonergodic empirical ground-motion models for predicting horizontal and vertical PGV, PGA, and 5% damped linear acceleration response spectra using data from the Iranian Plateau. *Bull Seismol Soc Am* 107:934–948

Stafford PJ (2014) Crossed and nested mixed-effects approaches for enhanced model development and removal of the ergodic assumption in empirical ground-motion models. *Bull Seismol Soc Am* 104:702–719

Strasser FO, Abrahamson NA, Bommer JJ (2009) Sigma: issues, insights, and challenges. *Seismol Res Lett* 80:40–56

Team RC (2013) R foundation for statistical computing. Team RC, Vienna, p 3

Thompson E, Wald DJ, Worden C (2014) A VS30 map for California with geologic and topographic constraints. *Bull Seismol Soc Am* 104:2313–2321

Venables W, Riley B (2002) Modern applied statistics (Fourth S., editor). Springer, New York

Vilanova SP et al (2018) Developing a geologically based VS30 site-condition model for portugal: methodology and assessment of the performance of proxies developing a geologically based VS30 site-condition model for Portugal. *Bull Seismol Soc Am* 108:322–337

Villani M, Abrahamson NA (2015) Repeatable site and path effects on the ground-motion sigma based on empirical data from southern California and simulated waveforms from the CyberShake platform. *Bull Seismol Soc Am* 105:2681–2695. https://doi.org/10.1785/0120140359

Wald DJ, Allen TI (2007) Topographic slope as a proxy for seismic site conditions and amplification. *Bull Seismol Soc Am* 97:1379–1395

Weatherill G, Crowley H, Lemoine A, Roullé A, Tourlière B, Kotha SR, Cotton F (2020a) Modelling seismic site response at regional scale for the 2020 European Seismic Risk Model (ESRM20). *Bull Eng (in-preparation)*

Weatherill G, Kotha SR, Cotton F (2020b) Re-thinking site amplification in regional seismic risk assessment. *Earthq Spectra*. https://doi.org/10.1177/8755293019899956

Weatherill G, Kotha SR, Cotton F (2020c) A regionally-adaptable “Scaled-Backbone” ground motion logic tree for shallow seismicity in europe: application in the 2020 European seismic hazard model. *Bull Earthq Eng (in-review)*

Weatherill G, Kotha SR, Cotton F, Bindi D, Danciu L (2020d) Updated GMPE logic tree and rock/soil parameterisation for ESHM18. vol Deliverable 25.4. Seismology and Earthquake Engineering Research Infrastructure Alliance for Europe (SERA)

Wickham H, Chang W, Henry L, Pedersen T, Takahashi K, Wilke C, Woo K (2019a) R package ‘ggplot2’ v. 3.1. 1. Cran R

Wickham H, François R, Henry L, Müller K (2019b) dplyr: a grammar of data manipulation. R package version 0.8. 0.1. ed

Woessner J et al (2015) The 2013 European seismic hazard model: key components and results. *Bull Earthq Eng* 13:3553–3596. https://doi.org/10.1007/s10518-015-9795-1

**Publisher’s Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.