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

HAL
open science

# A regionally adaptable ground-motion model for fourier amplitude spectra of shallow crustal earthquakes in Europe

Sreeram Reddy Kotha, Dino Bindi, Fabrice Cotton

► To cite this version:

Sreeram Reddy Kotha, Dino Bindi, Fabrice Cotton. A regionally adaptable ground-motion model for fourier amplitude spectra of shallow crustal earthquakes in Europe. Bulletin of Earthquake Engineering, 2022, 20 (2), pp.711-740. 10.1007/s10518-021-01255-1. hal-04502005

HAL Id: hal-04502005

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

Submitted on 18 Sep 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Regionally Adaptable Ground-Motion Model for Fourier Amplitude Spectra of Shallow Crustal Earthquakes in Europe

**Sreeram Reddy Kotha** (✉ sreeram-reddy.kotha@univ-grenoble-alpes.fr)

Institut des Sciences de la Terre https://orcid.org/0000-0002-4874-3730

**Dino Bindi**

Helmholtz-Zentrum Potsdam GFZ: Deutsches Geoforschungszentrum Potsdam

**Fabrice Cotton**

University of Potsdam: Universitat Potsdam

## Research Article

**Keywords:** Ground-motion model, Fourier spectra, robust mixed-effects regression, regionally adaptable, seismic hazard and risk, Europe

**Posted Date:** April 15th, 2021

**DOI:** https://doi.org/10.21203/rs.3.rs-377090/v1

**License:** © This work is licensed under a Creative Commons Attribution 4.0 International License.
Read Full License

**Version of Record:** A version of this preprint was published at Bulletin of Earthquake Engineering on October 30th, 2021. See the published version at https://doi.org/10.1007/s10518-021-01255-1.
# A Regionally Adaptable Ground-Motion Model for Fourier Amplitude Spectra of Shallow Crustal Earthquakes in Europe

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

## Abstract

Typical ground-motion models predict the response spectral ordinates (GMM-SA), which are the damped responses of a suite of single-degree-of-freedom oscillators. Response spectra represent the response of an idealized structure to input ground-motion, but not the physics of the actual ground-motion. To complement the regionally adaptable GMM-SA of Kotha et al. (2020), we introduce here a model capable of predicting Fourier amplitudes (GMM-FA); developed from the same Engineering Strong Motion (ESM) dataset for pan-Europe. This GMM-FA reveals the very high variability of high frequency ground-motions, which are completely masked in a GMM-SA. By maintaining the development strategies of GMM-FA identical to that of the GMM-SA, we are able to evaluate the physical meaning of the spatial variability of anelastic attenuation and source characteristics. We find that a fully data-driven geospatial index, Activity Index (AIx), correlates well with the spatial variability of these physical effects. AIx is a fuzzy combination of seismicity and crustal parameters, and can be used to adapt the attenuation and source non-ergodicity of the GMM-FA to regions and tectonic localities sparsely sampled in ESM. While AIx, and a few other parameters we touch upon, may help understand the spatial variability of high frequency attenuation and source effects, the high frequency site-response variability - dominating the overall aleatory variance - is yet unresolvable. With the rapid increase in quantity and quality of ground-motion datasets, we call for an upgrade of regionalization techniques, site-characterisation, and a paradigm shift to Fourier ground-motion models to complement the traditional response spectra prediction models.

**Keywords:** Ground-motion model, Fourier spectra, robust mixed-effects regression, regionally adaptable, seismic hazard and risk, Europe

Sreeram Reddy Kotha

sreeram-reddy.kotha@univ-grenoble-alpes.fr

https://orcid.org/0000-0002-4874-3730

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

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

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

Typical ground-motion models (GMMs) used in seismic hazard and risk assessments predict the random distribution of ground-motion in terms of spectral amplitude (SA), i.e. damped response of an elastic single-degree-of-freedom (SDOF) oscillator with fundamental resonance period T. In Kotha et al. (2020) we introduced the regionally adaptable GMM derived from the Engineering Strong-Motion (ESM) dataset (Bindi et al. 2019a; Lanzano et al. 2018), developed specifically to predict the 5% damped response spectra SA(T) of shallow crustal earthquakes in Europe and Mediterranean regions (GMM-SA from hereon). GMM-SAs are quite useful in engineering applications because they predict the response of idealized SDOF structures to input ground-motion. However, SA(T) integrate characteristics of Fourier amplitudes over a range of frequencies FA(f), and hence do not exactly represent the actual ground-motion itself. Therefore, a GMM capable of predicting the actual ground-motion FA(f) (GMM-FA from hereon) can be a useful alternative to GMM-SA (Bora et al. 2014).

Recollecting from Bayless and Abrahamson (2019b), there are several advantages of GMM-FA over GMM-SA: (1) behavior of a GMM-FA is easier to constrain using seismological theory and models such as Brune (1970) and Boore et al. (2014a); (2) linear site-response remains linear at all frequencies and does not depend on the spectral content of input motions, which isn't the case with SA(T) according to Bora et al. (2016); (3) calibration of input parameters and methods for finite-fault simulations is more appropriate using GMM-FAs than using GMM-SAs; and (4) inter-frequency correlation of FA(f), developed from residuals of a GMM-FA, can facilitate stochastic simulation of realistic ground-motions (Bayless and Abrahamson 2019a; Stafford 2017). In addition, using the random vibration theory (RVT) framework of Bora et al. (2016) one can conveniently generate SA(T) from a GMM-FA.

Alongside the GMM-SA of Kotha et al. (2020), the ESM dataset provides also an opportunity to develop a GMM-FA for shallow crustal earthquakes originating in the pan-European region; extending from Pyrenness in the west to North-Eastern Anatolian fault system of Turkey. As it is, there are very few GMM-FA developed for this region – the most recent one being that of Bora et al. (2014) derived from an older pan-European strong ground-motion dataset RESORCE (Akkar et al. 2014). Bindi et al. (2019a) developed a GMM-FA from the ESM dataset, but it was not intended for engineering application but to perform a sanity check of data through residual analysis. Essentially, the last update of a compendium of GMMs published since 1964 by Douglas (2011) features only 18 models that predict Fourier Amplitude Spectra (FAS). Therefore, given their practical relevance, we consider it is worth introducing a new GMM-FA to better explore the FAS ground-motion prediction epistemic uncertainty.

Our intention for this new GMM-FA is in fact more than to contribute to the list of usable models. We proposed a regionally adaptable GMM-SA in Kotha et al. (2020), whose partially non-ergodic median predictions [of SA(T)] can be adjusted using the various attenuating region, tectonic locality, and site-specific random-effects. The statistical significance of random-effect groups and
overall improvement in GMM-SA performance has already been demonstrated in that paper. However, what remained unevaluated was the physical meaning of these random-effects, and their adaptability to regions with no data in ESM. Correlating the GMM-SA random-effects estimated in response spectral domain to geophysical, geological indices did not seem appropriate. For instance, the short-period SAs integrate features of moderate-high frequency FAS, making it difficult to interpret physically the quantified random-effects at $T < 0.5$s. Towards this end, in this study we introduce a regionally adaptable GMM-FA, correlate the regional adjustments to an independently derived geospatial index, and suggest using these correlation to adapt the model to regions with little to no data in the ESM dataset.

Our best candidate geospatial index is the Activity Index developed by Chen et al. (2018). Although we have tried a suite of geological and geophysical parameters (e.g. Moho depth, 1Hz coda Q maps, etc) to serve as proxy for our regional adjustments (the random-effects), the spatial coverage and resolution, and the epistemic uncertainty of these parameters made them less viable than the globally available Activity Index. More about this parameter will be discussed in relevant sections.

## 2 GROUND-MOTION DATASET

The ground-motion dataset we used is the Fourier amplitude version of the ESM dataset. The Fourier version of ESM does not contain data from Iran, and is therefore a few dozen records smaller than that used to derive the Kotha et al. (2020) GMM-SA. The ground-motion data selection criteria, regionalisation models, functional form, and regression methods are identical to that of the Kotha et al. (2020) GMM-SA - except for a few minor changes noted subsequently:

1) 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 923 events with hypocentral depth $0 < D \le 39km$ are located in regions with $14 \le \text{Moho depth} \le 49km$, as per the Moho map of Grad et al. (2009). Note that the same selection criteria resulted in 927 events in case of GMM-SA

2) Only those events with $\ge 3$ records in the dataset are used in regression. Also, wherever available, the harmonised $M_W$ estimates according to the European Mediterranean Earthquake Catalogue [EMEC] (Grünthal and Wahlström 2012) are preferred over the ESM default values. As with GMM-SA, preference of EMEC $M_W$ is to maintain consistency with the seismic source models developed for the new European Seismic Hazard Map (ESHM20)

3) We keep data from all 1622 sites in the dataset irrespective of whether their $V_{S30}$ measured from geotechnical 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). Note that the GMM-SA used data from 1829 sites, including 197 sites from Iran whose data is unavailable in the Fourier version of ESM
1

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} = 477\,km$. In comparison, the GMM-SA was derived from data extending up to $R_{JB} = 545\,km$

2

3

4

5) The only additional record selection criterion relevant to the GMM-FA is the low-pass frequency filter limit. For the GMM-SA regression of $SA(T)$, we chose only the records with a high-pass filter frequency of both horizontal components $f_{hp} \le 0.8/T$. While for the GMM-FA regression of $FA(f)$, we chose the records with high-pass frequency $f_{hp} \le 0.8f$ and low-pass frequency $f_{lp} \ge f/0.8$, of both horizontal components.

5

6

7

8

9

6) The GMM-FA is developed to predict the Effective Amplitude Spectra [EAS as defined by Pacific Earthquake Engineering Research Centre (Goulet et al. 2018)] of the two horizontal $FA(f)$ components, for 25 values of $f$ in the range $0.15\,Hz - 20\,Hz$. EAS is an orientation independent horizontal-component $FA(f)$ of ground-motion that can be used with RVT framework to estimate $SA(T)$ (Bora et al. 2019)

10

11

12

13

<figure><img src="image_3.png" /><figcaption>**Fig.1** Number of records, events, and stations used in regression at each frequency</figcaption></figure>

14

15

16

Following the above criteria, the number of records from $3.1 \le M_W \le 7.4$ events at distances $0 \le R_{JB} < 477\,km$ available for GMM-FA regression is shown in Fig.1

17

## 3 REGIONALISATION DATASETS

18

19

20

21

22

23

24

25

26

The GMM-SA of Kotha et al. (2020) is a regionally adaptable model, wherein the event and path effects are regionalised. It is capable of predicting $SA(T)$ accounting the regional differences in distance decay through $\Delta c_{3,r}(T) \approx N(0, \tau_{c3}(T))$ and average of localised source effects through $\Delta L2L_l(T) \approx N(0, \tau_{L2L}(T))$, and site-specific effects through $\Delta S2S_s(T) \approx N(0, \phi_{s2s}(T))$ random-effect adjustments to the generic GMM-SA median ($\ln(\mu)$). We use the same regionalisation models for GMM-FA, and therefore, regionalise the data identically. Regional differences in locality-specific event and region-specific anelastic attenuation characteristics are estimated as random-effects in a mixed-effects regression. The GMM-FA random-effect values can be used in three ways: 1) in region- and locality-specific predictions accounting their epistemic uncertainties, 2) ignored while reintroducing
1

2

3

4

5

6

their random-variances into total aleatoric variability of ergodic GMM-FA predictions – i.e. no region- and locality-specific adaption, 3) investigated for correlation with geophysical parameters or proxies in order to extend their usage to regions with no data in ESM. Kotha et al. (2020) demonstrated the first two applications, while also evaluating the robustness of region- and locality-specific GMM-SA through a 10-fold cross validation exercise. In this study we discuss the third approach, i.e. evaluating the physical meaning and adaptability of the GMM-FA random-effects.

7

**3.1 ATTENUATION REGIONALISATION**

8

9

10

11

12

13

14

In a drive towards regionalising ground-motion predictions, a few recent GMM-SAs identified and quantified between-region variability of anelastic attenuation, and attributed it to spatial variability of crustal characteristics. These GMM-SAs were regionalised based on geopolitical boundaries, in absence of a more relevant geophysical regionalisation. However, spatial variability of attenuation captured by the 1 Hz coda Q maps of smaller regions within geopolitical boundaries is quite high; as shown for Italy and Turkey by Cong and Mitchell (1998), for mainland France by Mayor et al. (2018), and for Europe in general by Pilz et al. (2019).

<figure><img src="image_4.png" /><figcaption>Fig. 2 Distribution of records in each of the 45 attenuation regionalisation polygons</figcaption></figure>

15

16

17

18

19

20

21

22

In this study, as in Kotha et al. (2020), we chose instead of geopolitical boundaries, a geological-geophysical regionalisation model developed by Basili et al. (2019) in the TSUMAPS-NEAM project. Fig. 2 shows the TSUMAPS-NEAM regionalisation and the number of records in each attenuating region, as determined by the location of the recording sites. Assuming that anelastic attenuation is a dominant phenomenon at far-source distances (e.g. $R_{JB} \ge 80\,\text{km}$) and in near-surface crustal layers, we let the location of the receiving site decide to which attenuating region a particular ground-motion record should be allotted to.
The 45 regions in this map (Fig.2) divide Turkey into around eight distinct attenuating regions, East and West Anatolia being the largest, while Italy is divided into 12 regions. Fig.2 also shows that while some regions have a few thousand records assigned (e.g. the central and northern Apennines), there are few regions with only a few dozen records (e.g. only 13 in the Rhine Graben). In mixed-effects regression terminology, the random-variance of the *group* of attenuation regions quantifies the observed variability of anelastic attenuation between regions, where regions are *levels* of the group. If the random variance of this group is non-zero, it means that the chosen regionalisation scheme is able to capture the spatial variability of anelastic attenuation. Complementing the Kotha et al. (2020), in this study, we evaluate physical meaning of the region-specific random-effects.

## 3.2 EVENT LOCALISATION

Alongside the traditional earthquake-to-earthquake variability captured by the between-event random-effect group $\Delta B_e \approx N(0, \tau_e)$, Kotha et al. (2020) introduced an additional random-effect to quantify locality-to-locality variability of earthquake characteristics; analogues to the location-to-location variability $\Delta L2L_l \approx N(0, \tau_{L2L})$, of Al Atik et al. (2010). We assign the 923 shallow crustal events in the ESM dataset to 55 seismotectonic zones defined in the European Seismic Hazard Model 2020 (ESHM20). The seismotectonic source zonation (called "TECTO") was designed to distinguish regional tectonic features influencing the generation of crustal earthquakes, but not the very localised features. 'TECTO' was thus a good model for earthquake source localisation in Kotha et al. (2020), and distinct from the TSUMAPS-NEAM of Basili et al. (2019) used for regionalisation of anelastic attenuation.

<figure><img src="image_5.png" /><figcaption>**Fig.3** Distribution of records from shallow crustal events in each of the 55 event localisation polygons (tectonic zones)</figcaption></figure>

As in Kotha et al. (2020), this random-effects group is the between- locality $\Delta L2L_l \approx N(0, \tau_{L2L})$, wherein the levels are tectonic localities (locality *l*) containing their assigned ESM crustal
events and records. Fig.3 shows the distribution of records from active shallow crustal events within the event localisation polygons. As with the attenuating region group, this group is introduced to quantify the earthquake locality-to-locality variability of ground-motion in the dataset through the between-locality variance ($\tau_{L2L}^2$), which if zero indicates that the regionalisation model is not able to identify any significant spatial variability of earthquake characteristics.

## 3.3 ACTIVITY INDEX

Chen et al. (2018) introduced a fully data-driven global tectonic regionalisation model for selection of ground-motion models in seismic hazard applications. Based on a fuzzy logic workflow, they have rendered a regular grid with a spacing of 0.5°, wherein each cell is assigned a probability of being an active tectonic region – the Activity Index (AIx), as shown in Fig.4. AIx in a grid cell is calculated from the following fuzzy rules:

1) **High AIx** – IF [seismic moment rate density] is ‘High’, AND [1Hz Lg coda Q] is ‘Low’, AND [S-wave velocity variation at 175km] is ‘Low’, THEN the region is [‘Active’]

2) **Low AIx** – IF [seismic moment rate density] is ‘Low’, AND [1Hz Lg coda Q] is ‘High’ AND [S-wave velocity variation at 175km] is ‘High’, THEN the region is [‘Stable’]

AIx is derived as a combination of seismic moment rate density (Weatherill et al. 2016), 1Hz Lg coda Q (Mitchell et al. 2008), and shear-wave velocity variation at 175km (Ritsema et al. 2011). We use this dataset to evaluate the various regional variabilities quantified as the GMM random-effects.

<figure><img src="image_6.png" /><figcaption>**Fig.4** Activity Index map of pan-European region by Chen et al. (2018)</figcaption></figure>

## 4 FUNCTIONAL FORM

A mixed-effects GMM is composed of fixed-effects and random-effects. Fixed-effects part of the GMM is the continuous function built as a combination of predictor variables, which in this case are event
magnitude $M_W$ and event-to-site distance metric $R_{JB}$ (in km). As suggested by Fukushima (1996), and later in Bora et al. (2017) and Bora et al. (2019), the fixed-effects component of our GMM-FA remains the same as GMM-SA presented in Kotha et al. (2020). However, although recent GMM-FA models from NGA-West2 dataset (Ancheta et al. 2014) by Bayless and Abrahamson (2019b) and Bora et al. (2019) retained a magnitude-dependent geometrical spreading component, following the suggestion by Cotton et al. (2008) we have dropped this term from our GMM-FA; which in the Kotha et al. (2020) GMM-SA was the coefficient $c_2$.

$$ \ln(\mu) = e_1 + f_{R,g}(M_W, R_{JB}) + f_{R,a}(R_{JB}) + f_M(M_W) + \delta L2L_l + \delta B_{e,l}^0 + \delta S2S_s + \varepsilon \quad (1) $$

$$ f_{R,g} = c_1 \cdot \ln^2 \left( \frac{R_{JB}^2 + h_D^2}{R_{ref}^2 + h_D^2} \right) \quad (2) $$

$$ 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 (3) $$

$$ 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 (4) $$

In the mixed-effects formulation of equation (1), $f_{Rg}$, $f_{Ra}$, $f_M$ are the fixed-effects components; wherein $e_1$, $b_1$, $b_2$, $b_3$, $c_1$, $c_3$ are the fixed-effects regression coefficients (eq. 2 and 3). Regarding the random-effects components:

1) Between-region random-effects follow a Gaussian distribution $\Delta c_{3,r} = N(0, \tau_{c3})$. $\tau_{c3}$ quantifies the between-region variability of anelastic attenuations across the attenuating regions $r$ introduced earlier in Fig.2. Alongside a generic region-independent $c_3$, region-specific adjustments $\delta c_{3,r}$ are estimated as random-effects for the 45 attenuation regions. For region-specific GMM-FA predictions we replace the apparent anelastic attenuation term $c_3$ in equation (3) with a $c_{3,r} = c_3 + \delta c_{3,r}$

2) Between-locality variability of observed ground-motions are captured by the random-effect $\Delta L2L = N(0, \tau_{L2L})$. Locality-specific adjustments $\delta L2L_l$, estimated for the 55 localities $l$ shown in Fig.3, can be used to make locality-specific ground-motion predictions by replacing the offset $e_1$ in equation (1) with $e_{1,l} = e_1 + \delta L2L_l$

3) Event-to-event variability in this GMM is the traditional between-event random-effects $\Delta B_e = N(0, \tau_e)$ filtered for between-locality variability $\Delta L2L_l = N(0, \tau_{L2L})$, and is now captured by the $\Delta B_{e,l}^0 = N(0, \tau_0)$; where, for an event $e$ located in tectonic locality $l$, the event-specific term
can be seen as $\delta B_{e,l}^0 \approx \delta B_e - \delta L2L_l$. $\tau_0$ is the generic event-to-event variability corrected for
locality-to-locality variability, and does not vary with locality $l$

4) Site-to-site response variability is captured by the site-specific random-effects $\Delta S2S = N(0, \phi_{S2S})$. The potential of $\delta S2S_s$ in site-specific GMM-SAs is well-known, and are useful in studying regional differences in site-response scaling with $V_{S30}$ (time-averaged shear-wave velocity in 30m top-soil) as in Kotha et al. (2016), topographic slope and geology as in Weatherill et al. (2020a) or a variety of site parameters (Kotha et al. 2018; Pilz and Cotton 2019; Zhu et al. 2020), and even in site-specific seismic risk assessment (Kohrangi et al. 2020). Site-specific GMM-FA adjustments $\delta S2S_s$, estimated for the 1622 sites $s$ in the dataset Fig.3, can be used to make site-specific ground-motion predictions by replacing the offset $e_1$ in equation (1) with $e_{1,s} = e_1 + \delta S2S_s$

5) The left-over residuals $\epsilon = N(0, \phi)$ 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 at near-source distances (Kotha et al. 2019), anelastic attenuation effects in far-source distances (Sahakian et al. 2019), non-linear soil response (Loviknes et al. 2021), etc.

In all, this GMM has four random-effect groups. With this configuration of mixed-effects GMM, we run a robust linear mixed-effects regressions [rlmm by Koller (2016)] independently for the 25 EAS IMs $f = 0.15 - 20Hz$. Along with the fixed-effect regression coefficients and random-effect adjustments, we produced also the fixed-effects variance-covariance matrices needed to estimate the GMM epistemic uncertainty (Atik and Youngs 2014; Bindi et al. 2017a) and to update the GMM coefficients in a Bayesian framework (Kowsari et al. 2019). These resources can be disseminated on request.

# 5 FIXED-EFFECTS

The regression results comprise of fixed-effect coefficients and their covariance matrices, random-effect values and their rlmm weights and standard errors, residuals, and variances. Similar to the GMM-SA, outlier events, localities, regions, stations, and recordings are flagged as well. Fig.5 presents the predicted Fourier amplitude spectra for a few scenarios of interest. In this plot, we show fixed-effects prediction of FAS for:

i. [M4, M5.5, M7] implying events with $M_W = 4, 5.5, 7$

ii. [Shallow, Intermediate, Deep] to illustrate the effect of $h_D = 4, 8, 12km$ for events with depth $D < 10km, 10km \le D < 20km, 20km \le D$, respectively

iii. [Slower, Average, Faster] attenuating regions to illustrate the effect of replacing $c_3$ in equation (3) with $c_{3,r} = c_3 + \delta c_{3,r}$, with $\delta c_{3,r} = 0, +\tau_{c3}, -\tau_{c3}$, respectively
<figure><img src="image_7.png" /><figcaption>**Fig.5.** Predicted Fourier spectra of the GMM-FA for various scenarios, differentiated by color for Shallow ($D < 10km$), Intermediate ($10km \le D < 20km$), and Deep ($20km \le D$) events in the left panel. In the right panel are the Fourier spectra differentiated by regional anelastic attenuation: Slower and Faster than pan-European Average</figcaption></figure>

In the left panel, showing the near-source predictions, we notice that the depth-dependence has a non-negligible effect on the amplitudes at distances $R \le 5km$. For the (M7, 5km) scenario, the epistemic uncertainty (orange ribbon) on the median is wide enough to cover the variation with depth. A large part of this epistemic uncertainty is from the lack of near-source data from large magnitude events, which is necessary to constrain the magnitude-scaling $f_M(M_W)$ component of the GMM at $M_W \ge M_h = 5.7$.

In the right panel, FAS predictions at far-source distances are shown. Evidently, the $\delta c_{3,r} = 0, +\tau_{c3}, -\tau_{c3}$ adjustments become active at $f \ge 2Hz$ in far-source predictions. At lower frequencies, the differences are much smaller – which is to be expected since coefficient $c_3$ is meant to capture the (apparent) anelastic attenuation of moderate-to-high frequency ground-motions.

Along with the depth and anelastic attenuation dependencies, we notice that with increasing magnitude the Fourier spectra become flatter in low-moderate frequency range, as the apparent corner-frequency shifts towards lower frequencies for larger magnitude events. At near-source distances (left panel), the spectra decay rapidly beyond $f \ge 10Hz$, while at far source distances (right panel) this behavior is observed earlier at $f \ge 5Hz$; most likely from the enhanced anelastic decay of high frequency ground-motions.

# 6 RANDOM-EFFECTS AND RESIDUALS

Prior to discussing the random-effects and residuals in following sections, there are few inferences derived from the random-effect and residual variances in Fig.6. In this plot, $\phi_{S2S}, \tau_{L2L}, \tau_0, \tau_{C3}$ are the random-effect standard deviations of between-site, between-locality (tectonic zones), between-event (after between-locality correction), and between-region (anelastic attenuation) random-variables, respectively. $\phi$ is the residual standard deviation. The *total-sigma* of the ergodic version of GMM-FA,
1 $\sigma = \sqrt{\phi_{S2S}^2 + \tau_{L2L}^2 + \tau_0^2 + \phi^2}$, does not include $\tau_{c3}$ because it is intended for use as an epistemic uncertainty on the anelastic attenuation term $c_3$ of the GMM (Douglas 2018; Weatherill et al. 2020c).

<figure><img src="image_8.png" /><figcaption>**Fig.6 Random-effect, residual, and total variance estimates of the GMM for $f = 0.15 – 20Hz$ (solid lines), and $\tau_e = \sqrt{\tau_0^2 + \tau_{L2L}^2}$ (dashed line)**</figcaption></figure>

All the random-variances are comparable in size, which means the random-effect groups are meaningful for this dataset. The largest variability, however, is the site-to-site response variability - captured by the between-site standard deviation $\phi_{S2S}$. Increasing monotonically at $f \ge 3Hz$, $\phi_{S2S}$ suggests that site-response (in the dataset) is highly variable at moderate-high frequencies. Such rapid increase of $\phi_{S2S}$ towards higher frequencies has been reported for a various datasets (Bayless and Abrahamson 2019b; Bindi et al. 2017b; Bindi et al. 2021). For instance, ground-motion amplification at a site with $\delta S2S_s(f = 10Hz) = 1.5 * \phi_{S2S}(f = 10Hz)$ is 20 times larger than that at a site with $\delta S2S_s(f = 10Hz) = -1.5 * \phi_{S2S}(f = 10Hz)$. The large variability in site-response and the consequently large $\phi_{S2S}$ (the largest contributor to $\sigma$) suggests that partially non-ergodic site-specific ground-motion predictions may soon become indispensable (Faccioli et al. 2015; Kotha et al. 2017; Rodriguez - Marek et al. 2013).

The next largest random variance is that of between-event variability quantified into $\tau_0$. Note that a part of the overall between-event (spatiotemporal) variability $\tau_e$ is quantified in to the between-locality variability $\tau_{L2L}$; as in $\tau_e = \sqrt{\tau_0^2 + \tau_{L2L}^2}$. Apparently, variability of event-specific effects is the highest at $f \le 0.5Hz$. Seismic moment and moment-magnitude are the event-specific parameters estimated at these frequencies. For a GMM with $M_W$ (from EMEC catalog) as an explanatory variable, such large between-event variability at $f \le 0.5Hz$ suggests large differences in observed ground-motions between events of identical $M_W$. The most likely cause, to our understanding, is errors in $M_W$ in the dataset. A few studies (Holmgren and Atkinson 2018; Kuehn and Abrahamson 2017) demonstrated that $M_W$ uncertainty is a contributor to the between-event variability at long period spectral accelerations, which are analogues to low frequency Fourier amplitudes. Beyond $f \ge 1Hz$ the $\tau_0$ values are lower and remain almost constant.
The counter-part of $\tau_0$ is the between-locality variability $\tau_{L2L}$, which captures the average event-to-event variability when events are localized into seismotectonic zones ('TECTO'). $\tau_{L2L}$ values are much smaller than $\tau_0$ at $f < 1$Hz, and increase monotonically above $\tau_0$ at $f \ge 5$Hz. Assuming the errors in $M_W$ are contained in $\tau_0$, we discuss the physical meaning of $\tau_{L2L}$ in subsequent sections.

Region-to-region variability of anelastic attenuation is quantified into $\tau_{c3}$. Only the high frequency ground-motions are attenuated exponentially with distance. Therefore, $\tau_{c3}$ increases towards high frequencies in Fig.6. The residual standard deviation $\phi$, corrected for all parametric fixed-effects, non-parametric region, locality, and site-specific random-effects, remains almost constant across the frequency range.

Fig.6 illustrates the significance of the chosen random-effect groups, and the frequency dependence of their random-variances. It is necessary to validate the random-effects by correlating them to a physical parameter or an index. Since the random-effects are estimated only for the regions, localities, and sites (*levels* – as explained in Regionalisation Datasets) whose ground-motion data is used in GMM regression, new levels with no ground-motion data cannot benefit from the non-ergodic level-specific GMM predictions. However, correlating random-effects to a physical parameter or a geospatial index may allow, in a limited way, predicting region, locality, and site-specific adjustments for new levels outside the regressed dataset through the correlated parameter.

## 6.1 ANELASTIC ATTENUATION VARIABILITY

Anelastic attenuation of high frequency ground-motions comes into play at intermediate-far source distances (e.g. $R \ge 80$km). The coefficient $c_3$ in GMM median captures the average rate of exponential decay of ground-motion, while $c_1$ captures the linear decay. Substantial correlation between $c_3$ and $c_1$ estimates are to be expected because they together model the decay of ground-motions with distance (Abrahamson et al. 2014; Boore et al. 2014b; Campbell and Bozorgnia 2014). Therefore, it is more appropriate to refer to $c_3$ as a coefficient for *apparent* anelastic attenuation. $\delta c_{3,r}$ is meant to capture regional variability of this exponential decay. Fig.7 shows the frequency dependence of $\delta c_{3,r}$ for the 45 regions identified with subscript $r$. The region-to-region variability is largest at $f \ge 5$Hz, reflecting the large $\tau_{c3}$ in the Fig.6.
<figure><img src="image_9.png" /><figcaption>**Fig.7** $\delta c_{3,r}$ for $f = 0.15 – 20Hz$. Each line corresponds to one of the 45 attenuation regions, with colors indicating their weight in robust regression. Overlain red curves mark the $\pm\tau_{c3}$ values. Regions with $\delta c_{3,r}(f)$ beyond $\pm1.345\tau_{c3}(f)$ are given a lower than a unit weight</figcaption></figure>

Fig.8 maps the regional variability of $\delta c_{3,r}$ at $f \approx 0.3, 1, 3, 10Hz$ in the pan-European region; wherein the red colored polygons ($\delta c_{3,r} > 0$) represent regions with anelastic attenuation slower/weaker than the pan-European average, and the blue colored polygons with $\delta c_{3,r} < 0$ are those regions attenuating faster/stronger than average. Sparsely sampled regions with large epistemic uncertainty on their $\delta c_{3,r}$ are faintly colored. A few interesting remarks on these maps:

1) Weatherill et al. (2020c) showed that regions with similar attenuation characteristics, as estimated in Kotha et al. (2020), are spatial clustered. Apparently, high seismicity regions (e.g. Italy and Greece), perhaps with their highly fractured top crustal layers, attenuate the high frequency ground-motions faster/stronger than lower seismicity regions (e.g. central Europe)

2) The most well-sampled regions are located in central Italy; Northern and central Apennines W (West) with 5505 records, and Northern and central Apennines E (East) with 3199 records. While the difference is negligible at low frequencies ($f = 0.348Hz$ top-left panel), at $f \ge 0.991Hz$ (clockwise from top-right to bottom-left panels) the difference in $\delta c_{3,r}$ between these two regions is 0.2; which roughly translates to 10% larger ground-motion predictions at $R = 80km$ towards east than to west. In addition to attenuating the high frequency ground-motions faster than the pan-European average, there appears to be a frequency dependent contrast between these adjacent regions.

Recent 3D shear wave velocity ($V_s$) maps produced by Lu et al. (2018) reveal a relatively higher $V_s$ in the Apennines (west) compared to the Adriatic basin (east) at 10km and 30km crustal depths. If we assume that a higher $V_s$ means a higher crustal quality factor, these $V_s$ maps would imply a more efficient propagation in the crust, which translates to slower decay towards east compared to the west. However, the spatial trends of $\delta c_{3,r}$ indicate the contrary in this region.

Coincidentally, Apennines and region towards its west were reported to have shallower Moho (Grad et al. 2009) and significantly higher crustal temperatures due to submarine volcanic activity in the Tyrrhenian sea (Diaferia et al. 2019), compared to the Adriatic sea (east of Apennines); the latter phenomenon presumably better correlated to stronger attenuation of
mechanical waves than shear-wave velocity alone. Inopportunely, the large uncertainties in estimates of these parameters (Moho depth and thermal gradient) at pan-European scale did not motivate a more quantitative evaluation. Therefore, we prefered a globally available geospatial index that is designed specifically for regionalisation of GMMs, i.e. the Activity Index in this study

3) The region in Aegean Sea attenuating fastest the high frequency ground-motions ($f = 9.903Hz$ in bottom-left panel) is the Gulf of Corinth, where the sites appear to have received highly attenuated ground-motions from traversing the Aegean volcanic arc. This observation makes a strong case for using lithosphere temperature maps to guide regionalisation of ground-motion anelastic attenuation

<figure><img src="image_10.png" /><figcaption>**Fig.8** $\delta c_{3,r}$ variation across the 45 regions in ESM dataset for $f = 0.3, 1, 3, 10Hz$ (clockwise from top-left to bottom-left panels). Blue polygons locate regions with anelastic attenuation faster than the pan-European average, red polygons locate regions with slower attenuation, and grey polygons are regions close to the average. Regions with fewer ground-motion observations, thereby larger epistemic error on $c_{3,r}$, are more transparent and appear fainter. Note that the color scale is the same in all the maps in order to emphasize the frequency dependence of spatial variability of anelastic attenuation</figcaption></figure>

4) Contrast in high frequency attenuation ($f = 9.903Hz$ in bottom-left panel) is observed around the Alps regions, while minor differences are noticeable between west, north, and central Anatolia. Interestingly though, the regions attenuating the slowest are the Pyrenees, the western-Balkans, and the region north of Black sea to the east (east side of the maps here). Although not conclusive, these regions of very slow anelastic decay coincide quite well with the regions with deepest Moho (Grad et al. 2009); suggesting Moho depth maps could assist in regionalisation of ground-motion attenuation (see supplement Fig.S1)

We note that changing the resolution or geometry of the regions may change the estimated spatial variability and values of $\delta c_{3,r}$ as well. However, given the current configuration, we seek physical features that may correlate (qualitatively) with $\delta c_{3,r}$. A recent study by Sahakian et al. (2019) using a large data set of small-magnitude earthquakes in Southern California suggested that crustal
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

shear-wave velocity $V_s$ is only weakly correlated to anelastic attenuation, which also seems to be the case here with Central Italy (observation #2 above). Regional variability of anelastic attenuation may in fact be a combination of regional variability of crustal shear-wave velocity, crustal thermal gradient influencing the rigidity modulus of the crust, Moho depth, and other parameters that, at the time of this study, were not uniformly mappable across the pan-European region.

<figure><img src="image_11.png" /><figcaption>**Fig.9** $\delta c_{3,r}$ of the 45 regions versus Activity Index at site locations within each region, for $f \approx 0.3, 1, 3, 10$Hz (clockwise from top-left to bottom-left). The blue lines are loess fits between the two parameters. Marker colors indicate the weight assigned to $\delta c_{3,r}$ of each region in the rmml regression</figcaption></figure>

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

Activity Index (AIx) is a unique *composite* parameter that we preferred in this study. Activity Index is a data-driven continuous parameter inferred from a fuzzy combination of shear-wave velocity, seismic moment rate density, and crustal quality factors across the globe. A 0.5° gridded map of AIx was generated by Chen et al. (2018) for the sole purpose of regionalising GMMs or selecting suitable GMMs for a region with no region-specific ground-motion data. In this study, we extracted the AIx for every site location in the ESM dataset. A region with *n* sites will therefore have *n* values of AIx, which can serve as an epistemic uncertainty on the region-specific AIx. Fig.9 shows the loess fit [non-parametric moving average by Jacoby (2000)] between the $\delta c_{3,r}$ of the 45 regions and the AIx of sites located within each region. A strong negative correlation is evident at moderate-high frequencies (bottom panels), where regional variability $\tau_{c3}$ is the largest.

The negative correlation between $\delta c_{3,r}$ and AIx in Fig.9 suggests that regions with high seismic moment rate density, low shear-wave velocity, and low 1Hz coda Q (therefore, high AIx) attenuate significantly faster than regions more likely to be cratonic (low AIx). Chen et al. (2018) indicate that the regional variability of AIx is dominated by regional variability of seismic moment rate density in active crustal
regions ($AI \ge 0.7$), and to that of shear-wave velocity and 1Hz code Q in relatively stable regions ($AI < 0.7$). ESM dataset contains sites located in regions with $0.4 \le AI$ as seen in Fig. 9. The smooth transition of $\delta c_{3,r}$ between stable cratonic ($0.4 \le AI < 0.7$) to the more seismically active regions ($0.7 \le AI \le 1$) is an indication that it is a physically meaningful random-effect. However, these regionalisation models are different in nature: Activity Index is fully data-driven and gridded, while the regionalisation used in this GMM regression is based on expert elicitation and polygonised. In that sense, although there is a decent agreement, neither of the models may sufficiently replace the other.

<figure><img src="image_12.png" /><figcaption>**Fig. 10** $\delta c_{3,r}$ at $f \approx 10Hz$ of the 45 regions versus the range of Activity Index at site locations within each region</figcaption></figure>

Fig. 10 shows the ranges of AIx within each of the attenuation region. The regions are ordered in decreasing order of $\delta c_{3,r}$ ($f \approx 10Hz$) from top to bottom. This figure is to illustrate the exclusivity of the two GMM regionalisation models. For instance, the two best-sampled regions, Northern and central Apennines W (West) and Northern and central Apennines E (East), despite their $\delta c_{3,r}$ ($f \approx 10Hz$) values differing by 0.2 still have significant overlap of AIx ranges i.e. $0.68 \le AI \le 0.88$ and $0.62 \le AI \le 0.84$, respectively. Meaning, data-driven AIx by itself may not resolve the differences between these two adjacent regions with contrasting attenuation characteristics as efficiently as the more subjective Basili et al. (2019) regionalisation model. In lieu of more refined and unified regionalisation models, we foresee using both the models in correspondence to explain and predict attenuation characteristics ($\delta c_{3,r}$) for regions outside the ESM dataset. A study in this direction is anticipated with the recently published RESIF-RAP dataset of ground-motions recorded by the French accelerometric network (Péquegnat et al. 2008) by Traversa et al. (2020), and other low-moderate seismicity regions whose data is not integrated into ESM.
## 6.2 EVENT LOCALITY VARIABILITY

Source variability is divided into two components in this GMM-FA: variability across tectonic localities producing the earthquakes $\Delta L2L \approx N(0, \tau_{L2L})$, and the locality corrected between-event variability $\Delta B_{e,l}^0 \approx N(0, \tau_0)$. Since events are exclusively nested in their respective tectonic localities, $\delta L2L_l$ quantifies the *average* of the nested events' ground-motion characteristics as $\delta B_{e,l}^0 \approx \delta B_e - \delta L2L_l$; with notations from Al Atik et al. (2010). Fig.6 illustrates the frequency dependence of $\tau_{L2L}$ and $\tau_0$. It appears that the two random-variances capture disjoint frequency dependent earthquake characteristics; where $\tau_0 > \tau_{L2L}$ at low-moderate frequencies, and vice versa at high frequencies.

<figure><img src="image_13.png" /><figcaption>**Fig.11** $\delta L2L_l$ for $f = 0.15 - 20Hz$. Each line corresponds to one of the 55 tectonic localities, with colors indicating their weight in rlmm regression. Overlain red curves mark the $\pm\tau_{L2L}$ values. Localities with $\delta L2L_l(f)$ beyond $\pm1.345\tau_{L2L}(f)$ are given a weight lower than one</figcaption></figure>

Fig.11 shows the $\delta L2L_l(f = 0.15 - 20Hz)$ of the 55 tectonic localities in ESM dataset. Resembling Fig.6, the scatter of $\delta L2L_l$ significantly increases at $f \ge 5Hz$ in Fig.11. The epistemic uncertainty of $\delta L2L_l$, i.e. the standard-error $SE(\delta L2L_l)$, are generally smaller than $\tau_{L2L}$. In this regard, dropping $\tau_{L2L}$ from the aleatoric variability ($\sigma$) and using instead the $\delta L2L_l \pm SE(\delta L2L_l)$ adjustments to regionalise the GMM predictions is recommended. A database of $\delta L2L_l$, $SE(\delta L2L_l)$, and their rlmm weights indicating outliers can be provided on request for analyses and applications. For instance, Fig.11 suggests that the number of detected outliers increases towards higher frequencies, along with $\tau_{L2L}$. A few of these outliers are also well-sampled localities with a low $SE(\delta L2L_l)$; implying a more source specific study could be worthwhile.

Fig.12 maps the various tectonic localities (indexed $l$) color coded to their $\delta L2L_l(f \approx 0.3,1,3,10Hz)$ values. The colors in panel corresponding to $f \approx 0.3, 1Hz$ (in the top row) are fainter compared to those of $f \approx 3, 10Hz$ (bottom right and left, respectively); indicating the greater diversity in source characteristics and the larger $\tau_{L2L}$ in Fig.6. In the bottom-left panel showing spatial variability of $\delta L2L_l(f \approx 10Hz)$, a strong contrast in source characteristics between central Italy and western Anatolia is apparent. This means that, despite being of similar magnitude, earthquakes originating in central Apennine fault systems generated substantially lower high frequency ground-motions than those originating in north-western segment of Anatolian fault (near Istanbul, Turkey).
1

<figure><img src="image_14.png" /><figcaption>**Fig.12** Tectonic localities and the $\delta L2L_t$ for these locality polygons. Red colored localities produced events capable of higher GM($f$ = 0.3,1,3,10Hz) than the GMM median predictions, blue colored localities produced earthquakes weaker than dataset average, and grey colored localities are close to the dataset average. Localities with fewer ground-motion observations, thereby larger epistemic error on $\delta L2L_t$, are more transparent and appear fainter. Note that the color scale is the same in all the maps in order to emphasize the frequency dependence of spatial variability of event characteristics</figcaption></figure>

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

Similar but lower contrast is apparent when comparing the $\delta L2L_t$ ($f \approx 10Hz$) of central Apennines and Po-plain tectonic localities. The M6.5 Norcia earthquake and associated shocks (in 2016), and several well studied past earthquake sequences (e.g. L'Aquila sequence of 2009) originated in this central Apennines tectonic locality. On the other hand, the substantially stronger (with higher stress-drop) M6.45 Friuli earthquake (in 1976) and a few recent earthquakes are allotted to the Po-plain tectonic locality. Initially, we have considered that the spatial trends (of $\delta L2L_t$) in Fig.12 reflect the differences in observed ground-motions as governed by the predominant focal mechanisms in these localities. However, the within-locality diversity of focal mechanisms dissuaded us from this hypothesis.

It is important to note that the color scale ranges in Fig.12 are not frequency dependent. Despite, an interesting feature in Fig.12 is the inversion of $\delta L2L_t$ in the central Apennines from positive values at $f \approx 0.3Hz$ to negative values at $f \approx 10Hz$ (red to blue color). It is however inconclusive if the events in this region produced low frequency ground-motions stronger than pan-European average or if it is the inhomogeneity of $M_W$ estimates across the pan-European region that is being captured by the $\delta L2L_t$ at $f \le 0.3Hz$.

17

18

19

20

21

We have already seen in Fig.9 the clear correlation between $\delta c_{3,r}$ and the AIx values at station locations. Fig.13 is similar in description to Fig.9, but instead of AIx at each station location within an attenuating region, we use the AIx at event locations within a tectonic locality. Although not presented here, we found no correlation between $\delta B_{e,l}^0$ vs AIx at any frequency (see supplement Fig.S6). It means to say that when $\delta B_{e,l}^0$ variability is large within a tectonic locality, the much larger localities composed
1
2
3
4

of several 0.5º grid cells with similar AIx values are unlikely to resolve event-specific differences. On the other hand, the size of tectonic localities is comparable to the size of regions with distinguishable AIx values (Fig.4). As a result, Fig.13 shows an interestingly strong relationship between δL2Ll and AIx (at event locations).

<figure><img src="image_15.png" /><figcaption>**Fig.13** δL2L<sub>l</sub> of the 55 tectonic localities versus AIx at event locations within each locality, for *f* = 0.3, 1, 3, 10Hz (clockwise from top-left to bottom-left). The blue lines are loess fits between the two parameters. Marker colors indicate the weight assigned to δL2L<sub>l</sub> of each locality in the rlmm regression</figcaption></figure>

5
6
7
8
9
10
11
12

Up to *f* ≈ 1Hz, we observed no resolvable trends between δL2Ll and AIx. Moving towards higher frequencies, as τ<sub>L2L</sub> gains relevance, a significant negative correlation is observed. Essentially, the loess fits for *f* ≥ 2Hz suggest that the tectonic localities coinciding with regions with *AI* ≥ 0.7 are more likely to produce, on average, weaker high frequency ground-motions than those with *AI* < 0.7. However, since AIx is a fuzzy combination of three physical parameters, it is not obvious which one is responsible for the negative correlation with δL2Ll (see supplement Fig.S2 for correlation with seismic moment rate density).

13
14
15
16
17
18
19
20
A classical hypothesis has been that events in stable continental regions produce stronger high frequency ground-motions than those in active crustal regions by virtue of their larger stress-drops, e.g. Bommer et al. (2015). In this case, stable continental regions are those with lower AIx and appear to coincide with tectonic localities with higher δL2Ll values in Fig.13. Alternatively, the large τL2L value at high frequencies could be from regional variability of a high frequency source parameter, e.g. the κ<sup>source</sup>; which is the high frequency decay parameter of Brune's ω<sup>2</sup> source model (Bindi et al. 2019b). Bindi and Kotha (2020) estimated the κ<sup>source</sup> of several events in the ESM dataset, but the large variability of, and uncertainty in, deter us from associating these estimates with the already uncertain
1 $\delta L2L_t$. Although not reported here, we do observe a reasonable positive correlation between $\kappa^{\text{source}}$ and Aix (see supplementary figure).

3 Through an exploratory analysis, we deduced a third hypothesis by comparing the three best-sampled localities in: central Appenines of Italy (“PTTC007” - $\delta L2L_t(f = 10\text{hz}) = -0.6$), north-west segment (“SETC003” - $\delta L2L_t(f = 10\text{hz}) = 0.3$), and the north-east segment of Anatolian fault system (“AMTC001” - $\delta L2L_t(f = 10\text{hz}) = -0.3$). These three regions have similar seismic moment rate density, 1 Hz coda $Q_0$, and shear wave velocity variation at 175km depth, and therefore, similar AIx values (see Fig.4). These tectonic localities produced several $M_W \ge 6$ earthquakes as well. Bindi and Kotha (2020) showed that the stress-drop estimates for the large events in north-west Anatolia are a magnitude higher than those of the central Appenines; while those of north-east Anatolia were unconstrained due to lack of sufficient station coverage. Perrin et al. (2016) showed that the structurally immature north-west Anatolia produced large earthquakes with systematically lower rupture speed than those originating in the structurally mature north-east Anatolian segment. Chounet et al. (2018) followed up with a global database showing an anticorrelation between stress-drop and rupture speed. An earlier study by Radiguet et al. (2009) connected the lower near-field ground-motions (e.g. peak ground acceleration) to structural maturity of fault systems. Based on these studies, we infer that the structurally immature fault segments produced earthquakes capable of stronger high frequency ground-motions. In this study, we concur with the hypothesis that the younger, westward growing north-west Anatolian fault system produced stronger earthquakes than the older, north-east Anatolian and central Appenine fault systems.

21 In summary, $\delta L2L_t$ is introduced to capture partially, in a predictable way, the spatial variability of event dependent ground-motion variability. While the location corrected event-specific $\delta B_{e,l}^0$ retains a large party of event-to-event stress-drop variability as shown in Bindi et al. (2019b); $\delta L2L_t$ captures the average regional trends, which in-turn appear to be controlled by macroscopic fault system characteristics – and not necessarily by the stress-drop (see supplementary Fig.S4). Given its strong correlation with Activity Index, we hypothesize that the trends shown in Fig.13 may assist in *guessing* an appropriate $\delta L2L_t$ for low-moderate seismicity regions (e.g. France and Germany) whose ground-motion data is absent in ESM i.e., more positive $\delta L2L_t$ in regions with lower AIx values (in Fig.4).

29 **6.3 EVENT VARIABILITY**

30 Bindi and Kotha (2020) and a few earlier studies demonstrated a strong correlation between the traditional between-event random-effect $\delta B_e$ and stress-drop. We found a similarly strong correlation between the $\delta L2L_t$ corrected between-event random-effect $\delta B_{e,l}^0$ and stress-drop here as well (supplement Fig.S5) – but no correlation with AI (supplement Fig.S6). Stress-drop is not quite the spatiotemporally predictable parameter, and can be quite variable within any tectonic locality. However, it is commonly acknowledged that ground-motion variability among large magnitude events is smaller
1

2

3

than small events. This has been verified by Youngs et al. (1995) as a magnitude dependence of ground-motion variability – irrespective of sample size. Several subsequent GMM-SAs have since featured the so-called heteroskedastic between-event variance ($\tau_e$) models.

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

**Fig. 14** Magnitude and frequency dependent heteroskedasticity of between-event variability $\tau_0$. Top row panels show $\delta B_{e,l}^0 \sim M_W$ at $f \approx 0.3, 1, 3$ Hz, wherein the marker colors are rlmm weights, black dotted lines delimit $\pm\tau_0$, and red error bars indicate the Median Absolute Deviation (MAD) in $M_W$ bins (details in main text). Bottom row panels show the fitted magnitude dependent $\tau_0$ in red against the black dashed line of homoskedastic $\tau_0$

4

5

6

7

8

9

10

11

12

13

14

15

16

Magnitude-dependent heteroskedasticity of $\tau_e$ makes physical sense because larger ruptures tend to occur more-or-less on the same area of the fault plane, and periodically release a similar amount of accumulated elastic energy or stress. However, if the large events of similar magnitude originate in fault systems (or tectonic localities) with very different stress accumulation and release rates, the variability of stress-drop (and the $\delta B_e$) could be as large as that of smaller events. This means that, heteroskedasticity of $\tau_e$ in a global dataset consisting of events from Italy, Turkey, Japan, China, etc., (e.g. NGA-West2) may not be as significant as that of events within a tectonic locality within any of these geopolitical boundaries. In our case, we remove the locality-to-locality variability of event characteristics through the random-effect $\Delta L2L_l = N(0, \tau_{L2L})$, leaving us a large sample of $\Delta B_{e,l}^0 = N(0, \tau_0)$ that can be examined for a generic heteroskedasticity with $M_W$, independent of events' localities. Based on this understanding, we have modeled a frequency and $M_W$ dependent heteroskedastic $\tau_0$, as illustrated in Fig. 14
In the top panel of Fig.14, we show the $\delta B_{e,l}^0 \sim M_W$ color coded by their rlmm weights. Overlain on the scatter plot are error-bars (in red), which bound $\pm MAD$ (median absolute deviation) of $\delta B_{e,l}^0$ binned according to $M_W \in [3.4, 4), [4,5), [5,5.5), [5.5,6), [6,6.5), [6.5,7.4]$. As in most GMM-SAs, our GMM-FA between-event variability decreases with increasing $M_W$. At $M_W < 5$, the MAD is comparable to the $\tau_0$, while at $M_W \ge 6.5$ the MAD is on average 50% smaller than $\tau_0$ – depending on the frequency. In the lower panels of Fig.14, we show the estimated $M_W$ binwise MAD (red markers), and the piece-wise linear model of heteroskedastic between-event variability, shown in equation (5).

$$ \tau_0(M_W) = \begin{cases} \tau_{0,M_1} = \tau_0 & M_W < M_1 = 5 \\ \tau_{0,M_1} + (M_W - M_1) \left( \frac{\tau_{0,M_2} - \tau_{0,M_1}}{M_2 - M_1} \right) & M_1 \le M_W < M_2 = 6.5 \\ \tau_{0,M_2} & M_W \ge M_2 \end{cases} \quad (5) $$

In equation (5), for each frequency ($f$, not shown in equation), $\tau_{0,M_1}$ is the between-event variability for events with $M_W < M_1 = 5$, which is set equal to $\tau_0$. For large events with $M_W \ge M_2 = 6.5$, the between-event variability is a smaller $\tau_{0,M_2}$. For $M_1 \le M_W < M_2$, the between-event variability decreases linearly from $\tau_{0,M_1}$ to $\tau_{0,M_2}$. Values of these parameters are provided in the supplementary coefficient table.

## 6.4 SITE-RESPONSE VARIABILITY

The next, and by far the largest, random variability is the site-response component $\Delta S2S_s = \mathcal{N}(0, \phi_{S2S})$. Fig.6 shows that $\phi_{S2S}$ is consistently the largest random-variance at all frequencies, when measured as site-to-site ground-motion variability across the 1622 sites in ESM dataset. Site-specific ground-motion predictions are more accurate and precise (smaller aleatory variability) than ergodic or region-specific predictions, but are only possible when site-specific ground-motion data are available. In absence of site-specific observations, site-response proxies are necessary to extrapolate spatially the site-specific terms $\delta S2S_s$ (Kotha et al. 2018; Weatherill et al. 2020b). For such studies, we disseminate a database of $\delta S2S_s(f = 0.15 - 20Hz)$ derived from the ESM dataset.

Fig.15 shows the relation of $\delta S2S_s(f \approx 0.3,1,3,10Hz)$ with measured $V_{s30}$ (left column), and topographic slope (right column). While only 400 sites are provided with measured $V_{s30}$ in the ESM dataset, topographic slope is available at all site locations. In Fig.15, each marker corresponds to a site with an estimated $\delta S2S_s$ (irrespective of number of records), color coded to their rlmm weight. The error-bars (red) illustrate the mean and MAD within each Eurocode 8 site-class. The blue curve represents the proposed linear soil-response model, derived as a quadratic function of $V_{s30}$ or slope.
<figure><img src="image_17.png" /><figcaption>**Fig.15** Plot showing the trend of δS2S<sub>s</sub> (f ≈ 0.3, 1, 3, 10Hz) (from top to bottom rows) with V<sub>s30</sub> (left column) and topographic slope at site location (right column). Note that only 400 sites have measured V<sub>s30</sub> in the ESM dataset, while all 1622 sites have topographic slope values.</figcaption></figure>

Looking at the δS2S_s trends with V_s30, it is evident that sites in EC8 class D (V_s30 ≤ 180m/s) and class C (180 < V_s30 ≤ 360m/s) significantly amplify low frequency ground-motions compared to the average of the dataset. In addition, the within class site-to-site variability (error-bar) is apparently larger than that of EC8 class B (360 < V_s30 < 800m/s) and A (800m/s < V_s30). The site-response of class A sites does not appear to scale with V_s30 at low frequencies. The converse is observed at high frequencies, wherein, class A and B sites exhibit higher site-to-site variability, and scale steeply with
$V_{s30}$. Interestingly, the flattening of $f \approx 3, 10Hz$ site-response function (blue curve) towards lower $V_{s30}$ suggests that class C and D sites exhibit high frequency amplifications lower than that a linear $V_{s30}$ scaling function would predict. Although this is likely from nonlinear behavior of soft soils when subjected to strong input ground-motion, a further record-to-record investigation is necessary. Alternatively the left-over residuals at these sites can be examined for nonlinear soil response, as in Loviknes et al. (2021).

<figure><img src="image_18.png" /><figcaption>Fig.16 Between-site and total variance estimates ($\phi_{S2S}$, $\sigma$) of the GMM-FA for $f = 0.15 - 20Hz$ (in left-panel) and of the Kotha et al. (2020) GMM-SA for $T = 0.01 - 8s$ (in right-panel), along with the reduction in ($\phi_{S2S}$, $\sigma$) to ($\phi_{S2S}^{V_{s30}}$, $\sigma^{V_{s30}}$) when using $V_{s30}$, and to ($\phi_{S2S}^{slope}$, $\sigma^{slope}$) when using Slope as site-response proxies</figcaption></figure>

For now, we only provide the site-response ($SR^{V_{s30}}$, $SR^{slope}$) as quadratic functions of $V_{s30}$ and slope, along with a database of site-response terms, as in equations (6) and (7). Here, the measured $V_{s30}$ is in $m/s$ and slope in $m/m$, the regression coefficients $g_0, g_1, g_2$ are different for $SR^{V_{s30}}$ and $SR^{slope}$, and change with frequency. Robust linear fits using an M estimator (Venables and Ripley 2002), at each of 25 frequencies between $f = 0.15Hz - 20Hz$, PGA ($T = 0s$) and PGV, are derived for $\delta S2S_s \sim V_{s30}$ correlation of 400 sites with measured $V_{s30}$ available, and $\delta S2S_s \sim$ slope of the 1622 sites with slope derived from digital elevation models provided by Shuttle Radar Topography Mission (Jarvis et al. 2008). The residuals from $\delta S2S_s \sim V_{s30}$ (measured) and $\delta S2S_s \sim$ slope regressions are $\delta S2S_s^{V_{s30}}$ and $\delta S2S_s^{slope}$, with robust standard-deviations $\phi_{S2S}^{V_{s30}}$ and $\phi_{S2S}^{slope}$, respectively. All the coefficients and standard-deviations are provided in the supplementary coefficient table.

$$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)$$

Slope is a poorer explanatory parameter of site-response than measured $V_{s30}$, but is relatively easier to obtain at any site location. The left panel of Fig.16 shows the reduction in $\phi_{S2S}$ from using
1

2

3

4

5

6

$V_{s30}$ and slope as site-response proxies in the GMM-FA. Although the reduction at lower frequencies is substantial, an explanation for the large variability of high frequency site-response is still evasive. A few studies, e.g. Hollender et al. (2020), Mucciarelli et al. (2017), and Stewart (2000), discussed the impact of instrument housing and soil-structure interaction on ground-motion recordings; which might explain the relatively large site-to-site response variability at high frequencies, even among sites with similar $V_{s30}$, compared to that at lower frequencies.

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

Also shown in this figure, in the right panel, is the reduction in $\phi_{S2S}$ from using $V_{S30}$ and slope as site-response proxies in the GMM-SA of Kotha et al. (2020). While the variances at $f \le 3Hz$ and $T \ge 0.3s$ are comparable, the ground-motion variability at short periods (e.g. $T < 0.3s$) captured by GMM-SA is a severe underestimate of the actual variability of high frequency ground-motions (e.g. $f > 3Hz$). Differences as such were also observed while comparing the between-locality and between-region random-variances $\tau_{L2L}$ and $\tau_{c3}$, respectively. This observation reinforces the importance of GMM-FAs over GMM-SAs in the various practical applications we highlighted in the Introduction section, and also our motivation to evalute the physical meaning of between-region $\Delta c_{3,r}(f)$ and between-locality $\Delta L2L_l(f)$ random-effects in Fourier domain instead of in response spectral domains – especially at high frequencies, where their variances are largest (see Fig.6). In case of the between-locality and between-region random-effect groups, the correlations $\delta L2L_l \sim AI$ in Fig.13 and $\delta c_{3,r} \sim AI$ in Fig.9 explain to a good extent the large $\tau_{L2L}$ and $\tau_{c3}$ values at high frequencies. But the large $\phi_{S2S}$ at high frequencies could not be resolved using Activity Index. As stated above, the variability of site-response at high frequencies could be from instrument housing type or more physical; such as the influence of deeper soil layers, plasticity of soil column, weathering of bedrock, seasonal changes in shear-wave velocity in shallow layers (Alexis et al. 2021; Roumelioti et al. 2020), etc. However, in-lieu of resolving such issues, and given the practical importance of high frequency site-response, the most efficient solution yet is to collect more site-specific data (Bard et al. 2019).

25

26

27

28

29

30

31

32

33

Apart from the random-effects analyses presented in the previous section, customary checks for magnitude and distance dependencies showed no peculiarities – therefore, we skip presenting them here. Finally, the aleatoric ‘left-over’ residuals can be examined for evidences of some secondary phenomenon which are record-specific and are not captured by the mixed-effects such as: shear-wave radiation pattern in near-source distances (Kotha et al. 2019), SmS reflections at intermediate- and far-source distances (Bindi and Kotha 2020; Bindi et al. 2006), nonlinear soil response (Loviknes et al. 2021), etc. However, each of these phenomena require information that are not currently available in the ESM dataset (e.g. centroid-moment-tensor solutions and crustal shear-wave velocity profiles). For sake of brevity of the manuscript, we anticipate in-depth residual analyses to follow-up studies.
# 7 DISCUSSION

We presented in this study a ground-motion model capable of predicting Fourier amplitudes in the frequency range $f = 0.15 – 20Hz$. This model is developed as complementary to the regionally adaptable ground-motion model for response spectra predictions by Kotha et al. (2020). Hence it’s applicable range and application strategy are the same as those described in Kotha et al. (2020) and Weatherill et al. (2020c). However, unlike Bayless and Abrahamson (2019b), we have not smoothed the fixed-effects coefficients to extend the applicability range of our model. Therefore, we advice against extrapolating the model to higher or lower frequencies. The random-effect and residual variances of our model are large, but are still comparable to the recent models by Stafford (2017) and Bora et al. (2019). Interestingly, all these models’ variances show striking similarity in the rapid increase of estimated ground-motion variability at the high frequencies – a feature that is completely masked in the variance estimates of response spectra ground-motion model.

The three types of partial non-ergodicity are captured by three random-effect groups: between-locality, between-region, and between-site. Random-variances of all three groups increase rapidly at $f > 3Hz$. This intriguing feature was not apparent in the random-variances of Kotha et al. (2020) at periods $T < 0.3s$. Therefore, we found it more appropriate to evaluate the physical meaning of between-locality and between-region random-effects in Fourier domain, which is more closely related to actual physics of seismic ground-motion than the response spectral domain. We have attempted correlating these random-effects to a suite geological and geophysical parameters, but have not reported them here for various reasons – and of course for the sake of brevity. The most convincing correlation, however, was with Activity Index of Chen et al. (2018). This fuzzy logic based transparent and data-driven geospatial index is able to resolve the high between-locality and between-region random-variances at high frequencies. Our interpretation is that regions with high Activity Index (high seismic moment rate density, low 1HZ coda Q, low crustal shear-wave velocity at 175km) attenuate the high frequencies much faster/stronger than regions with low Activity Index (low seismic moment rate density, high 1HZ coda Q, high crustal shear-wave velocity at 175km). Similarly, tectonic localities characterised by low Activity Index (e.g. moderate seismicity regions like France) produced more energetic earthquakes, on average, compared to localities with high Activity Index (e.g. high seismicity regions like central Italy). Activity Index as a composite parameter correlates better to the spatial variability of between-locality and between-region random-effects, than its components individually.

Nevertheless, we find a strong potential in Activity Index as regionalisation parameter in our future ground-motion models. There are several directions we foresee: (1) using Activity Index as *regionalizing* fixed-effect in ground-motion models, (2) using Activity Index to adapt the between-locality and between-region random-effects to sparsely sampled low seismicity areas, and (3) upgrading the Activity Index with newer/updated maps of Moho depth, lithospheric temperature, crustal shear-
1

2

3

4

5

wave velocity maps, etc. The direction (2) is already on course using the recently published ground-motion dataset for France (Traversa et al. 2020) and a weak motion dataset for central and north Europe (Zaccarelli et al. 2019). We strongly encourage the various ground-motion seismology communities to collaborate in testing our inferences and in adapting our model to their datasets; especially those not integrated into ESM.

6

7

8

9

10

11

12

13

Despite our efforts however, the very high site-response variability at high frequencies remains a glaring challenge. Site-specific ground-motion prediction using the site-specific random-effects seems to be the most reasonable approach yet – if we were to avoid the very large ergodic aleatory variability in seismic hazard and risk assessments. It is however quite encouraging to see several research groups attempting to investigate and parametrize site-specific characteristics (e.g. shear-wave velocity profiles, sediment thicknesses, topography, etc) for use in ground-motion models. To aid these efforts and many other practical applications, we consider ground-motion prediction models in Fourier domain are extremely useful; especially in quantifying the high frequency ground-motion variabilities.

14

15

16

17

18

19

20

21

22

23

24

One important aspect we have skipped discussing in this study, for the sake of brevity, is the analyses of aleatory residuals. These residuals are record-specific and also contain useful information on secondary phenomena, such as the shear-wave radiation pattern in near-source distances, arrival of secondary phases at intermediate-far source distances, nonlinear response of soft soils, etc. All these phenomena can be included into a ground-motion model for more realistic (e.g. anisotropic) predictions. Of course, to quantify these effects would require much more information (e.g. centroid moment tensor solutions and incidence angles) in ground-motion datasets. Currently, we are attempting to populate the ESM dataset with as much of relevant information as possible, in order to facilitate further realism of ground-motion models. Based on the experience from and utility of this Fourier amplitude ground-motion model, we plan to pursue such improvements preferably in Fourier domain rather than in response spectra domain.

25
# 8 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. All the outputs of the robust linear mixed-effects regressions are available for dissemination upon request.
1

9 ACKNOWLEDGMENTS

2

3

4

5

6

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 feedback provided by the collaborators in Horizon 2020 “Seismology and Earthquake Engineering Research Infrastructure Alliance for Europe (SERA)” project.
10 REFERENCES

1

2 Abrahamson NA, Silva WJ, Kamai R (2014) Summary of the ASK14 Ground Motion Relation for

3 Active Crustal Regions Earthquake Spectra 30:1025-1055 doi:10.1193/070913eqs198m

4 Akkar S et al. (2014) Reference database for seismic ground-motion in Europe (RESORCE) Bulletin

5 of earthquake engineering 12:311-339 doi:10.1007/s10518-013-9506-8

6 Al Atik L, Abrahamson N, Bommer JJ, Scherbaum F, Cotton F, Kuehn N (2010) The variability of

7 ground-motion prediction models and its components Seismological Research Letters

8 81:794-801

9 Alexis R, Sokos E, Lefils V, Briole P (2021) Seasonal variations in amplitudes and resonance

10 frequencies of the HVSR amplification peaks linked to groundwater Geophysical Journal

11 International

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

13 Atik LA, Youngs RR (2014) Epistemic uncertainty for NGA-West2 models Earthquake Spectra

14 30:1301-1318

15 Bard P-Y, Bora SS, Hollender F, Laurendeau A, Traversa P (2019) Are the Standard V S-Kappa Host-

16 to-Target Adjustments the Only Way to Get Consistent Hard-Rock Ground Motion

17 Prediction? Pure and Applied Geophysics:1-20

18 Basili R et al. (2019) NEAMTHM18 Documentation: the making of the TSUMAPS-NEAM Tsunami

19 Hazard Model 2018

20 Bayless J, Abrahamson NA (2019a) An Empirical Model for the Interfrequency Correlation of

21 Epsilon for Fourier Amplitude SpectraAn Empirical Model for the Interfrequency

22 Correlation of Epsilon for Fourier Amplitude Spectra Bulletin of the Seismological Society

23 of America 109:1058-1070

24 Bayless J, Abrahamson NA (2019b) Summary of the BA18 Ground-Motion Model for Fourier

25 Amplitude Spectra for Crustal Earthquakes in California Bulletin of the Seismological

26 Society of America 109:2088-2105

27 Bindi D, Cotton F, Kotha SR, Bosse C, Stromeyer D, Grünthal G (2017a) Application-driven ground

28 motion prediction equation for seismic hazard assessments in non-cratonic moderate-

29 seismicity areas Journal of Seismology 21:1201-1218 doi:10.1007/s10950-017-9661-5

30 Bindi D, Kotha S-R, Weatherill G, Lanzano G, Luzi L, Cotton F (2019a) The pan-European

31 engineering strong motion (ESM) flatfile: consistency check via residual analysis Bulletin

32 of Earthquake Engineering 17:583-602

33 Bindi D, Kotha S (2020) Spectral decomposition of the Engineering Strong Motion (ESM) flat file:

34 regional attenuation, source scaling and Arias stress drop Bulletin of Earthquake

35 Engineering:1-26

36 Bindi D, Parolai S, Grosser H, Milkereit C, Karakisa S (2006) Crustal attenuation characteristics in

37 northwestern Turkey in the range from 1 to 10 Hz Bulletin of the Seismological Society of

38 America 96:200-214

39 Bindi D, Picozzi M, Spallarossa D, Cotton F, Kotha SR (2019b) Impact of Magnitude Selection on

40 Aleatory Variability Associated with Ground-Motion Prediction Equations: Part II—

41 Analysis of the Between-Event Distribution in Central Italy Bulletin of the Seismological

42 Society of America 109:251-262

43 Bindi D, Spallarossa D, Pacor F (2017b) Between-event and between-station variability observed

44 in the Fourier and response spectra domains: comparison with seismological models

45 Geophysical Journal International 210:1092-1104

46 Bindi D, Zaccarelli R, Kotha SR (2021) Local and Moment Magnitude Analysis in the Ridgecrest

47 Region, California: Impact on Interevent Ground-Motion Variability Bulletin of the

48 Seismological Society of America 111:339-355

49 Bommer JJ et al. (2015) A SSHAC level 3 probabilistic seismic hazard analysis for a new-build

50 nuclear site in South Africa Earthquake Spectra 31:661-698
1 Boore DM, Di Alessandro C, Abrahamson NA (2014a) A generalization of the double-corner-
2 frequency source spectral model and its use in the SCEC BBP validation exercise Bulletin
3 of the Seismological Society of America 104:2387-2398

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

7 Bora SS, Cotton F, Scherbaum F (2019) NGA-West2 Empirical Fourier and Duration Models to
8 Generate Adjustable Response Spectra Earthquake Spectra 35:61-93
9 doi:10.1193/110317eqs228m

10 Bora SS, Cotton F, Scherbaum F, Edwards B, Traversa P (2017) Stochastic source, path and site
11 attenuation parameters and associated variabilities for shallow crustal European
12 earthquakes Bulletin of Earthquake Engineering 15:4531-4561 doi:10.1007/s10518-017-
13 0167-x

14 Bora SS, Scherbaum F, Kuehn N, Stafford P (2014) Fourier spectral-and duration models for the
15 generation of response spectra adjustable to different source-, propagation-, and site
16 conditions Bulletin of Earthquake Engineering 12:467-493

17 Bora SS, Scherbaum F, Kuehn N, Stafford P (2016) On the relationship between Fourier and
18 response spectra: Implications for the adjustment of empirical ground-motion prediction
19 equations (GMPEs) Bulletin of the Seismological Society of America 106:1235-1253

20 Brune JN (1970) Tectonic stress and the spectra of seismic shear waves from earthquakes Journal
21 of geophysical research 75:4997-5009

22 Campbell KW, Bozorgnia Y (2014) NGA-West2 Ground Motion Model for the Average Horizontal
23 Components of PGA, PGV, and 5% Damped Linear Acceleration Response Spectra
24 Earthquake Spectra 30:1087-1115 doi:10.1193/062913eqs175m

25 Chen Y-S, Weatherill G, Pagani M, Cotton F (2018) A transparent and data-driven global tectonic
26 regionalization model for seismic hazard assessment Geophysical Journal International
27 213:1263-1280

28 Chounet A, Vallée M, Causse M, Courboulex F (2018) Global catalog of earthquake rupture
29 velocities shows anticorrelation between stress drop and rupture velocity Tectonophysics
30 733:148-158

31 Cong L, Mitchell B (1998) Lg coda Q and its relation to the geology and tectonics of the Middle
32 East. In: Q of the Earth: Global, Regional, and Laboratory Studies. Springer, pp 563-585

33 Cotton F, Pousse G, Bonilla F, Scherbaum F (2008) On the discrepancy of recent European ground-
34 motion observations and predictions from empirical models: Analysis of KiK-net
35 accelerometric data and point-sources stochastic simulations Bulletin of the Seismological
36 Society of America 98:2244-2261

37 Diaferia G, Cammarano F, Faccenna C (2019) Thermal structure of a vanishing subduction system: an example of seismically-derived crustal temperature along the Italian peninsula
38 Geophysical Journal International 219:239-247

40 Douglas J (2011) Ground-motion prediction equations 1964-2010. Pacific Earthquake
41 Engineering Research Center Berkeley, CA,

42 Douglas J (2018) Calibrating the backbone approach for the development of earthquake ground
43 motion models Best Practice in Physics-based Fault Rupture Models for Seismic Hazard
44 Assessment of Nuclear Installations: Issues and Challenges Towards Full Seismic Risk
45 Analysis

46 Faccioli E, Paolucci R, Vanini M (2015) Evaluation of probabilistic site-specific seismic-hazard
47 methods and associated uncertainties, with applications in the Po Plain, northern Italy
48 Bulletin of the Seismological Society of America 105:2787-2807

49 Fukushima Y (1996) Scaling relations for strong ground motion prediction models with M 2 terms
50 Bulletin of the Seismological Society of America 86:329-336

51 Garnier S (2019) Viridis: Default color maps from "matplotlib". 2018 URL https://github.com/sjmgarnier/viridis R package version 03 4:27
1 Goulet C et al. Effective amplitude spectrum (EAS) as a metric for ground motion modeling using Fourier amplitudes. In: 2018 Seismology of the Americas Meeting, 2018.

2
3 Grad M, Tiira T, Group EW (2009) The Moho depth map of the European Plate Geophysical Journal
4 International 176:279-292

5 Grünthal G, Wahlström R (2012) The European-Mediterranean earthquake catalogue (EMEC) for
6 the last millennium Journal of seismology 16:535-570

7 Hollender F, Roumelioti Z, Maufroy E, Traversa P, Mariscal A (2020) Can We Trust High-Frequency
8 Content in Strong-Motion Database Signals? Impact of Housing, Coupling, and Installation
9 Depth of Seismic Sensors Seismological Research Letters 91:2192-2205

10 Holmgren JM, Atkinson GM (2018) Effect of Uncertainty in Source Parameters on Ground-Motion
11 Variability for Potentially Induced Earthquakes in the Central United States Seismological
12 Research Letters 89:702-711 doi:10.1785/0220170208

13 Jacoby WG (2000) Loess:: a nonparametric, graphical tool for depicting relationships between
14 variables Electoral Studies 19:577-613

15 Jarvis A, Reuter H, Nelson A, Guevara E (2008) Hole-Filled Seamless SRTM Data V4: International
16 Centre for Tropical Agriculture (CIAT): http srtm csi cgiar org, accessed 31

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

18 Kohrangi M, Kotha SR, Bazzurro P (2020) Impact of Partially Non-Ergodic Site-specific
19 Probabilistic Seismic Hazard on Risk Assessment of Single Buildings Earthquake spectra
20 (in-review)

21 Koller M (2016) robustlmm: An R package for robust estimation of linear mixed-effects models
22 Journal of statistical software 75:1-24

23 Kotha SR, Bindi D, Cotton F (2016) Partially non-ergodic region specific GMPE for Europe and
24 Middle-East Bulletin of Earthquake Engineering 14:1245-1263

25 Kotha SR, Bindi D, Cotton F (2017) From ergodic to region- and site-specific probabilistic seismic
26 hazard assessment: Method development and application at European and Middle Eastern
27 sites Earthquake spectra 33:1433-1453 doi:10.1193/081016EQS130M

28 Kotha SR, Cotton F, Bindi D (2018) A new approach to site classification: Mixed-effects Ground
29 Motion Prediction Equation with spectral clustering of site amplification functions Soil
30 Dynamics and Earthquake Engineering doi:10.1016/j.soildyn.2018.01.051

31 Kotha SR, Cotton F, Bindi D (2019) Empirical models of shear-wave radiation pattern derived from
32 large datasets of ground-shaking observations Scientific reports 9

33 Kotha SR, Weatherill G, Bindi D, Cotton F (2020) A Regionally Adaptable Ground-Motion Model
34 for Shallow Crustal Earthquakes in Europe Bulletin of Earthquake Engineering 18:4091-
35 4125 doi: 10.1007/s10518-020-00869-1

36 Kowsari M, Halldorsson B, Hrafnkelsson B, Snæbjörnsson JÞ, Jónsson S (2019) Calibration of
37 ground motion models to Icelandic peak ground acceleration data using Bayesian Markov
38 Chain Monte Carlo simulation Bulletin of Earthquake Engineering 17:2841-2870

39 Kuehn NM, Abrahamson NA (2017) The Effect of Uncertainty in Predictor Variables on the
40 Estimation of Ground-Motion Prediction Equations Bulletin of the Seismological Society
41 of America 108:358-370 doi:10.1785/0120170166

42 Lanzano G et al. (2018) The pan-European engineering strong motion (ESM) flatfile: compilation
43 criteria and data statistics Bulletin of Earthquake Engineering:1-22

44 Loviknes K, kotha SR, Cotton F, Schorlemmer D (2021) Testing ground-motion models non-linear
45 amplification factors Bulletin of Seismological Society of America

46 Lu Y, Stehly L, Paul A, Group AW (2018) High-resolution surface wave tomography of the
47 European crust and uppermost mantle from ambient seismic noise Geophysical Journal
48 International 214:1136-1150

49 Mayor J, Traversa P, Calvet M, Margerin L (2018) Tomography of crustal seismic attenuation in
50 Metropolitan France: implications for seismicity analysis Bulletin of Earthquake
51 Engineering 16:2195-2210
1 Mitchell BJ, Cong L, Ekström G (2008) A continent-wide map of 1-Hz Lg coda Q variation across

2 Eurasia and its relation to lithospheric evolution Journal of Geophysical Research: Solid

3 Earth 113

4 Mucciarelli M, Gallipoli MR, Lizza C, Vignola L (2017) Site condition at accelerometric stations in

5 Southern Italy: observations gained from peculiar cases Bulletin of Earthquake

6 Engineering 15:2367-2376

7 Péquegnat C, Guéguen P, Hatzfeld D, Langlais M (2008) The French accelerometric network (RAP)

8 and national data centre (RAP-NDC) Seismological Research Letters 79:79-89

9 Perrin C, Manighetti I, Ampuero JP, Cappa F, Gaudemer Y (2016) Location of largest earthquake

10 slip and fast rupture controlled by along-strike change in fault structural maturity due to

11 fault growth Journal of Geophysical Research: Solid Earth 121:3666-3685

12 Pilz M, Cotton F (2019) Does the One-Dimensional Assumption Hold for Site Response Analysis?

13 A Study of Seismic Site Responses and Implication for Ground Motion Assessment Using

14 KiK-Net Strong-Motion Data Earthquake Spectra 35:883-905

15 Pilz M, Cotton F, Zaccarelli R, Bindi D (2019) Capturing Regional Variations of Hard-Rock

16 Attenuation in Europe Bulletin of the Seismological Society of America

17 Radiguet M, Cotton F, Manighetti I, Campillo M, Douglas J (2009) Dependency of near-field ground

18 motions on the structural maturity of the ruptured faults Bulletin of the Seismological

19 Society of America 99:2572-2581

20 Ritsema J, Deuss aA, Van Heijst H, Woodhouse J (2011) S40RTS: a degree-40 shear-velocity model

21 for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-

22 mode splitting function measurements Geophysical Journal International 184:1223-1236

23 Rodriguez-Marek A et al. (2013) A model for single-station standard deviation using data from

24 various tectonic regions Bulletin of the seismological society of America 103:3149-3163

25 Roumelioti Z, Hollender F, Guéguen P (2020) Rainfall-Induced Variation of Seismic Waves Velocity

26 in Soil and Implications for Soil Response: What the ARGONET (Cephalonia, Greece)

27 Vertical Array Data Reveal Bulletin of the Seismological Society of America 110:441-451

28 Sahakian V, Baltay A, Hanks T, Buehler J, Vernon F, Kilb D, Abrahamson N (2019) Ground-Motion

29 Residuals, Path Effects, and Crustal Properties: A Pilot Study in Southern California Journal

30 of Geophysical Research: Solid Earth

31 Stafford PJ (2017) Interfrequency Correlations among Fourier Spectral Ordinates and

32 Implications for Stochastic Ground-Motion Simulation Interfrequency Correlations among

33 Fourier Spectral Ordinates and Implications Bulletin of the Seismological Society of

34 America 107:2774-2791

35 Stewart JP (2000) Variations between foundation-level and free-field earthquake ground motions

36 Earthquake Spectra 16:511-532

37 Team RC (2013) R foundation for statistical computing Vienna, Austria 3

38 Traversa P et al. (2020) RESIF RAP and RLBP dataset of earthquake ground motion in mainland

39 France Seismological Research Letters

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

41 Weatherill G, Crowley H, Lemoine A, Roullé A, Tourlière B, Kotha SR, Cotton F (2020a) Modelling

42 Seismic Site Response at Regional Scale for the 2020 European Seismic Risk Model

43 (ESRM20) Bulletin of Earthquake Engineering in-preparation

44 Weatherill G, Kotha SR, Cotton F (2020b) Re-thinking Site Amplification in Regional Seismic Risk

45 Assessment Earthquake spectra doi:10.1177/8755293019899956

46 Weatherill G, Kotha SR, Cotton F (2020c) A Regionally-Adaptable "Scaled-Backbone" Ground

47 Motion Logic Tree for Shallow Seismicity in Europe: Application in the 2020 European

48 Seismic Hazard Model Bulletin of Earthquake Engineering

49 Weatherill G, Kotha SR, Cotton F, Bindi D, Danciu L (2020d) Updated GMPE logic tree and rock/soil

50 parameterisation for ESHM18 vol Deliverable 25.4. Seismology and Earthquake

51 Engineering Research Infrastructure Alliance for Europe (SERA),
1 Weatherill G, Pagani M, Garcia J (2016) Exploring earthquake databases for the creation of

2 magnitude-homogeneous catalogues: tools for application on a regional and global scale

3 Geophysical Journal International 206:1652-1676

4 Wickham H, Chang W, Henry L, Pedersen T, Takahashi K, Wilke C, Woo K (2019a) R Package

5 'ggplot2' v. 3.1. 1 Cran R

6 Wickham H, François R, Henry L, Müller K (2019b) dplyr: A Grammar of Data Manipulation. R

7 package version 0.8. 0.1. ed,

8 Youngs R, Abrahamson N, Makdisi F, Sadigh K (1995) Magnitude-dependent variance of peak

9 ground acceleration Bulletin of the Seismological Society of America 85:1161-1176

10 Zaccarelli R, Bindi D, Strollo A, Quinteros J, Cotton F (2019) Stream2segment: An Open-Source

11 Tool for Downloading, Processing, and Visualizing Massive Event-Based Seismic

12 Waveform Datasets Seismological Research Letters 90:2028-2038

13 Zhu C, Weatherill G, Cotton F, Pilz M, Kwak DY, Kawase H (2020) An Open-Source Site Database of

14 Strong-Motion Stations in Japan: K-NET and KiK-net

15
# Figures

<figure><img src="image_19.png" /><figcaption>Figure 1</figcaption></figure>

Number of records, events, and stations used in regression at each frequency
<figure><img src="image_20.png" /><figcaption>Figure 2</figcaption></figure>

Distribution of records in each of the 45 attenuation regionalisation polygons Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.

<figure><img src="image_21.png" /></figure>
<caption>**Figure 3**</caption>

Distribution of records from shallow crustal events in each of the 55 event localisation polygons (tectonic zones) Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.

<figure><img src="image_22.png" /><figcaption>**Figure 4**</figcaption></figure>

Activity Index map of pan-European region by Chen et al. (2018) Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.
<figure><img src="image_23.png" /><figcaption>Figure 5</figcaption></figure>

Predicted Fourier spectra of the GMM-FA for various scenarios, differentiated by color for Shallow (M<10<100), Intermediate (100≤M<200), and Deep (200≤M) events in the left panel. In the right panel are the Fourier spectra differentiated by regional anelastic attenuation: Slower and Faster than pan-European Average
<figure><img src="image_24.png" /><figcaption>Figure 6</figcaption></figure>

Please see the Manuscript PDF file for the complete figure caption.
<figure><img src="image_25.png" /><figcaption>Figure 7</figcaption></figure>

3,¤ for ¤ = 0.15–20¤. Each line corresponds to one of the 45 attenuation regions, with colors indicating their weight in robust regression. Overlain red curves mark the ±¤¤3 values. Regions with 3,¤(¤) beyond ±1.345¤¤3(¤) are given a lower than a unit weight
<figure><img src="image_26.png" /><figcaption>Figure 8</figcaption></figure>

□□3,□ variation across the 45 regions in ESM dataset for □=0.3,1,3,10□□ (clockwise from top-left to bottom-left panels). Blue polygons locate regions with anelastic attenuation faster than the pan-European average, red polygons locate regions with slower attenuation, and grey polygons are regions close to the average. Regions with fewer ground-motion observations, thereby larger epistemic error on □3,□, are more transparent and appear fainter. Note that the color scale is the same in all the maps in order to emphasize the frequency dependence of spatial variability of anelastic attenuation Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.
<figure><img src="image_27.png" /><figcaption>Figure 9</figcaption></figure>

|||3,|| of the 45 regions versus Activity Index at site locations within each region, for ||≈0.3,1,3,10||||
(clockwise from top-left to bottom-left). The blue lines are loess fits between the two parameters. Marker
colors indicate the weight assigned to |||3,|| of each region in the rlmm regression
<figure><img src="image_28.png" /><figcaption>Figure 10</figcaption></figure>

条形图中，条的粗细代表不同地区的活动强度，颜色代表不同的活动类型。其中， Activity Index 表示活动强度，从浅蓝色代表低强度活动到红色代表高强度活动。不同颜色的虚线代表不同的活动范围。
<figure><img src="image_29.png" /><figcaption>**Figure 11**</figcaption></figure>

22 for $\mathbf{X}$ = 0.15-20$\mathbf{X}$. Each line corresponds to one of the 55 tectonic localities, with colors indicating their weight in rlmm regression. Overlain red curves mark the ±2$\mathbf{X}$ values. Localities with 22($\mathbf{X}$) beyond ±1.3452$\mathbf{X}$($\mathbf{X}$) are given a weight lower than one
<figure><img src="image_30.png" /><figcaption>Figure 12</figcaption></figure>

Tectonic localities and the ⌜2 Emergence Emergence⌝ for these locality polygons. Red colored localities produced events capable of higher ⌜2 Emergence⌝( Emergence = 0.3,1,3,10 Emergence ) than the GMM median predictions, blue colored localities produced earthquakes weaker than dataset average, and grey colored localities are close to the dataset average. Localities with fewer ground-motion observations, thereby larger epistemic error on ⌜2 Emergence⌝, are more transparent and appear fainter. Note that the color scale is the same in all the maps in order to emphasize the frequency dependence of spatial variability of event characteristics. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.
<figure><img src="image_31.png" /><figcaption>Figure 13</figcaption></figure>

条形图中的x轴和y轴分别代表什么？

在本例中，我们比较了55个取样地的岩浆喷发事件与A1x（事件时间）之间的关系。为了比较这些事件的岩浆喷发时间，我们使用了不同的时间步长（即时间步进）。我们使用了0.3, 1, 3, 10等不同的时间步长（clockwise from top-left to bottom-left）。这些时间步长的代表颜色为蓝色，从浅蓝色（时间步长最短）到深蓝色（时间步长最长）。每个事件的岩浆喷发时间都被映射到了时间轴上的相应位置。不同事件的岩浆喷发时间被用不同颜色的柱形表示出来， colors indicate the weight assigned to each locality in the r1mm regression。
<figure><img src="image_32.png" /><figcaption>Figure 14</figcaption></figure>

Magnitude and frequency dependent heteroskedasticity of between-event variability ≲0. Top row panels show ≲0.01, ≲0.1 at ≲0.3, 1, 3 ≲0.5, wherein the marker colors are rlmms weights, black dotted lines delimit ± ≲0, and red error bars indicate the Median Absolute Deviation (MAD) in ≲0.5 bins (details in main text). Bottom row panels show the fitted magnitude dependent ≲0 in red against the black dashed line of homoskedastic ≲0
<figure><img src="image_33.png" /><figcaption>Figure 15</figcaption></figure>

Plot showing the trend of ♂♂2♀♀(♀ ≈0.3,1,3,10♀♀) (from top to bottom rows) with ♂♂30 (left column) and topographic slope at site location (right column). Note that only 400 sites have measured ♂♂30 in the ESM dataset, while all 1622 sites have topographic slope values.
<figure><img src="image_34.png" /><figcaption>Figure 16</figcaption></figure>

Please see the Manuscript PDF file for the complete figure caption.

**Supplementary Files**

This is a list of supplementary files associated with this preprint. Click to download.

* Supplementaryfigures.pdf