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

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


Introduction
Ground-Motion Models (GMMs) characterise the random distribution of ground-motions at site for a combination of earthquake source, wave travel-path, and the affected site's geological properties. Typically, GMMs are regressed over a compendium of strong groundmotion recordings collected from several earthquakes recorded at multiple sites scattered across a variety of geographical regions. The necessity of compiling such large datasets is to expand the range of magnitude and distance, and diversity of site-types, in order to derive a GMM capable of predicting realistic ground-motions for prospective seismic risk scenarios, e.g. large magnitudes at short distances from a reference rock site. NGA-West2 (Ancheta et al. 2014) is one such dataset compiled of ground-motion observations recorded around the globe-primarily from Western US, and in smaller fractions from Alaska, China, Italy, Japan, Taiwan, Turkey, etc. Several GMMs have been derived from this dataset for application in probabilistic seismic hazard (PSHA) and risk assessments. Given the clear tectonic and geological diversity of the data, possible regional differences in observed ground-motions needed to be quantified (Douglas 2004). Four of the NGA-West2 GMMs (Abrahamson et al. 2014;Boore et al. 2014;Campbell and Bozorgnia 2014;Chiou and Youngs 2014) accounted for regional differences in ground-motions through region-specific regression coefficients. Through region-specific adjustments, these GMMs were able to capture and predict, for example, the faster attenuation of short-period ground-motions with distance in Japan compared to Western US, and the relatively slower attenuation in China. Earlier GMMs were not capable of such predictions simply owing to the lack of sufficient data from individual regions to quantify the differences.
In case of Europe and the Middle-East, RESORCE (Akkar et al. 2014b) is one such dataset compiled of data from Italy, Turkey, Greece, and other active seismic regions in pan-Europe. Using mixed-effects regression algorithms (Abrahamson and Youngs 1992), a few GMMs (Akkar et al. 2014a;Bindi et al. 2014) were derived and used in regional (Grünthal et al. 2018) and continental PSHA (Giardini et al. 2018;Woessner et al. 2015). However, these GMMs were not regionalised despite known geological differences between Italy and Turkey. In fact, even the regionalised GMMs (Bora et al. 2017;Kale et al. 2015;Kotha et al. 2016;Kuehn and Scherbaum 2016) limited the distinction to geopolitical boundaries, yet the geological diversity within these regions is far more complex. In essence, quantification of regional differences requires, first, a regionalisation scheme, and then, sufficient data from each region. With the arrival of the new European Strong-Motion (ESM) dataset (Bindi et al. 2018b;Lanzano et al. 2019a) and regionalisation models, and ongoing efforts to update the pan-European PSHA, a revision of the regionalised pan-European GMMs is proposed.
In this study, we present an upgrade of the RESORCE dataset based region-and sitespecific GMM (Kotha et al. 2016 for shallow crustal earthquakes using the new ESM dataset. We often see GMMs evolve with progressively larger datasets, and supersede their older versions in terms of applicability ). However, with increasing data and complex parametrization, a reduction in the apparent aleatory variability (σ) of GMMs could not be achieved (Douglas and Edwards 2016;Strasser et al. 2009). Of course, increasing amount of data comes with increasing the spatiotemporal diversity of ground-motion observations, and thus an increasing σ. One approach would be to introduce new predictor variables into the GMM, but then, uncertainty or unavailability of predictor values then becomes an issue during application Kuehn and Abrahamson 2017). With this in mind, this revision of the Kotha et al. (2016) GMM will attempt to regionalise and refine the aleatory variability, while maintaining its original parametrization. In addition, assuming a possible contamination of data and deviation from the assumption of log-normality, instead of the usual ordinary least-square estimates of GMM median and variances we compute their robust counterparts; while flagging outlier events, stations, and records in the dataset. Table 1 summarizes the additional features we introduced in this new GMM (abbreviated K20) with respect to the recent pan-European GMMs: AK14 (Akkar et al. 2014a), B14 (Bindi et al. 2014), K16 (Kotha et al. 2016), L19 (Lanzano et al. 2019b). Figure 1 compares the data distribution between the RESORCE and ESM datasets. The increase in amount of data between 2014 and 2018 for GMM development is dramatic. While the K16 GMM was regressed over 1251 records, the proposed revision (K20) derives from 18,222 records. One striking feature of the ESM dataset is the number of stations with ≥ 3 ground-motion recordings. RESORCE dataset had about 150 stations with ≥ 3 records, while ESM has 1077. This increase is highly sought in empirical site-specific GMM, PSHA, and seismic risk applications (Faccioli et al. 2015;Kohrangi et al. 2020;Kotha et al. 2017;Rodriguez-Marek et al. 2013).

Ground-motion data and selection criteria
Another feature is the remarkable increase in the number of small earthquakes with M W ≤ 4.5 . This could imply an increase in σ of the revised GMM over K16, similar to the reported increase for NGA-West2 GMMs with respect to their superseded NGA-West counterparts. Nevertheless, ground-motions from frequent small-moderate sized earthquakes drive the hazard in low-moderate seismicity regions of Europe, such as France (Drouet et al. 2020) and Germany (Grünthal et al. 2018). It is necessary that the GMM is well behaved in these small-moderate magnitude ranges. Moreover, if the site-specific terms S2S s were to be estimated with low uncertainty, data from several small-moderate sized earthquakes is indispensable.
The distance (Joyner-Boore metric, R JB ) range and density of data is also superior to RESORCE. The region-specific anelastic attenuation terms of K16, NGA-West2 and other GMMs (Sedaghati and Pezeshk 2017) were estimated from records at R JB ≥ 80 km. This revised GMM aims to refine the regionalisation of K16, and therefore, such increase in farsource data is quite useful. In addition, advanced studies on spatial and temporal variability of attenuation are now possible (Bindi et al. 2018a;Dawood and Rodriguez-Marek 2013;Kotha et al. 2019;Landwehr et al. 2016;Piña-Valdés et al. 2018;Sahakian et al. 2019).
The data visualized in Fig. 1 are from shallow crustal earthquakes in the ESM dataset. The full dataset contains ground-motions from other tectonic regimes as well, such as subduction interface, subduction in-slab, Vrancea, etc. To filter out these and other records not suitable for a shallow crustal GMM development, we adopt the following selection criteria: 1. To keep only the shallow crustal earthquakes, we select only those events classified as non-subduction events by Weatherill et al. (2020d). The selection removes inslab, interface, outer-rise, and upper-mantle events from the regression dataset.  (Grünthal and Wahlström 2012) are preferred over the ESM default values 3. We keep all sites in the dataset irrespective of whether their V S30 , measured from geotechnical investigations, is provided or not in ESM. This is to estimate the site-specific terms ( S2S s ) at as many sites as possible, and then explore various site-response proxies to characterize them (Kotha et al. 2018;Weatherill et al. 2020b) 4. Choice of distance metric is R JB where available, otherwise the epicentral distance R epibut only for events with M W ≤ 6.2 . The distance range is not truncated and extends up to R JB = 545 km 5. Only those records with high-pass filter frequency f hp ≤ 0.8∕T , where T is the period, for both horizontal components are kept in the regression of spectral accelerations SA(T). This is to ensure that the filter does not significantly affect the response spectral values. As a result the dataset varies with the T (Abrahamson and Silva 1997). Further detail on waveform processing can be found in Lanzano et al. (2019a).
Following the above criteria, the number of records available for GMM regression is 18,222, from 927 events ( 3.1 ≤ M w ≤ 7.4 ) recorded at 1829 stations ( 0 ≤ R JB < 545 km ). These numbers decrease to 9698, 491, and 1341 respectively due to the high-pass filter criteria (#5) when T = 8 s.

Regionalisation datasets
In order to attempt GMM regionalisation, we first need to define regions, classify the data into regions, and quantify the regional variabilities through regression. The mixed-effects regression estimates first the group random-variances, and then the random-effects for individual levels in the group (Bates et al. 2014). For instance, in estimating the well-known between-event term ( B e ), the preliminary step is to quantify observed event-to-event ground-motion variability as the between-event random-variance (τ 2 ). Following this, the between-event random-effects for each event are estimated from the total-residuals. The procedure is similar for the between-site terms ( S2S s ), and other random-effects in a mixed-effects regression. Essentially, within the event group, individual events are the levels. It is important to note that random-effect values for the levels (in a group) can be used in three ways: (1) in level-specific predictions accounting the epistemic uncertainty of level-specific adjustment, (2) investigated for any physical phenomena and proxy parameters or (3) ignored, while treating the group random-variance as an aleatory variability or an epistemic uncertainty in a GMM logic-tree (Douglas 2018a, b;Weatherill et al. 2020c).
We emphasize that an exploratory analysis on random-effects (of levels in a group) is to evaluate if they concur with known geophysical properties. If yes, the random-effect group remains in the regression; otherwise the grouping needs revision. In this study, we aim to capture variabilities among events, event sources (e.g. fault systems), attenuating regions, and recording sites. Anticipating significant random-variances, we first describe the various regionalisation datasets used to formulate the random-effect groups.

Attenuation regionalisation
Recent GMMs have demonstrated a strong between-region variability of anelastic attenuation, which is partially attributed to spatial variability in crustal characteristics, e.g. the 1 Hz Lg-coda Q values (Cong and Mitchell 1998), shear-wave velocity (Lu et al. 2018), etc. The first generation of regionalised GMMs, however, relied on national administrative boundaries for regionalisation of anelastic attenuation (e.g. Italy, Turkey, Japan, California, etc.)-which could yield incongruous estimates of random-variances. For example, spatial variability of attenuation imaged by 1 Hz coda Q maps of regions within Italy, Turkey (Cong and Mitchell 1998), France (Mayor et al. 2018), and Europe in general (Pilz et al. 2019) is rather high. Such spatial variability cannot be contained within political boundaries. Therefore, in this study, we adopt a more geological-geophysical regionalisation developed by Basili et al. (2019) under the purview of the TSUMAPS-NEAM project. Figure 2 shows the regionalisation and the number of records within each region polygon, as decided by the recording site location. Assuming that the anelastic attenuation is Fig. 2 Number of records in each of the 46 attenuation regionalisation polygons a phenomenon dominant at far-source distances and in near-surface crustal layers, we let the recording site location decide to which region a particular record is to be assigned. Alternatively, regionalising records based on event location would cause some inconsistencies concerning event depth. Seismic waves from shallow and deep events sample different depths of the crust and upper mantle. TSUMAPS-NEAM regions used in this study are surficial only, and cannot accommodate depth dependence of attenuation. Therefore, records are allotted to the regions based on the recording surface site location.
There are 46 regions in this map (Fig. 2), which is substantially larger than K16 with Italy, Turkey, and rest of Europe as the three attenuating regions. This regionalisation scheme splits Turkey into at least eight regions, East and West Anatolia being the largest. Italy is highly fragmented as well. Figure 2 indicates that there are regions with a few thousand records (e.g. Central and Northern Apennines), and a few with less than a hundred records (e.g. only 13 records from Rhine Graben). Thus, random-variance of the attenuation regions group will be a quantification of the observed region-to-region variability in anelastic attenuation across these regions (each region is a level). A finer or coarser regionalisation may yield a different estimate of between-region variance. Customarily, only if the random-variance for this group is significant it would mean the regionalisation is appropriate, and that there is a significant regional variability in anelastic attenuation. Subsequently, we will inspect whether the region-specific random-effect values resemble any physical phenomenon.

Event localisation
Traditionally, earthquake-to-earthquake variability is captured by the between-event random-variance (τ 2 ), and the event-specific random-effect values ( B e ) are estimated for each event in this group. Based on earlier works on fault-maturity (Bohnhoff et al. 2016;Manighetti et al. 2007;Radiguet et al. 2009), we hypothesized that events associated to particular earthquake source (or a fault system) show systematic differences in their ground-motions, and that event-to-event ground-motion variability is similar across various fault-systems.
In this study, we introduced an additional random-effect to quantify earthquake locality-to-locality variability, similar to the location-to-location variability ( L2L l ) defined in Al Atik et al. (2010). We assign the ESM shallow crustal events to seismotectonic zones defined in the European Seismic Hazard Model 2020 (ESHM20) towards development of seismic source models. The seismotectonic source zonation (referred to as 'TECTO') is designed to be large-scale and does not attempt to resolve smaller scale seismogenic features. As such it is intended to reflect the local tectonics influencing the seismic source generation but not to resolve highly localised features, which makes it a good candidate for our purpose-and more appropriate for source localisation than the TSUMAPS-NEAM model used for attenuation regionalisation.
The random-effects group will be referred to as between-locality ( L2L l ) from hereon, with each level (locality l) being a tectonic locality with allotted shallow crustal events. Figure 3 shows the distribution of recordings across the tectonic localities with at least one shallow crustal event associated to them. Once again, this group is introduced to quantify the earthquake locality-to-locality variability of ground-motion in the dataset through the between-locality variance ( 2 L2L ), which if close to zero indicates no regional variability of earthquake characteristics.

Functional form
A mixed-effects GMM is composed of fixed-effects and random-effects. Fixed-effects part of the GMM predicts generic ground-motions as a continuous function of predictor variables, which in this case are event magnitude M W and distance metric R JB (in km). Random-effects can serve as adjustments to the generic fixed-effects regression coefficients in order to produce ground-motion predictions specific to an event, tectonic locality, region, and site. Equation (1) illustrates the decomposition of observed groundmotion values into mixed-effects. We note that Eq. (1) is not an exact mixed-effects formula, but is used here only as an aid to explain the approach.
In the LHS of Eq. (1), ln GM e,l,s,r is the natural-log of ground-motion ( GM e,l,s,r , which are rotD50 measures of SA(T) here) produced by an event (e) originating in a tectonic locality (l in Fig. 3), and recorded by a site (s) located in region (r in Fig. 2). These observed ground-motions are decomposed into mixed-effects. Since events are exclusive to tectonic localities, and sites to regions, the formulation contains nested mixed-effect grouping. Alongside, since events in a locality can be recorded by sites in multiple surrounding regions, the formulation contains crossed mixed-effect grouping as well. Therefore, this GMM is developed using a crossed and nested mixed-effects approach (Stafford 2014).
In the RHS of Eq. (1), ln( ) is the generic fixed-effects formula without any specificities related to event, locality, site, and region. The B 0 e,l , L2L l , S2S s , and c 3,r are random-effect adjustments to the coefficients of the fixed-effect ln( ) ; and are aimed to be specific to the event (e), locality (l), site (s), and region (r), respectively. ϕ is the standard-deviation of the left-over residuals = N 0, 2 , which contain phenomenon that could not be contained in the crossed and nested mixed-effects. Subsequent sections will explain these terms in more detail.

Fixed-effects
Equations (2) through (5) present the fixed-effects formula of the GMM. In Eq.
(2), ln( ) is the generic median prediction as a function of M W and R JB ; where, e 1 is the generic offset, f R,g M W , R JB represents the geometric spreading, f R,a R JB represents the apparent anelastic attenuation, and f M M W represents the magnitude scaling of ground-motion prediction.
The ln( ) is the predicted natural-log of RotD50 (Boore 2010) measures of PGA (in gal) and PGV (in cm s −1 ), and the 5% damped elastic response spectral ordinates in acceleration (SA, in gal) at 34 periods ranging from 0.01 to 8 s. Therefore, the fixed-effects coefficients (2-5) change with the period, PGA and PGV, and are generic to all events, tectonic localities, attenuating regions, and sites.
The fixed-effects component of this GMM remains similar to that of K16, but with a few minor changes based on non-parametric analyses. In Fig. 4, the top row of the plot shows non-parametric scaling of SA(T = 0.1s) with R JB ; wherein, for clarity, the data is split into magnitude bins and plotted in separate panels. The bottom row of the plot shows non-parametric scaling of SA(T = 0.1 s) with M W ; where, the data is split into distance bins and plotted in separate panels. The smooth (coloured) curves are the loess fits (Jacoby 2000) illustrating non-parametric scaling of SA(T = 0.1 s) with magnitude and distance. Given the variety of hypocentral depths in the dataset, and that we are using the depth-insensitive R JB , we categorized the events into three depth bins ( D < 10 km, 10 km ≤ D < 20 km, 20 km ≤ D ), and plotted unique non-parametric attenuation curves for each bin. We note that, although the non-parametric plots presented here are for SA(T = 0.1 s) and the shape of loess fits varies gradually with periods, the differences were not large enough to make changes to the GMM functional form. Based on such plots we introduced some changes into the GMM with respect to K16: 1. In the SA(T = 0.1 s) versus R JB panels, it is evident that deeper events ( 20 km ≤ D ) have a longer near-source saturation plateau-a feature controlled by the so-called h-parameter in GMMs (h D in Eqs. 3 and 4). Shallower events have a shorter near-source saturation plateau, and a steeper decay with distance. However, these differences are only prominent at R JB ≤ 30 km. ( The marginally larger SA(T = 0.1 s) values of 20 km ≤ D event at R JB > 30 km are also reasonable (Derras et al. 2012). Deeper events are closer to Moho; producing weaker surface ground-motions in the epicentral zone due to longer travel paths of direct arrivals, and relatively stronger ground-motions at far-source distances from the more efficient propagation in deeper crustal layers. Based on these observations, we allowed the h-parameter to vary with the depth bins, as indicated by the subscript D in h 2 D of Eqs. (3) and (4). In K16 GMM, this parameter was depth-independent, and had no subscript. Moreover, instead of making this parameter a regression coefficient, we keep the regression linear by assigning a priori values (based on preliminary non-parametric analyses and nonlinear regression trials), which are independent of magnitude and period: h D = 12 km for deep 20 km ≤ D events, h D = 8 km for events of intermediate depth 10 km ≤ D < 20 km , and for shallow events with D < 10 km we assign h D = 4 km 2. In addition to the above, since ground-motions appear reasonably depth independent at R JB ≥ 30 km , we set (3) is the reference magnitude, and remains the same as in K16 i.e., M ref = 4.5 3. In the SA(T = 0.1s) versus M W panels, although the depth dependence of near-source attenuation is not evident, we do observe saturation towards large magnitudes. In the panel showing non-parametric ground-motion scaling with M W , the evidence suggests that saturation begins at M W ≥ 6.2. However, this is a feature most noticeable at R JB ≤ 30 km and at short periods; here SA(T = 0.1 s) . Towards longer periods and at longer distances, the saturation is less pronounced. Therefore, unlike in K16 where the hinge-magnitude was set as M h = 6.75 , in this GMM we set M h = 6.2 . M h is period independent; saturation or otherwise beyond M h is captured by b 3 in Eq. (5).

Random-effects
In the mixed-effects formulation of Eq. (1), ln( ) is the fixed-effects component. Equations (2) through (5) describe its constituents, wherein e 1 , b 1 , b 2 , b 3 , c 1 , c 2 , c 3 are the fixedeffects regression coefficients. These generic fixed-effect coefficients can be adjusted using random-effect estimates to predict event, tectonic locality, attenuation region, and site specific ground-motions. Regarding the random-effects components: 1. c3 quantifies the between-region variability of anelastic attenuations across the attenuation region group described earlier (in Fig. 2). This means that, along with a generic c 3 and random-variance 2 c3 , region-specific adjustments c 3,r are estimated as randomeffects. These random-effects follow a Gaussian distribution Δc 3 = N 0, 2 c3 . The generic fixed-effect coefficient c 3 in Eq. (4) can be adjusted with a region's c 3,r to achieve region-specific coefficient for that region r, as in c 3,r = c 3 + c 3,r . Unlike the generic c 3 , c 3,r is the apparent anelastic attenuation term that varies with region (subscript r identifies regions in Fig. 2). In K16, r identified the regions Italy, Turkey, and the rest of pan-Europe as Other. This grouping ensured that each of the regions had sufficient data to estimate statistically reliable c 3,r . Accordingly, c3 quantified the regional variability of anelastic attenuation when RESORCE dataset is grouped into Italy, Turkey, and Other. ESM contains much more data from Italy and Turkey, and also several other nations. Since the number of regions is 15 times (46 regions) that of K16 (three regions), the quantified regional variability in anelastic attenuation, in terms of c3 , is also larger than that of K16. Regions with sufficient data also benefit from a well-constrained region-specific adjustment c 3,r 2. Between-locality variability of observed ground-motions are captured by the randomeffect ΔL2L = N 0, 2 L2L , where the mixed-effects regression quantifies the variability as L2L , at each T, and earthquake locality-specific terms as L2L l (subscript l identifies regions in Fig. 3). This random-effect can be used to adjust the fixed-effect coefficient e 1 in Eq. (2), as in e 1,l = e 1 + L2L l , to achieve ground-motion predictions specific to the tectonic locality of the events' origin 3. Event-to-event variability in this GMM is the traditional between-event random-effects ΔB e = N 0, 2 filtered for between-locality variability ΔL2L = N 0, 2 L2L , and is now captured by the ΔB 0 e = N 0, 2 0 ; where, for an event e located in tectonic locality l, the event-specific term can be seen as B 0 e,l ≈ B e − L2L l . 0 is the generic event-to-event variability corrected for locality-to-locality variability, and does not vary with tectonic locality l 4. Site-to-site response variability is captured by the site-specific random-effects ΔS2S = N 0, 2 S2S . The potential of S2S s in site-specific GMMs is well-known, and are useful in studying regional differences in site-response scaling with V S30 (timeaveraged shear-wave velocity in 30 m top-soil) as in K16 or other site-response proxies (Kotha et al. 2018;Weatherill et al. 2020b) 5. The left-over residuals = N 0, 2 contain the unexplained natural variability of ground-motion observations, and thus represent the apparent aleatory variability of the model. These residuals can be investigated for less dominant phenomenon, such as the anisotropic shear-wave radiation pattern  In all, this GMM has four random-effect groups, i.e. one degree-of-freedom more than K16, to explain more than 15 times the data. Those common with K16 are refined with a more physical regionalisation scheme, and greater geographical coverage of shallow crustal events and recording sites. With this configuration of mixed-effects GMM, we run a robust linear mixed-effects regression (Koller 2016) independently for the 36 RotD50 IMs of 5% damped SA for T = 0.01 − 8 s , PGA and PGV. Along with the regression coefficients, we estimate and provide also the fixed-effects variance-covariance matrices needed to estimate the GMM epistemic uncertainty (Atik and Youngs 2014;Bindi et al. 2017) and to update the GMM coefficients in a Bayesian framework (Kowsari et al. 2019(Kowsari et al. , 2020Kühn and Scherbaum 2015).

Regression method adapted to large datasets
Unlike K16, where the GMM is regressed using an ordinary least-squares mixed-effects regression algorithm (Bates et al. 2014), in this study we employ a robust mixed-effects regression algorithm (Koller 2016). As any other real-life data, ground-motion data may contain outliers and other contaminations. Even minor contamination may drive the classical ordinary least-square mixed-effect estimates away from those without contamination. Robust linear mixed-effects (rlmm) regression is then quite useful in limiting the influence of outlier events, sites, and records on the GMM median and variances. It is important to note that, these outlier data points are not entirely removed from the GMM regression but are simply down-weighted.
A feature of rlmm relevant to this study is that the random-effects (attenuation regions, events, event localities, and sites) and residuals (records) with values beyond ± 1.345 standard-deviations of their respective normal random-distributions are assigned progressively lower weights (< 1); whereas, in ordinary least-squares all data is assigned unit weight. Any event, site, and record with non-unit weight are considered a possible outlier, and needs to be examined for its peculiarity. Although it is possible to tune the regression parameters, e.g. those of the Huber loss function (Huber 1992) in robust regressions, we chose to remain with the suggested default values (Koller 2016) optimized for efficiency and robustness in detecting outliers.

Results and discussion
The regression results comprise of fixed-effect coefficients and covariance matrices, random-effect values including weights and standard-errors, residuals, and variances. It is customary to check the behaviour of GMM fixed-effects component, and its epistemic uncertainty in various magnitude and distance ranges. Random-effects and residuals are checked for any noticeable biases or trends with predictor variables. We discuss them separately. Figure 5 shows the GMM's median SA(T = 0.1 s) and SA(T = 1 s) (in gal) predictions over magnitude and distance ranges. Along with the median prediction (lines), its epistemic uncertainty in terms of asymptotic standard-deviation ( ± ) is shown as well (coloured ribbons).

Distance scaling
First, we discuss the predicted SA(T = 0.1 s) and SA(T = 1 s) scaling with R JB in the left panel of Fig. 5. Looking at the curves for M4, we notice the impact of depth-dependent h D in rendering a shorter near-source saturation (plateau ~ 0 to 3 km) for shallower events, compared to the intermediate depth events (plateau ~ 0 to 5 km) and deeper events (plateau ~ 0 to 10 km). The three curves merge at about 30 km, which is our R ref . Beyond R ref = 30 km , the depth-dependence of distance scaling is non-existent. The predictions show good agreement with the non-parametric trends in Fig. 4.
Beyond the R ref , we see the region-dependent anelastic adjustments coming into play, which show the impact of adjusting c 3 with c 3,r = 0, − c3 , + c3 . A region with faster than average attenuation will have a c 3,r more negative (with c 3,r = − c3 ) than generic average, and vice versa for slower attenuation (with c 3,r = + c3 ). Which of the 46 regions in pan-Europe apparently attenuate faster or slower than pan-European average (with c 3,r = 0 ) will be discussed in the following sections, and elaborated in a follow-up study. Effect of In Fig. 5, we also show the epistemic uncertainty on median predictions. The orange ribbon is almost too thin to be noticeable for M4 and M5.5 predictions. Only for M7 and larger events, the ribbon is visibly wide because of the limited data from large magnitude events at near-source distances in the ESM dataset, and the uncertainty on the coefficient b 3 , which describes the magnitude scaling beyond M h = 6.2.

Magnitude scaling
In the right panels of Fig. 5 we show the SA(T = 0.1 s) and SA(T = 1 s) scaling with M W . Here as well we show 27 curves, but this time for R JB = 10, 50, 150 km instead of the three M W values. The features in scaling with distance discussed in reference to the left panels also prevail here; c 3,r is effective at R JB = 150 km , while h D is effective at R JB = 5 km , and neither are at R JB = 50 km . More important in this context, is the difference in scaling with magnitude at R JB = 10 km compared to R JB = 50, 150 km . Evidently, the scaling of SA(T = 0.1 s) is more gradual (less steep) at near-source distance for M W < M h and oversaturates for M W ≥ M h . This is a known physical behaviour wherein, near-source groundmotions, especially the short period SAs, are less sensitive to M W (Campbell 1981; Schmedes and Archuleta 2008). A few previous GMMs observed the same with various datasets. Boore et al. (2014) and K16 allowed oversaturation at large magnitudes ( M W ≥ 6.75) , but whether this is realistic or not needs to be verified with specially compiled near-source ground-motion datasets, e.g. the one recently made available by Pacor et al. (2018). The oversaturation effect, however, disappears at longer periods, as shown here for SA(T = 1 s). Figure 6 shows the predicted response spectra for the same scenarios shown in Fig. 5. The effect of event depth and regionalized anelastic attenuation are the same as in Fig. 5. An interesting feature is that the peak of the response spectra shifts towards longer periods for larger magnitude events, reflecting the decreasing corner-frequency with increasing magnitude of a Brune (1970) source spectra. An oddity in the response spectra is the kink appearing between T = 3.5 s and T = 4.5 s , which is most likely from the sudden loss of approximately 3000 usable records (in regression) due to the high-pass filtering frequency criterion #5 in Ground-Motion Data and Selection Criteria. Similar kinks can be seen in a few other subsequent plots as well. Traditionally, GMM fixed-effect regression coefficients are smoothened in a post-processing step, in order to smoothen the predicted response spectra and to constrain the GMM behaviour in sparsely sampled magnitude-distance ranges, e.g. Abrahamson et al. (2014). However, in this study we chose not to alter the fixed-effect coefficients, and to maintain full consistency with their variance-covariance matrices; which will be necessary in testing and updating of the GMM. Figure 7 shows the random-effect and residual standard-deviations of the GMM. The total ergodic standard-deviation of the GMM is = √ 2 L2L + 2 0 + 2 S2S + 2 , when considering all source and site variabilities as aleatory. The solid lines in this plot correspond to the variance estimates of the GMM from ESM dataset (this study), while the dashed lines indicate those of the K16 GMM from RESORCE dataset. Note that the K16 GMM does not regionalise the events, and thus between-locality variability component L2L is absent. For the sake of comparison between the two models, we combine the two source related random-variances of this study into the traditional between-event variability =

Random-effects and residuals
We show that the total standard-deviation σ of the new GMM is considerably larger than that of K16 GMM. The largest increase in σ is at short-moderate periods, and is from the increase in between-event τ and between-site S2S variabilities. The increase in betweenevent variability can be attributed to the increased regional diversity of earthquakes, the increase of number of 4 ≤ M W ≤ 5.5 events from 164 in RESORCE to 778 in ESM, and the additional 76 events with 3 ≤ M W < 4 in ESM. Between-site variability is clearly from the increase in number of recording sites from 385 in RESORCE to 1829 in ESM. The residual variability of the new GMM, however, is equal or smaller than that of K16, despite the increase in regional diversity, the sample size, and the recording distance range from 300 to 545 km. In our view, this plot emphasizes the need to move from ergodic to partially non-ergodic region-and site-specific ground-motion predictions to inhibit the effect of increased σ on PSHA (Anderson and Brune 1999;Bommer and Abrahamson 2006;Kotha et al. 2017;Rodriguez-Marek et al. 2013). Figure 8 shows the region-specific c 3,r adjustments (coloured lines and ribbons) of the 45 regions for T = 0.01 − 8 s . Most of these curves lay within the ± c3 (red lines) bounds. As indicated by c3 , the region-to-region variability of apparent anelastic attenuation is the highest at short periods, and decreases gradually towards longer periods. High-frequency ground-motions are more susceptible to strong anelastic decay in the crust, which could be related to the crustal properties. Long period SAs are not effected as much by anelastic attenuation, therefore regional differences are relatively smaller, and c3 is smaller by ~ 50% at T ≥ 1 s.

Anelastic attenuation variability
In RESORCE dataset, K16 observed that short period SAs attenuate faster in Italy than in Turkey, which was observed earlier in NGA-West2 dataset by Boore et al. (2014), and confirmed later in Fourier domain by Bora et al. (2017). However, these observations were based on distinguishing regions by administrative boundaries and not geological or geophysical features. Since c3 estimated from the new regionalisation is not zero, we can assert that regional variability of anelastic attenuation exists. Of course, regions with fewer data have larger epistemic uncertainty (standard-error) on their δc 3,r , but the largest epistemic uncertainty is always lower than the aleatory τ c3 , and decreases with increasing data. Figure 9 shows the spatial variability of apparent anelastic attenuation as captured by the regionalisation model of Fig. 2. In this figure, red polygons identify regions with slower than pan-European average of anelastic attenuation c 3 i.e. δc 3,r > 0 , and δc 3,r < 0 for the faster attenuating blue regions. Regions with insufficient data, and thereby large epistemic uncertainty on their δc 3,r , are white in color, and given δc 3,r = 0 . There are few interesting features in these maps: 1. Regions with similar attenuation characteristics are spatially clustered, although this is period dependent. In general, regions characterized with high seismic activity (e.g. Italy and Greece) show strong attenuation compared to those with lower seismic activity (e.g. central Europe). A Bayesian clustering of regions with similar attenuation characteristics is proposed in the companion paper by Weatherill et al. (2020c); wherein, similar regions (r) are clustered (cl) and assigned cluster-specific δc 3,cl with lower standard-error and smaller cluster-specific τ c3,cl . The epistemic uncertainty (standard-error) on cluster-specific δc 3,cl is smaller than that of δc 3,r due to the accumulation of data from multiple regions; while the cluster-specific τ c3,cl is smaller than the overall τ c3 due to a smaller within-cluster variability of region-specific δc 3,r 2. The best-sampled regions are in central Italy with 5703 records from Northern and central Apennines W (West), 3438 records from Northern and central Apennines E (East). While attenuating faster than the pan-European average, there appears to be a strong contrast between these adjacent regions Regions with fewer ground-motion observations, thereby larger epistemic error on c 3,r , are more transparent and appear white. Note that the color scale is limited to vary between ±3 * c3 (T) for each period T

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

Source variability
Earthquake source variability is split in two random-effect groups: between-event ΔB 0 e = N 0, 2 0 and between-locality ΔL2L = N 0, 2 L2L . Between-event terms can be estimated for the recorded earthquakes, but are difficult to predict for prospective earthquakes because of their spatiotemporal variability. Even when correlated with stress-drop or a source parameter that can explain the relative differences in ground-motions, prediction of stress-drop for the next event is not yet possible. Traditionally therefore, betweenevent variability is considered purely aleatory. Between-locality random-effect is intended to quantify a portion of the between-event spatial variability into tectonic localities. Figure 10 shows the L2L l (T = 0.01 − 8 s) for 56 tectonic localities. As always, the standard-error (or epistemic uncertainty) of L2L l values is smaller than the L2L . L2L values are non-negligible, and the GMM fit improves (based on analysis of variance) with introduction of between-locality random-effect group. Therefore, we consider it an effective random-effect grouping. Although it is tempting to relate L2L l to some source parameter, it is preferably done in the Fourier domain in a follow-up study. Figure 11 maps the various tectonic localities (indexed l) color coded to their δL2L l (T = 0.1 s) and δL2L l (T = 1 s) values. In the top panel, corresponding to δL2L l (T = 0.1 s) , a clear difference between central Italy and north-western Anatolia can be seen. Apparently, earthquakes located in central Apennine region produce substantially lower short period ground-motions than those in north-western Anatolia. Similarly, there is an apparent distinction between the central Apennines and Po-plain earthquakes. The central Apennines tectonic locality contains the recent M6.5 Norcia earthquake (2016) and associated shocks, and several historic earthquakes in that region. The Po-plain locality contains data from the substantially stronger M6.45 Friuli earthquake (1976) and a few recent earthquakes. At a first glance, it may appear as if the spatial patterns are due to the predominant focal mechanisms in the regions, but the diversity of focal mechanisms within each region (especially among smaller events) dissuades this hypothesis. Spectral decomposition of the ESM dataset by Bindi and Kotha (2020) revealed a much larger stress-drop of the M6.45 Friuli earthquake compared to M6.5 Norcia earthquake-the δL2L l of regions containing these earthquakes also contrast similarly . Once again, correlating the random-effects obtained in response spectra domain (as here) to the earthquake parameters (such as stress-drop) Fig. 11 Variation of L2L l across the 56 tectonic locality (l) for T = 0.1 s (top panel) and T = 1 s (bottom panel). Blue polygons identify tectonic localities that produced earthquakes less energetic than the pan-European average, red polygons identify localities with more energetic earthquakes, and white polygons are localities close to the average. Tectonic localities with fewer ground-motion observations, thereby larger epistemic error on L2L l , are more transparent and appear white. Note that the color scale is limited to vary between ±3 * L2L (T) for each period T 1 3 may not hold well at short periods. We reserve further elaboration of random-effects with the Fourier version of the ESM dataset.
Between-event variability, now partially corrected for spatial variability through ΔL2L = N 0, 2 L2L , is quantified in the distribution ΔB 0 e = N 0, 2 0 . Customary checks for B 0 e,l include dependencies on magnitude, depth, and style-of-faulting. A few points to note here: 1. EMEC estimates of M W are used in the GMM regression. However, uncertainties in M W are ignored despite their know impact on the τ 0 estimates (Kuehn and Abrahamson 2017) 2. Depth in our case is hypocentral depth of the event. Buried ruptures are likely to produce stronger ground-motions than the exposed ruptures. This phenomenon is introduced in some of the NGA-West2 GMMs e.g., Abrahamson et al. (2014) modelled the events with depth-to-top-of-rupture z tor ≥ 20 km to have 0.75-3 times larger SAs than exposed ruptures, depending on the period 3. In this study, we use the typical style-of-faulting classification provided in the ESM dataset: Normal (NF), Thrust (TF), Strike-slip (SS), and Unknown (U). However, instead of introducing a style-of-faulting random-effect in the GMM regression, we queried the ΔB 0 e and ΔL2L distributions and found no systematic differences. This result is consistent with the fact that style-of-faulting factors cannot be quantified using isotropic factors with no azimuthal variations . Centroid Moment Tensor (CMT) solutions are however available for too few events to treat style-of-faulting with azimuth dependent factors as in Kotha et al. (2019). Even if available, the diversity of crustal structure across pan-Europe makes it difficult, without substantially more metadata (e.g. take-off angles) in the ESM dataset The left column of Fig. 12 shows B 0 e,l (T = 0s, 0.1s, 1s) trends with M W . No significant offsets or trends imply the magnitude scaling of the GMM sufficiently captures the magnitude dependence of SAs. The error bars show the median absolute deviance (MAD) of B 0 e,l within bins of size M W = 0.5 . MAD is a robust estimate of variance when normality of distribution is not necessarily satisfied within each bin. The MAD estimates appear to be magnitude dependent, indicating heteroscedasticity of τ 0 . However, we do not provide a heteroscedastic variance model now without investigating first its significance; given the small number of large magnitude events compared to small-moderate magnitude event sample (see Fig. 1 for data distribution).
The right column of Fig. 12 shows the B 0 e,l (T = 0.01s, 0.1s, 1s) trends with depth. We see no significant trend with depth. Although deeper events are likely to produce stronger ground-motions, there exists several shallower events with similarly large B 0 e,l , and correcting B 0 e,l for depth shows no remarkable reduction in τ 0 -as indicated by the bin-wise τ 0 values (error bars in right column of Fig. 12).

Site-response variability
As in K16 GMM, the KiK-net (Dawood et al. 2016) based Kotha et al. (2018) GMM, and Brooks et al. (2020) study on North-sea ground-motion data, we did not introduce a siteresponse scaling parameter in the fixed-effects for two reasons: (1) Only 419 of the 1829 sites have measured (and not inferred) V S30 available in the ESM, and using data from only the sites with measured V S30 leads to a large reduction in data for GMM regression; (2) investigation of site-to-site variability can be performed in a subsequent step. For instance, depending on the application, one can regress a relation between S2S s ∼ V S30 measured or inferred from topographic slope and geology (Thompson et al. 2014;Vilanova et al. 2018;Wald and Allen 2007) or relate S2S s directly to topographic slope and geology at a regional scale (Crowley et al. 2019;Weatherill et al. 2020a, b).
Given the importance of site-response in hazard and risk assessments, the complexities in finding a compromise between a site-response proxy, its availability at different regional scales, and the propagation of uncertainties from GMMs to risk assessments, we intend to perform a separate investigation. Nevertheless, a database of S2S s (T = 0.01 − 8s) for the 1829 ESM sites along with their robust regression weights and standard-errors is provided.
Note that we always use the measured V S30 in the analysis shown in this study, and never the inferred values. Early analyses of S2S s (T = 0.0, 0.1, 1.0 s) trends with measured V S30 and topographic slope are presented in Fig. 13. The left column of Fig. 13 plots the mean and MAD (robust standard-deviation estimate) of S2S s (T = 0.0, 0.1, 1.0 s) within ranges of V S30 coinciding with Eurocode8 site classes A, B, C, and D, and slope divided into 9 bins of equal width between 0.001 and 1.000 m m −1 . Trend with V S30 (not measured but) inferred from a V S30 ∼ slope correlation model of Wald and Allen (2007) is comparable to the trend with slope from which it is inferred in the right column of Fig. 13.
The correlation S2S s ∼ V S30 (measured only) at short-periods T = 0.0, 0.1 s is rather poor, as indicated by similar mean and MAD (the robust standard-deviation) for V S30 < 800 m/s in the left column of Fig. 13. It appears that short-period site-responses of soft rock, stiff, and soft soils (EC8 class B, C, D) in this dataset are not adequately distinguishable. However, it is interesting to see that the mean of EC8 class A 'rock' sites with V S30 ≥ 800 m/s is much lower than the rest, along with a considerably larger variability. Short-period linear soil-response of rock sites is known to be highly variable compared to softer soils, whose nonlinear soil-response may suppress the otherwise high variability from linear-only amplification (Bazzurro and Cornell 2004). On the contrary, the long-period site-response of rock sites is less variable than that of softer soils (e.g. 180 m/s ≤ V S30 < 800 m/s) . Similar inferences can be drawn from the S2S s ∼ slope plots in the right column of Fig. 13, where higher slopes are (usually) indicative of rock sites on steep hillsides, and lower slopes at softer sites located on flatter sediments. and topographic slope at site location (right column). Note that, only 419 sites have measured V S30 , compared to 1829 sites with slope. At short-moderate periods (T = 0, 0.1 s), the binned mean and MAD (red error bars) show a poor non-parametric trend, which improves at longer periods (T = 1 s). The blue line represents the parametric forms (Eqs. 6 and 7) fit to the data For completeness, along with a database of S2S s (T = 0.0 − 8.0 s) for the 1829 ESM sites, we derive a continuous empirical models for both S2S s ∼ V S30 (measured) and S2S s ∼ slope correlations seen in Fig. 13. We chose a quadratic functional form instead of the traditional piecewise linear function; as shown in Eqs. (6) and (7). Here, the measured V S30 is in m s −1 and slope in m m −1 , the regression coefficients g 0 , g 1 , g 2 are different for SR V s30 and SR slope , and change with period. The residuals from S2S s ∼ V S30 (measured) and Robust linear fits using an M estimator (Venables and Ripley 2002), at each of 34 periods between T = 0.01s − 8 s, PGA(T = 0 s) and PGV, are derived for δS2S s ∼ V S30 correlation of 419 sites with measured V S30 available, and δS2S s ∼ slope of the 1829 sites with slope derived from digital elevation models provided by Shuttle Radar Topography Mission.
Although heteroscedastic models of V S30 S2S and slope S2S appear reasonable in Fig. 13, we chose not to propose one without testing its significance; given the uneven distribution of sites in different bins. Figure 13 shows the reduction in between-site variance from using V S30 and slope as site-response proxies (Eqs. 6 and 7) in the GMM. For comparison, the variances of K16 GMM are also shown in this plot. Note that the K16 GMM comes in two variants: one without a site-response component, and another with measured V S30 as a proxy for linear site-response. Note that, the K16 standard-deviations shown in Fig. 7 are those when not using V S30 as a parameter, while those in Fig. 14 are when using V S30 as site-response proxy-hence, lower in the latter.
14 Between-site and total variance estimates ( S2S , ) of the GMM for T = 0.01-8 s (solid lines) compared with those of K16 GMM from RESORCE dataset (dashed lines). Reduction of ( S2S , ) to ( V S30 S2S , V S30 ) using V S30 , and to ( slope S2S , slope ) using Slope as site-response proxies is shown here. Note that the ( S2S , ) of K16 GMM are those using V S30 as site-response proxy, and are smaller than those shown in Fig. 7, which are from the K16 GMM version without a site-response parameter 1 3 A significant reduction in between-site standard-deviation can be achieved using an efficient site-response proxy or a combination of proxies. Since only a few sites are provided with measured V S30 values, V S30 S2S is substantially smaller than slope S2S . For a new site with V S30 or slope available, Eq. (6) or (7) can be appended to the ln( ) in Eq. (2), while replacing S2S (when estimating σ) with V S30 S2S or slope S2S , respectively. For site with neither site-response proxy available, but with sufficient strong ground-motion recordings, the site-specific S2S s term can be estimated and used for site-specific ground-motion predictions.

Aleatory variability
The last component of the GMM is the apparent aleatory variability, quantifying the natural randomness of the ground-motion data-that which is not captured by the mixedeffects. The aleatory residuals, = N 0, 2 , are tested for event depth and recording site distance dependencies in Fig. 15. For all periods ( T = 0.0, 0.1, 1.0 s ), we observe no significant trends in binned means and MAD, which implies the f R,g M W , R JB and f R,a R JB of the GMM (Eqs. 2 and 3) explain the distance scaling of ground-motion sufficiently well. Heteroscedasticity of is not evident either in these plots.
The moderate-long period residuals of this dataset show a clear evidence of anisotropic shear-wave radiation pattern in the near and intermediate distance ranges of R JB ≤ 80 km , similar to that reported by Kotha et al. (2019) with the KiK-net (Dawood et al. 2016;Kotha et al. 2018) and NGA-W2 (Ancheta et al. 2014;Boore et al. 2014) datasets and their GMMs. In addition, there is a strong evidence of SmS phases from Moho reflection in the 60 km < R JB ≤ 200 km distance ranges, e.g. Bindi et al. (2006), especially from events deeper than 10 km and close to Moho boundary. However, we chose not to discuss these features in the response spectral domain in which the GMM residuals are estimated, but instead with those in the Fourier domain as a follow-up.

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

Ergodic application
The first approach is by ignoring all repeatable effects, i.e. the region, locality, and site-specific adjustments. The between-locality, between-event, between-site, residual standard-deviations can be combined into an ergodic, total standard-deviation = √ 2 L2L + 2 0 + 2 S2S + 2 , as shown in Fig. 7. In this case, the regional differences in anelastic attenuation, quantified by c3 , will be treated as an epistemic uncertainty on far-source distance scaling. The epistemic uncertainty on the regionalised anelastic attenuation coefficient c 3 in Eq. (4) is c3 . This uncertainty can be handled with a GMM logic tree consisting of a slower c 3,r = c 3 + . c3 , average c 3,r = c 3 , and faster c 3,r = c 3 − . c3 attenuating branches. Following the Miller III and Rice (1983) three point approximation of a Gaussian distribution, with ≈ 1.732 one can use logic tree branch weights of 0.167, 0.666, and 0.167, respectively for slower, average, and faster branches. Consequently, the ground-motion prediction is a weighted mixture of three Gaussian distributions N ln( ), 2 , where ln( ) is estimated from Eq. (2) for three values of c 3,r ∈ c 3 , c 3 + ⋅ c3 , c 3 − ⋅ c3 .
Within the context of an ergodic application, if a site has either the measured V S30 or slope information available, but no site-specific ground-motion recordings, then the S2S in the estimation of = √ 2 L2L + 2 0 + 2 S2S + 2 can be replaced with V S30 S2S and slope S2S , whilst appending the fixed-effects of Eq. (2) with Eqs. (6) and (7), respectively. Figure 14 shows the consequent reduction of σ to V S30 and slope values, when using site-response proxies V S30 and slope, respectively. Consequently, the ground-motion predictions follow the mixed Gaussian distribution N ln( ) + SR V s30 , V S30 2 or N ln( ) + SR slope , slope 2 , where ln( ) is estimated from Eq. (2) for three values of c 3,r ∈ c 3 , c 3 + ⋅ c3 , c 3 − ⋅ c3 , and SR V s30 or SR slope are estimated from Eq. (6) or (7). The reduced aleatory variability from using a site-response proxy is beneficial until when enough ground-motion data can be collected at a site, and then S2S s for the new site can be estimated using equations provided in Kotha et al. (2017);Rodriguez-Marek et al. (2013); Sahakian et al. (2018); Villani and Abrahamson (2015), etc.

Region-specific application
For region-specific applications, the predictions can be upgraded with the region-specific anelastic attenuation and tectonic locality specific adjustments. Anelastic attenuation is regionalised by adjusting the generic coefficient c 3 in Eq. (4) with a region-specific value c 3,r , as in c 3,r = c 3 + c 3,r , where region r is decided by the location of the site. Since c 3,r are estimated from a smaller region-specific sample of ground-motion recordings, it should be treated as epistemically uncertain. For this purpose, standard-error on c 3,r are provided as well, and these are always smaller than c3 . Treating the SE c 3,r as uncertainty on mean of a normally distributed sample, the 95% confidence interval of c 3,r would be c 3,r ± 1.6SE c 3,r . Similarly, depending on the tectonic locality of the earthquake, the GMM predictions can be further regionalised by adjusting e 1 in Eq. (2) to e 1,l = e 1 + L2L l , where l identifies locality tectonic locality of the earthquake. Since L2L l is estimated from a smaller locality-specific ground-motion sample, the 95% confidence interval is bounded by L2L l ± 1.6.SE L2L l . A reduction of up to 10% in σ is achieved by dropping the L2L from aleatoric variance, resulting in a smaller r = √ 2 0 + 2 S2S + 2 , as shown in Fig. 16. Further reductions to Slope r are V s30 r possible if the site-response scaling with slope and V s30 (Eqs. 6 and 7) are used in predictions, respectively.
It is interesting to compare the total standard-deviation σ (dashed line in Fig. 16) of the regionalised, V s30 based K16 GMM with V s30 r of this study. At T ≤ 0.3s , the V s30 r is about 6% larger than its K16 counterpart; while at T > 0.3 s, V s30 r is as much as 14% smaller. Despite the enormous increase in data, it appears that regionalisation has helped in substantially curtailing the increase in apparent aleatory variability.

Region-and site-specific application
Partially non-ergodic region-and site-specific ground-motion predictions are possible for those sites with S2S s provided with this GMM or for new sites with sufficient groundmotion data. S2S s for the 1829 sites in the ESM dataset are provided, along with the c 3,r Fig. 16 Reduction in total standard deviation estimates (σ) of the GMM for T = 0.01-8 s (solid lines) from ergodic with and without site-response proxy, region-specific with and without site-response proxy, and region-and site-specific values. Standard deviation estimates when using site-response proxies are indicated by the annotations with corresponding superscript, i.e. ( V S30 ) with V S30 , and to ( slope ) with slope. Annotations with subscript r correspond to the variances for regionalized predictions, i.e. from discounting region-to-region variability L2L from σ. Values annotated with subscript r, s correspond to region-and sitespecific predictions discounting also the S2S . The black dashed line is the total standard-deviation (same as in Fig. 14) of the regionalized and V s30 based Kotha et al. (2016) GMM, and compares with the V S30 r of this study 1 3 of their region, and L2L l of nearby tectonic localities (earthquake sources). Since betweensite and between-locality variabilities do not apply for site-specific predictions, the reduction in apparent aleatory variability is enormous, i.e. r,s = √ 2 0 + 2 is about 40% smaller than = √ 2 l2l + 2 0 + 2 S2S + 2 , as shown in Fig. 16. However, the reduction in aleatory variability will be accompanied by additional epistemic uncertainty. In addition to those in region-specific predictions, uncertainty on (the mean) S2S s should be accounted with ±1.6.SE S2S s . Region-and site-specific predictions therefore are a mixture of 27 Gaussian distributions N ln r,s , 2 r,s , where ln r,s is estimated from Eq. (1) for three values of c 3,r ∈ c 3 + c 3,r , c 3 + c 3,r + 1.6SE c 3,r , c 3 + c 3,r − 1.6SE c 3,r , three values of e 1,l,s ∈ e 1,l + S2S s , e 1,l + S2S s + 1.6SE S2S s , e 1,l + S2S s − 1.6SE S2S s , wherein e 1,l ∈ e 1 + L2L l , e 1 + L2L l + 1.6SE L2L l , e 1 + L2L l − 1.6SE L2L l .
In complement, a more practical application is presented in the companion study by Weatherill et al. (2020c), where a more thorough comparison of this GMM with contemporary models is provided. Weatherill et al. (2020c) also provides a more complete discussion on how the various epistemic uncertainties in this GMM can be handled in the logic tree framework of a PSHA.

Towards non-ergodic ground-motion predictions
For tectonic localities, attenuating regions, and sites with sufficient amount of recordings the epistemic uncertainty on the random-effect adjustments are negligible with respect to the random-effect and standard-deviations. Collecting more ground-motion recordings is principal in moving towards non-ergodic predictions. The benefits in resolving the ergodic assumption and progressing towards region-and site-specific in ground-motion prediction is demonstrated Fig. 17. In this plot, predictions for the M6.5 Norcia event of the central Italy sequence, which occurred on 30th October, 2016, are compared to the response spectra recorded at three sites covered by the Italian strong motion network (Gorini et al. 2010). These sites are identified by the network code IT in the ESM dataset: (1) permanent, freefield station LSS (Leonessa) with V S30 = 1091 m/s located 25 km from the event epicentre, (2) permanent, free-field station MVB (Marsciano Monte Vibiano) with V S30 = 1046 m/s located 65 km from the event epicentre and, (3) permanent, free-field station PSC (Pescasseroli) with V S30 = 1000 m/s located 110 km from the event epicentre. The three columns in Fig. 17 correspond to the three stations.
These event and stations are selected to demonstrate progressively (in Fig. 17) the impact of moving from ergodic prediction relying on V S30 as site-response proxy (top row), through region-specific predictions (middle row) considering regional (Northern and central Apennines West) anelastic attenuation ( c 3,r = c 3 + c 3,r ) and adjustment specific to the tectonic locality ( e 1,l = e 1 + L2L l ) containing the event (locality ID: "PTTC007"), to region-and site-specific predictions (bottom row) from an additional site-specific adjustment ( e 1,l,s = e 1 + L2L l + S2S s ). Both the median prediction and standard-deviation change in process, which is reflected by the width of the coloured ribbon in Fig. 17. The S2S s (T = 0.01 − 8s) of these sites are estimated from 29, 15, and 20 records from predominantly small-moderate earthquakes (details in the Fig. 17 panels). A few comments on this figure: 1. The ergodic median predictions (central line) and one V S30 interval (ribbon) are systematically above the observed response spectra at the three rock sites, located at near 1 3 (25 km), intermediate (65 km), and far (110 km) source distances. This is likely because the M6.5 Norcia event produced relatively weaker ground-motions compared to other large magnitude events recorded in Greece and Turkey-as quantified by their respective B 0 e,l values. Since the ergodic predictions consider all event, region, site, and record variabilities as aleatory, the V S30 (Fig. 16) is large yet not large enough to contain the M6.5 event observations within ± V S30 boundaries. 2. Region-specific ground-motion predictions for these sites are achieved by adjusting the GMM with the c 3,r of the Northern and central Apennines (West), in which the sites are located; and the L2L l of the tectonic locality PTTC007, in which the M6.5 event occurred (along with a few other prominent events and aftershocks). The epistemic uncertainties of these adjustments are relatively very small given the large number of recordings. In the middle row of Fig. 17, we notice the observed response spectra are closer to region-specific predictions than to the ergodic predictions. c 3,r and L2L l of Fig. 17 Comparison of ergodic (top row), region-specific (middle row), and sites-specific (bottom row) ground-motion median (solid colored line) and variability (colored ribbon) predictions for the M6.5 Norcia earthquake with the observed response spectra (solid red line) at three sites in Italy (column wise) these regions are both lower than the pan-European average (which is zero), meaning the region attenuates short-period ground-motions faster and the events on average produce weaker ground-motions than elsewhere in pan-European region. It is interesting to note that, at short distances (site IT_LSS) the c 3,r has no effect on region-specific predictions, and the shift is mostly from L2L l -and so is at intermediate distance (site IT_MVB). At far-source distances (site IT_PSC), the combined effect of c 3,r and L2L l worked well to capture the observed response spectra within the narrower ± V S30 r ( Fig. 16) range about the region-specific median. 3. Region-and site-specific predictions (site-specific in short) for the three sites are shown in the bottom row of Fig. 17. Along with the curves, details on the number of recordings, magnitude and distance ranges (1st and 3rd quantile) of the recordings are provided. The additional adjustment to the preceding region-specific predictions is through e 1,l,s = e 1 + L2L l + S2S s . While most of the data in estimating the S2S s of these sites is from small-moderate sized earthquakes, the site-specific predictions fit quite well with the observations for the large M6.5 event. Since S2S s are used to adjust the region-specific ground-motions, the V S30 becomes irrelevant, and resulting in a smaller r,s .
In the above example, we demonstrated that applying region-specific adjustments noticeably improved the match between observations and predictions. The best agreement was clearly from using region-and site-specific adjustments. To substantiate this claim, we performed a 10-fold cross-validation exercise to verify if the introduction of the various random-effects into GMM functional form indeed improves its prediction capabilities. In doing so, we re-run the regression with three functional forms: • Ergodic model with no regionalisation of anelastic attenuation, no localisation of event terms, and no site-specific adjustments, i.e. no random-effects c 3,r , L2L l , S2S s • Regional model with regionalisation of anelastic attenuation and tectonic localisation of events, leaving out the random-effect to capture site-to-site variability, i.e. only the random-effects c 3,r and L2L l , but no S2S s • Site-specific model identical to the GMM presented here i.e., regionalisation of anelastic attenuation, localisation event terms, and the site-to-site random-effect To perform the cross-validation, the dataset is split into 10 parts with non-overlapping events. Meaning, earthquakes (and records) are exclusive to their subsets and do not feature in any other subset. We perform regression of the three models on any nine subsets combined, and test the predictions on the tenth subset. Root-mean-squarederrors are estimated for each trial and then averaged over the ten trials. This exercise is repeated for IMs:PGA, PGV, SA(T = 0.1, 0.2, 0.5, 1, 2, 4 s). Figure 18 shows a histogram of RMSE for each IM, from the three models. The reduction in RMSE from ergodic to region-specific GMMs is clear and most prominent at short periods. Towards longer periods, the improvement is less pronounced but is still substantial. This is because the regionalisation random-effects are focused on the capturing variabilities in short-period ground-motions. In Fig. 7, we notice the overall variability σ peaks at short periods. In Figs. 8 and 10, the largest regional variability of anelastic attenuation c3 and tectonic localities L2L is also at short periods T ≤ 0.3 s . Without regionalisation, along with the increased the short period ground-motion prediction variability and reduced precision of an ergodic median, the predictive capability (measured as RMSE) of the GMM is reduced as well.
Across all the periods, i.e. the entire response spectra, the best predictive capabilities are those of region-and site-specific GMM. In Fig. 16, the largest reduction in σ is achieved not from using V S30 or slope as the site-response proxy, but from using S2S s as the sitespecific adjustment. However, it is unlikely that every site has sufficient ground-motion data to estimate it's site-specific S2S s . In that case, alternative site-response proxies are sought to predict the S2S s , as in Kotha et al. (2018); Weatherill et al. (2020b). However, even in these studies, while the long period site-response could be partially explained using some geotechnical parameters, short-period site-response is much more variable-even among the so-called reference rock sites (Bard et al. 2019;Pilz et al. 2020).

Summary
In this study, we present the update of the Kotha et al. (2016) region-and site-specific GMM, using the recently published European Strong Motion dataset (ESM) by Lanzano et al. (2019a) and Bindi et al. (2018b). The update derived from a dataset 15 times larger than the older RESORCE dataset (Akkar et al. 2014b), with only minor changes to fixedeffects formula, and one additional random-effect grouping. Leveraging on the recently developed regionalisation datasets, TSUMAPS-NEAM by Basili et al. (2019) and the tectonic zonation defined in the purview of European Seismic Hazard Model 2020 project, we have regionalised the ground-motion data to capture spatial variation of anelastic attenuation and earthquake characteristics, respectively.
Due to the increased spatial and temporal diversity of the data, the random-effect and residual variances are larger than the older pan-European GMMs. This increased variance motivates the development of region-dependent ground-motion models since a larger ergodic σ will severely affect the probabilistic seismic hazard and risk assessments. We propose making use of the attenuation region, tectonic locality, and site-specific terms (random-effects), along with their epistemic uncertainty, to develop partially non-ergodic ground-motion predictions with a reduced σ. The improved predictive power of the regionand site-specific version of the GMM is substantiated with a 10-fold cross-validation. All the random-effect values will be provided on request, along with the fixed-effects Fig. 18 10-fold cross validation statistics comparing the predictive capabilities of ergodic, regionalized, and site-specific GMM versions at various periods coefficients, and the variance-covariance matrices needed to update the model in a Bayesian framework.
Only the customary random-effect and residual analyses are presented here, leaving out their physical evaluation for a follow-up study. A similar exercise is actually carried out on the Fourier amplitude version of the ESM dataset. Random-effects and residuals of the Fourier GMM are easier to associate to geophysical and geological parameters that are not involved in the GMM regression. If indeed the spatial variability of source, path, and site random-effects can be attributed to a globally available parameter, e.g. a crustal velocity map, we can attempt migrating the GMM developed from seismically active regions to the less active regions with sparse ground-motion recordings. The GMM is developed with the intention of perpetual updates, when and wherever new datasets become available. Further elaboration on random-effect analyses, GMM update and application procedures are ongoing, and will complement this study. Meanwhile, this new GMM has the widest applicability yet for the pan-European region; as it is derived from a dataset stretching from Pyrenees in the west to Iran in the east, Rhine-Graben in the north to Hellenic arc in the south, and constituted of manually processed ground-motion data shallow crustal earthquakes of 3 ≤ M W ≤ 7.4 recorded at 0 ≤ R JB ≤ 545km.

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