Checking the site categorization criteria and amplification factors of the 2021 draft of Eurocode 8 Part 1–1

In this paper the site categorization criteria and the corresponding site amplification factors proposed in the 2021 draft of Part 1 of Eurocode 8 (2021-draft, CEN/TC250/SC8 Working Draft N1017) are first introduced and compared with the current version of Eurocode 8, as well as with site amplification factors from recent empirical ground motion prediction equations. Afterwards, these values are checked by two approaches. First, a wide dataset of strong motion records is built, where recording stations are classified according to 2021-draft, and the spectral amplifications are empirically estimated computing the site-to-site residuals from regional and global ground motion models for reference rock conditions. Second, a comprehensive parametric numerical study of one-dimensional (1D) site amplification is carried out, based on randomly generated shear-wave velocity profiles, classified according to the new criteria. A reasonably good agreement is found by both approaches. The most relevant discrepancies occur for the shallow soft soil conditions (soil category E) that, owing to the complex interaction of shear wave velocity, soil deposit thickness and frequency range of the excitation, show the largest scatter both in terms of records and of 1D numerical simulations. Furthermore, 1D numerical simulations for soft soil conditions tend to provide lower site amplification factors than 2021-draft, as well as lower than the corresponding site-to-site residuals from records, because of higher impact of non-linear (NL) site effects in the simulations. A site-specific study on NL effects at three KiK-net stations with a significantly large amount of high-intensity recorded ground motions gives support to the 2021-draft NL reduction factors, although the very limited number of recording stations allowing such analysis prevents deriving more general implications. In the presence of such controversial arguments, it is reasonable that a standard should adopt a prudent solution, with a limited reduction of the site amplification factors to account for NL soil response, while leaving the possibility to carry out site-specific estimations of such factors when sufficient information is available to model the ground strain dependency of local soil properties.


Introduction
In modern seismic codes, the actions for design are generally based on a seismic hazard assessment at the reference ground. With few exceptions, the latter is assumed at the standard outcropping seismic bedrock, usually defined by shear wave velocity V s larger than a threshold value and flat topography. Site amplification factors are then prescribed to modify the elastic design spectrum, based on a categorization of ground types. In case ground conditions are not compliant with the standard ground types, seismic codes often require numerical analyses of seismic site response, which are typically based on 1D models of the specific soil profile.
To take into account the variety of natural soils in a simplified but comprehensive way, the site categorization of a seismic norm should target several, partly conflicting, objectives: (i) it should group soils into clusters as homogeneous as possible to achieve similar site amplification levels; (ii) it should include the largest variety of realistic soil profiles; (iii) it should be based on a limited number of proxies, easy to be measured and capable to identify with no ambiguities the different categories. Obviously, the best candidate proxies are those mostly related to seismic site amplification, such as: the shear wave velocity of the soil deposit; its thickness; the impedance ratio with the bedrock; or a combination of the previous parameters such as the fundamental frequency of the soil deposit f 0 .
Part 1 of Eurocode 8 (CEN 2004, denoted in the following by EC8-1) introduced a relatively simple and effective site categorization (Table 1), mainly based on V s,30 , i.e., the time-averaged shear wave velocity in the top 30 m, but with the possibility to consider, as alternative parameters, the number of blows in the Standard Penetration Test (N SPT ) and the undrained shear strength (c u ). V s,30 was firstly proposed by Borcherdt (1994), after the study on site effects during the Loma Prieta earthquake and adopted later as a proxy for site amplification by several seismic codes, besides EC8-1, such as the BSSC (2015) provisions, the ASCE/SEI 7 (ASCE 2010) and the New Zealand Standard (NZS 2004). Despite its extensive use, the ability of V s,30 to capture the complexity of site effects has been questioned by many authors (e.g., Castellaro et al. 2008;Gallipoli and Mucciarelli 2009;Lee and Trifunac 2010). As a consequence, several alternative classification schemes have been suggested, based on a combination of different proxies (e.g., Rodriguez-Marek et al. 2001;Bouckovalas and Papadimitriou 2003;Luzi et al. 2011). More recently, Pitilakis et al. (2019) proposed a refined version of their multi-parameter classification scheme (Pitilakis et al. 2013), along with the associated amplification factors, that adopts as main proxies the fundamental period of the soil profile, the thickness and the average shear wave velocity of the soil deposit down to the "seismic bedrock", together with a description of the soil deposit.
Several limitations have often been pointed out referring to the EC8-1 classification scheme: 1) the ambiguity of the qualification of "several tens of m in thickness" of the soil deposit to identify ground types B and C (see Table 1); 2) the classification may not be unique, depending on whether different proxies are used (for example, the value of V s,30 may identify a given ground type, that may be in contradiction with the identification based on N SPT ); 3) in case the variation of the soil profile does not fit "a gradual increase of mechanical properties with depth" (see ground type B description   Table 1), such as when velocity inversions are present, it is not clear what to do; 4) the alternative proxies to V s,30 (N SPT and the undrained strength c u ) are related to soil strength, and only poorly related to the dynamic soil response; as a consequence, a soil profile might easily be misclassified by using such parameters alone.
To cope with such limitations, the introduction of alternative (or complementary) proxies to V s,30 , such as f 0 , was often proposed, especially taking advantage of the improved non-invasive techniques for its estimation (Comina et al. 2011). However, all of these arguments should face the trade-off between accuracy of the approach and easeof-use, a typical dilemma for a technical norm. This is especially true for the problem at hand, since the seismic site amplification effects do not depend on the site properties alone, but also on the way that such properties interact with the input motion, including its frequency content, amplitude and direction of propagation. Therefore, the site amplification factors (SAF) associated to a given site category are inevitably subjected to large dispersion. Reducing such dispersion may require a larger number of site categories, but this would be in contrast with the necessary "ease-of-use" and with the principle of Occam's razor.
The empirical estimation of SAFs from earthquake records is getting more and more reliable thanks to several factors such as: the continuously increasing number of ground motion records; the improved site characterization of recording stations; and the development of the techniques for statistical analysis of the station-to-station residuals (Stafford 2014). However, some "extreme" ground types such as very hard or very soft soils are still relatively poorly sampled in the worldwide record datasets, especially in the case of large intensity ground motions recorded in the near-source region of large earthquakes.
Being rooted on research studies dated back to twenty years ago, EC8-1 site amplification factors were established based on a relatively limited amount of records from European earthquakes (Rey et al. 2002). Specifically, only few tens velocity profiles were available, with limited sampling especially for reference rock outcrops and for very soft soil conditions. According to current EC8-1, site amplification is determined based on the coupled effect of the ground type dependence of the site parameter S (affecting the whole spectral period range) and of the corner periods T B and T C (Fig. 1). More specifically, the dependence of T C on the site category has the role of widening the spectral shape towards longer periods as the ground type becomes softer and thickness of the soil deposit increases. Besides, discrimination of high seismicity (Type 1) and low seismicity (Type 2) spectral shapes, and of the corresponding S-factors and corner periods, allows to implicitly account for the intensity dependence of site amplification.
The EC8-1 site categorization scheme and the corresponding SAFs proved to be a rationale and effective solution for the simplified quantification of seismic site amplification effects on the elastic design spectra. However, because of the previous limitations, the EC8-1 criteria were reassessed within the comprehensive work in progress by the Committee CEN/TC250/SC8 "Eurocode 8; Earthquake resistance design of structures" for drafting the new version of Eurocode 8. This implied a new definition of elastic design spectra, where spectral amplitudes are introduced based on two parameters related to seismic hazard, i.e. the constant acceleration branch at short periods (S α,RP ); and the spectral acceleration at T = 1 s (S β,RP ). These parameters are to be evaluated for the reference outcropping rock conditions (shear wave velocity V s > 800 m/s) and for the return period RP. Correspondingly, the site amplification factors F α and F β were also introduced to multiply S α,RP and S β,RP , related to the site amplification at short and intermediate periods, namely: F α = SAF of the constant acceleration spectral plateau at short periods; F β = SAF of the T = 1 s spectral ordinate.
After introduction of the criteria for site categorization and for the evaluation of the corresponding site amplification factors implemented in the current draft of EC8 (in the following denoted by "2021-draft", CEN/TC250/SC8, 2020), in the subsequent parts of the paper a validation of the proposed SAFs is illustrated, taking advantage of two complementary approaches.
First, an empirical approach is adopted, in which a wide dataset of strong motion records is built. Records are downloaded from the European Engineering Strong Motion database (ESM, https:// esm. mi. ingv. it/;Luzi et al. 2016a, b), and from the Japanese KiK-net database (http:// www. kyosh in. bosai. go. jp). Using the available site information, the recording stations are classified according to the 2021-draft criteria, while the empirical spectral amplifications are estimated computing the site-to-site residuals with respect to the reference rock conditions according to the ground motion prediction model used to derive the new European seismic hazard maps .
Second, the SAFs of the 2021-draft are compared with the results of a comprehensive parametric numerical study of 1D site amplification based on randomly generated shear-wave velocity profiles, classified according to the 2021-draft criteria.
Finally, the non-linear correction factors to site amplification, introduced within the 2021-draft, are checked through the detailed analysis of ground response at several Kik-Net stations, i.e., some of the few worldwide stations that recorded both weak and strong motions, at both ground surface and borehole levels. In this way, it is possible to obtain a quantitative estimation of the non-linear reduction effects on the response spectral ordinates.

Site categorization
Since it was recognized that V s,30 alone is not sufficient to encompass the variety of site amplification effects in real geological profiles, in the 2021-draft the decision was taken to explicitly introduce two proxies for site categorization. The two selected parameters are H 800 and V s,H (Table 2): H 800 is the depth of the seismic bedrock formation, which is identified by V s higher than 800 m/s; Comparing Tables 1 and 2, the 2021-draft introduces explicitly a further "F" category for deep soil deposits (H 800 > 100 m) and slightly changes the ranges of variability of V s,H associated to stiff, medium-stiff and soft ground classes. Furthermore, the maximum sediment thickness associated to the E category increases from 20 to 30 m. Guidance is provided in the 2021-draft for those cases where a direct evaluation of V s,H and H 800 is not available and for sites where, in the absence of direct measurements of V s , empirical correlations with other geotechnical or dynamic parameters can be used. More specifically a categorization criterion, alternative to Table 2, is introduced for sites where H 800 is not known and instrumental measurements of the fundamental site frequency f 0 are available (Table 3). In Fig. 2, the performance of this alternative scheme is tested against a set of accelerometric stations of the Italian and Swiss seismic networks, for which both the V s profiles and the instrumentally estimated f 0 values are available Swiss Seismological Service (SED) at ETH Zurich 2015). As shown in Fig. 2, the performance of the f 0 -based classification scheme is in reasonable agreement with that based on H 800 , given that both schemes require V s,H as the additional classification proxy.

Site amplification factors
Aiming at smoothening the strong jumps of SAFs between one category and the adjacent one, a continuous variability for F α and F β is introduced in the 2021-draft, as shown in Table 4. Alternatively, default values are introduced for those cases where a simplified    Table 3, checked according to the available V s profiles and the instrumentally calibrated f 0 values of 57 accelerometric stations of the Italian (43 stations) and Swiss (14 stations) seismic networks (Aimar 2018) 1 3 identification of the site category is adopted because of insufficient information. F α and F β mainly depend on the ratio V s,H /800, according to a functional form that is often used to quantify site effects in empirical ground motion prediction models (see e.g., Cauzzi et al. 2015;Lanzano et al. 2019a). Besides, the empirical factors r α and r β introduce a dependence of F α and F β on the amplitude of the hazard related spectral ordinates S α,RP and S β,RP . An example of the SAFs obtained according to Table 4 is given in Fig. 3, for two pairs of S α,RP and S β,RP , representative of high and low seismicity conditions. The continuous expressions in Table 4, including the r α and r β factors, were set to approach, albeit in a simplified way, the period-dependent SAFs according to the Ground Motion Model (GMM) by Akkar et al. (2014), based on European and Middle-East records. This GMM was considered because it provides amplitude-dependent SAFs, so that it complies with the objective of the 2021-draft to differentiate SAFs according to ground motion intensity.
In Fig. 4 (left side) a comparison is shown between the SAFs introduced by Akkar et al. (2014), for two representative values of PGA, and those of two further GMMs (Lanzano et al. 2019a;Kotha et al. 2020) recently introduced in the framework of the research studies to support the new seismic hazard map for Europe, in progress within the European Project SERA (http:// sera-cms. ethz. ch/). Note that the latter GMMs do not account for amplitude variability of SAFs. For all considered GMMs the reference is given by V s,30 = 800 m/s. In the right-hand side of Fig. 4 the ratio of elastic response spectra with respect to the ground type A, both according to EC8-1 (Type 1 and Type 2) and to 2021-draft are also shown. Some comments on such comparison follow: • The trend of variability of SAFs with both period and V s,30 according to the GMM is reasonably well reproduced, taking into account the simplified form of the 2021-draft factors; also, application of the r α and r β factors in Table 4 establishes a good agreement with the amplitude-dependent SAFs of Akkar et al. (2014); • There is an overall good agreement of response spectral ratios between EC8-1 and 2021-draft for the high seismicity areas (Type 1 spectra of EC8-1), with the single exception of V s,30 = 200 m/s. As a matter of fact, this corresponds to a value slightly above the boundary between the adjacent ground categories C and D of EC8-1 (the  Kotha et al. (2020) and Lanzano et al. (2019a). Right: ratios of elastic design spectra with respect to ground type A according to EC8-1 (Type 1 and Type 2) and to 2021-draft, considering typical values for low and high seismicity regions. Values of 2021-draft refer to a high seismicity (S α,RP = 6 m/s 2 , S β,RP = 2 m/s 2 ) and a low seismicity site (S α,RP = 2 m/s 2 , S β,RP = 0.4 m/s 2 ). Site amplification factors of 2021-draft were computed for the intermediate depth class threshold being at 180 m/s). Therefore, while according to EC8-1 the site belongs to class C, according to the continuous formulation of Table 4 the corresponding SAF gets closer to class D; • For Type 2 spectra, the long period amplification observed from strong motion records and reflected by the GMMs is underestimated by the EC8-1 spectral ratios. Instead, 2021-draft provides a better agreement; • The most recent GMMs by Kotha et al. (2020) and Lanzano et al. (2019a) tend to provide larger amplification factors for low values of V s,30 , especially for short periods, probably because of a more careful identification of site conditions at recording stations.
Note that the comparison in Fig. 4 is valid for intermediate soil depth 30 m < H 800 < 100 m, since the considered GMMs do not distinguish shallow or deep soil sites (category E or F). For ground category F, the standard values in Table 4 for intermediate depth (categories B, C, D) are modified by a reduction factor 0,9 at short periods and an amplification factor 1,25 at long periods. These values are in reasonable agreement with the empirical site amplification factors by Lanzano et al. (2016), that are based on records from the Po Plain, Northern Italy, and with numerical modelling of site response of deep deposits in the Central Us, presented by Hashash and Moon (2011).
For ground category E, the SAFs in Table 4 were obtained by introducing a simple continuous expression, as a function of V s,H and H (coincident in this case with H 800 ), that could provide amplification at short periods with no jumps with the adjacent site categories. Validation against records and comparison with 1D simulations in the next sections provides a check of such approximation.

Datasets and site categorization
To make an empirical validation of the site amplification factors in the 2021-draft, we exploited records, stations and events information contained in two accelerometric databases, i.e., the ESM and the KiK-net. Note that, while the ESM records are in the  Table 4, KiK-net records are not, so that a proper independent validation can be achieved in the second case. For both databases, flat-files of ground-motion parameters and metadata are publicly available. Namely, the ESM flat-file (Lanzano et al. 2019b) consists of metadata and intensity measurements of about 23,000 waveforms of 2,179 pan-European events with moment magnitude greater than 4.0, recorded in the time span 1969-2016 by 2,080 accelerometric stations. The KiK-net flat-file (Dawood et al. 2016) contains metadata and intensity measurements of about 108,000 recordings of 3,110 Japanese events with moment magnitude greater than 3.5 from 1997 to 2011, recorded by 692 accelerometric stations. The magnitude-distance distribution of the two datasets is reported in Fig. 5: in both cases, data are well sampled in the magnitude range 3.5-6.5 and for epicentral distances up to 200 km. Although the bulk of data is related to moderate and low earthquakes, there is also a significant number of near-source records (e.g., with epicentral distance less than 15-20 km) related to strong events with magnitude in the range between 6.0 and 7.6, that is expected to constrain reliably results in that magnitude and distance range (Lanzano et al. 2019a, b).
To classify the recording sites according to the 2021-draft (Table 2), we selected only stations for which both V s,H and H 800 are available. This corresponds to 212 stations, out of 2,080, for the ESM and to 655 stations, out of 692, for the KiK-net flat-file. Figure 6 shows the distributions of the characterized stations in the six soil categories of the new EC8-1 proposal for the two flat-files.
Tables 5 and 6 show the migration matrices between EC8-1 and 2021-draft site categorization for the ESM and the KiK-net datasets, respectively. Among the total 212 stations of the ESM flat-file almost half of them are classified as B (47%), 50 stations as A (23.5%), while only the remaining 63 are classified in softer site categories such as C (7.5%), E (10%), F (10%) and D (less than 2%). This result confirms that accelerometric stations installed in Europe are mainly located in soft rock sites or stiff soils, that are the most representative site conditions in European countries with high seismicity and with the largest number of installed stations with availability of records, such as Italy (Forte et al. 2019). The situation is different for the Japanese KiK-net stations, the site categorization of which denotes the dominant presence of either shallow (class E, 37%) or deep soil deposits (F, 16%), or stiff soils (B, 27%), while there is a poor presence of class C (10%), class D (3%) and especially of class A soils (5%) that makes the evaluation of site amplification with respect to reference site conditions difficult to be properly constrained.
Finally, it is worth noticing that site category A of the 2021-draft is, with few exceptions, a subset of site category A for EC8-1: indeed, when H 800 < 5 m, sites with V s,H < 250 m/s are now classified as E (as it is the case for many Japanese stations, formerly classified as A), while, when 5 m < H 800 < 30 m, sites that were previously classified as A because V s,30 > 800 m/s may be classified differently because of V s,H < 800 m/s. This confirms that the definition of the reference class A in the 2021-draft is, in general, more restrictive than in the EC8-1.

Residual analysis with respect to reference site conditions
The selected approach is based on the residual analysis of ground motion records, aiming at a statistical evaluation of the distance of records, belonging to a given site category, from the predictions obtained by a ground motion model (GMM) at the reference site conditions for the corresponding magnitude and source-to-site distance. More specifically, following the Al Atik et al. (2010) scheme, the total residuals (R es ) with respect to the median GMM are decomposed as follows: where the subscripts e and s refer to events and stations, respectively. δB e represents the event-term (between-events residual), which corresponds to the average misfit of recordings of one particular earthquake with respect to the median ground-motion model; δS2S s (1) R es = B e + S2S s + W 0,es is the site term, that quantifies the average misfit of recordings from one particular site with respect to the event-corrected median ground-motion; δW 0,es is the remaining residual after site-and event-terms are subtracted from total residuals (Lanzano et al. 2017). To have a more robust calculation of the residual terms, the between-event, site-to-site and event-and site-corrected components are computed by the random-effect model (Stafford 2014).
In this study, we focused on δS2S s term, that is directly related to the site response of the recording station. In our application, we computed these site terms with respect to reference site condition of the selected GMM (re-named δS2S A ), since they represent the empirical amplification function of the station (see e.g., Priolo et al. 2019;Lanzano et al. 2020).
For the aim of this work, comparison of δS2S A with the SAFs of the 2021-draft is meaningful if two conditions apply: 1. The reference site conditions should be the same for both the GMM model and the SAFs; 2. Both GMM and SAFs should be consistent with the reference conditions of the seismic hazard model that the seismic norms are based upon.
For this reason, we considered, as the reference for the residual analysis, the GMM introduced by Kotha et al. (2020) for Europe (denoted by KOT20) and employed as a backbone predictive model of shallow crustal earthquakes for the update of the seismic hazard maps of Europe (ESHM20, Weatherill et al. 2020). Although each European member country will be required to provide its own seismic hazard map in terms of a National Annex, the ESHM20 maps will be in any case an important reference in the background of the upcoming new version of Eurocode 8.
Note that the ESHM20 maps are defined for the reference condition V s,30 = 800 m/s. This is consistent with the 2021-draft that, as shown in Table 4, prescribes SAFs larger than 1 to be applied only for site categories other than A. Therefore, all sites belonging to category A are conventionally assumed by the 2021-draft, as by the EC8-1 as well, to have the same unity amplification level as the rock site with V s,30 = 800 m/s. As a consequence, for those sites of category A harder than V s,30 = 800 m/s, the 2021-draft makes an "implicit" over-estimation of seismic hazard, and consequently of seismic actions for design, essentially aiming at ease-of-use. A quantification of such overestimation will be illustrated in the next section (Fig. 7).

Comparison with the site amplification factors of 2021-draft
Consistent with the selected GMMs, records from the ESM and KiK-net datasets were selected from shallow active crustal earthquakes with focal depth less than 35 km and source-site distance within 200 km. This caused a reduction of the number of records to 17,420 and 23,307, corresponding to 1,726 and 933 events for ESM and KiK-net, respectively. The response spectral ordinates were then computed as the geometric mean of the two horizontal components.
Furthermore, when information was available, we removed stations with potential topographic amplification effects. Finally, for a more robust estimation of δS2S A , only stations with at least 5 records are considered. At the end, the number of stations used for the comparison decreased to 114 and 599 for ESM and KiK-net, respectively. First of all, we computed the median of the site terms δS2S A (μ δS2S-A ) associated with stations classified as site category A for both ESM and KiK-net (Fig. 7).
The δS2S A was estimated for 37 and 27 recording stations, classified as A in ESM, according to EC8-1 and 2021-draft, respectively; while for KiK-net, the recording stations are 59 (category A in EC8) and 27 (category A in 2021-draft). The trends of median sites terms are very similar for both A classes in ESM: they are almost flat with amplitude around -0.5 (corresponding to a de-amplification of a factor equal to 0.6 with respect to   Table 5. However, it also implies that the level of seismic motion at the "median" site category A is over-estimated with respect to the level provided by the V s,30 = 800 m/s seismic hazard, as discussed in the previous section. In the case of KiK-net (right-hand side of Fig. 7), the differences between EC8-1 and 2021-draft are remarkable at short periods: in particular, the EC8-1 median trend exhibits slight or even zero de-amplification at short periods. To justify this trend, it should be recalled that, as shown in Table 6, many KiK-net stations migrated from site category A (EC8-1) to E (2021-draft). Therefore, while the median trend of site category A (EC8-1) KiK-net stations shows a peak at short periods, likely due to the shallow-soil induced resonance, this peak disappears because of the migration to site category E. As a consequence, the median trend of site category A (2021-draft) for KiK-net stations becomes similar to that observed for ESM. Figures 8 and 9 compare the exponential values of the averaged δS2S A (e δS2S ) with the 2021-draft amplification functions, for the ESM and KiK-net datasets, respectively. The empirical amplifications are computed averaging δS2S A in the interval T = 0.07-0.4 s for short periods and in the interval T = 0.7-2.0 s for intermediate periods. The short (F α ) and intermediate (F β ) periods amplifications are computed for a representative low seismicity level, namely, F α is evaluated for S α,RP = 0.2 g and F β for S β,RP = 0.05 g.
The amplification functions introduced by the 2021-draft are in reasonably good agreement with the results obtained by the ESM dataset, both at short periods and intermediate periods. A similar agreement is also observed for Japanese data, and such finding is strengthened thanks to both the larger amount of experimental data and to the fact that Japanese records are absent from the calibration dataset of the reference GMM of KOT20.
As far as variability is concerned, it is generally large for both datasets, especially at short periods. Large variability can also be observed for class F at intermediate periods, probably depending on the reduced number of available data.
As regards class E (Fig. 10), we considered different groups of stations according to three ranges of H 800 , namely, H 800 = 0-10 m, 10-20 m and 20-30 m. The comparison is shown only for the KiK-net subset of class E sites, because it consists of a relevant number of 226 stations. While there is a reasonably good agreement of the F α factor with δS2S A in the 20-30 m thickness range, in the other ranges, and especially at the intermediate periods, a general trend is found where the F α and F β factors overestimate the average δS2S A values. Note however that at intermediate periods and low thickness the average δS2S A is lower than 1, showing that these KiK-net sites tend to de-amplify also with respect to the reference V s,30 = 800 m/s of the KOT20 GMM.

Check of the 2021-draft site amplification factors against 1D stochastic numerical simulations
A large set of stochastic 1D numerical simulations of seismic ground response under vertical wave propagation was constructed, based on 3,202,500 equivalent linear viscous-elastic analyses (EQL, Idriss and Seed 1968), carried out with the SHAKE91 code (Sun and Idriss 1992). The dataset provides a large numerical benchmark for comparison with the 2021draft SAFs, including an assessment of the statistical dispersion within each site category.

Generation of subsoil models
The 1D ground models were generated through a Monte-Carlo procedure, starting from 252 real sites covering the V S,H -H 800 domain shown in Fig. 11a. Most data belong to publicly available databases, such as ITACA  . The collection also includes results from specific site characterization studies (Comina et al. 2011;Minarelli et al. 2016;Capilleri et al. 2009). The population of 1D ground models was then generated using the geostatistical model by Passeri et al. (2020). Two criteria were used to select the 1D ground models from the database, according to the site categories defined in the 2021-draft. On one side, the same number of models was assigned to each site category, to guarantee uniform inter-category representativeness. On the other side, the selection sought to obtain a homogenous intra-category distribution of models, to give equal weight to each specific site condition. For these purposes, every site category -represented by a region in the V S,H -H 800 domain -was discretized into 100 blocks with equal size, as shown in Fig. 11d. The region of interest was limited down to 200 m depth. Each block contained a maximum of 200 models. The limit was set to 20 in the upper triangular portion of the region corresponding to the site category B. This area, indeed, clusters thick deposits with stiff layers in the superficial portion (e.g. deep stratifications of altered rock), that are rather rare in common practice.
In the generation of V S profiles, each real soil model was taken as base-case soil profile, coincident with a median V S profile, from which layers' thicknesses and S-wave velocities were generated with a suitable number of realizations (Fig. 11b). The generated V S profiles were able to cover a region of the reference V S,H -H 800 domain (Fig. 11c). By merging the models generated from each real profile, the V S,H -H 800 domain under examination was gradually filled (Fig. 11d). The result is a collection of 91,500 1D ground models, which exhibit realistic features and represent with equal probability different subsoil conditions of engineering interest. For consistency between the ground response database and the 2021draft site categories, the bedrock V S was limited within the range between 800 m/s and 900 m/s, therefore closely approaching the reference condition that the ESHM20 European seismic hazard model is based on, as discussed in the previous section.
A preliminary analysis of the results led to remove soil profiles implying shear strains larger than 1% for at least one analysis, considered to be a threshold above which nonlinear phenomena become too large for the EQL model to provide reliable results (e.g., Matasovic and Hashash 2012; Kaklamanos et al. 2013), as well as those where convergence was  (Papaspiliou et al. 2012). The operation led to a final population for each site category of at least 3,000 profiles.
EQL analyses require a description of the soil nonlinear behaviour, in terms of modulusreduction and damping (MRD) curves. For this purpose, this study assumed four possible soil materials, i.e. clay, sand, gravel and weathered rock. The MRD curves were assigned based on the following models: Darendeli (2001) for sandy and clayey materials; Rollins et al. (1998) for gravels; Sun and Idriss (1992) for rock-like materials. Further details about the assignment of the material type and the definition of the remaining soil parameters for the ground response analysis -i.e., unit weight, plasticity index, over-consolidation ratio, "at-rest" lateral pressure coefficient, groundwater depth -are available in Aimar et al. (2020).

Seismic input
The input motions consist of a collection of acceleration time histories, selected to represent different levels of seismic hazard, in terms of S α,475 and S β,475 . For this purpose, five reference sites in Italy were considered (Table 7). At each site a suite of 7 accelerograms recorded in flat rock-like outcropping conditions (i.e., on category A sites) was selected, based on the seismic hazard disaggregation in Italy (Barani et al. 2009) and on the spectrum compatibility criteria implemented in the software Inspector (Acunzo et al. 2014). In the EQL analyses, the input motions were applied at the rock outcrop. Figure 12 shows the average response spectra of the input motions for each reference site. Further details about the selected acceleration time histories are available in Aimar et al. (2020).

Fig. 12
Average response spectra of the input motions for each reference site (see Table 7)

Intra-category dispersion of results
An effective site categorization system should be able to cluster various soil conditions into homogeneous groups, i.e. classes of elements sharing similar response with limited dispersion. Therefore, the intra-category variability of results in the site categorization system introduced by the 2021-draft is compared with the one resulting from the EC8-1, to verify whether the new scheme leads to an improvement in this respect.
For this purpose, the numerical amplification factors were firstly computed assuming the same period intervals (i.e. T = 0.07-0.4 s for F α and T = 0.7-2.0 s for F β ) adopted for the empirical amplification factors introduced in the previous section. The results were then clustered within each site category and the corresponding standard deviation was computed. The estimate assumed a lognormal distribution of the results within each site category (Li and Assimaki 2010;Aimar et al. 2020). Moreover, the computation accounted for the dependence of the SAFs on V S,H (and H 800 , in site category E), by estimating the logarithmic standard deviation across narrow ranges of V S,H (and H 800 ) and averaging it over each site category. In this way, the effect of the novel formulation for the SAFs on data dispersion was also investigated.
As shown in Fig. 13, the logarithmic standard deviation ranges from around 0.35 for F α down to 0.2 for F β . At long vibration periods, in fact, simulations provide lower scattering because long-period waves sample a large portion of the soil profile and their propagation mostly depends on global features of the deposit, rather than local variations (Stewart et al. 2014). The 2021-draft scheme significantly reduces the variability with respect to EC8-1, for multifold reasons. On one side, the new parameter V S,H overcomes the limitations of the V S,30 and it allows more effective separation of models with different stratigraphy-and different response-reducing the intra-category dispersion. Furthermore, the new site category F takes a part of the uncertainty from deformable models in site categories C and D of EC8-1, as it isolates very deep deformable soil models from shallower ones. Finally, the continuous formulation of the SAFs helps in reducing data variability because it accounts for the dependency of stratigraphic amplification on soil parameters.

Comparison with 2021-draft amplification factors
The distribution of results of 1D numerical simulations is compared to the SAFs of 2021draft. This check focuses both on the proposed continuous formulation of the SAFs as a function of the site parameters (i.e., V S,H , H 800 , S α,RP and S β,RP ) and on the default values. Figure 14 compares the distribution of results, represented by the interval defined by one logarithmic standard deviation, with the distribution of F α and F β values obtained for all possible V S,H and H 800 combinations, in each site category. In addition, a symbol is added to represent the default values of SAFs of 2021-draft. For almost all cases, the 2021-draft SAFs are in good agreement with the distribution of results at short vibration periods. For stiffer soil models, the estimated F α is close to or slightly below the data mean for low seismicity and it increases for high-seismicity areas. This effect is relevant in site categories C and F, where the predicted F α according to the continuous formulation gradually shifts from the lower part to the upper part of the distribution of results. As for site category D, instead, the 2021-draft gives an estimate on the safe side for almost all cases, as an effect of the prudential way with which nonlinearity is accounted for in the 2021draft. Conversely, in site category E, the SAFs estimated through the continuous formulation mostly lie in the lower part of the distribution, whereas the default value is larger than the data mean, independently of the intensity of the ground motion.
At intermediate vibration periods, the continuous formulation provides a F β estimate almost coincident with the data distribution in site categories B and C, whereas the value computed according to the default SAFs is on the safe side. In soft soil models, both the estimates tend to be on the safe side with respect to the results. A similar finding is observed also in site category E, even though the predicted value based on the continuous expression spans over a wide range of values, partly superimposed with the data distribution and partly exceeding it.
Focusing on the continuous formulation for the SAFs, Figs. 15 to 17 present a more detailed comparison with the distribution of the simulation results, as a function of V S,H and H 800 . These results refer to Godrano site (S α,475 = 2.79 m/s 2 ), but findings are similar for the other sites.
The distribution of F α obtained from the analyses (Fig. 15) exhibits a gradual variation with respect to V S,H . This is significant in the site category B, where F α decreases for growing V S,H down to unity, due to the strong dependence of the amplification on the impedance ratio. The 2021-draft formulation accounts for this behaviour, thus ensuring agreement with the simulations. On the contrary, the simulation-based F α decreases at smaller V S,H in deformable models due to the strong attenuation induced by large strain levels. A significant increase in data variability is also observed, probably due to nonlinear behaviour and to the increased relevance of additional parameters on the response, such as the plasticity index (Ciancimino et al. 2018;Aimar et al. 2020). The effect of nonlinearity is accounted for in the 2021-draft, but the decrease of the predicted F α is smaller compared to the numerical simulations. Therefore, the 2021-draft provides an estimate generally on the safe side to cope with the large uncertainties affecting the ground response in soft soils. Additionally, EQL analyses may provide excessive de-amplification for deep and deformable soil profiles, due to MRD curves poorly constrained at large depths.
As for site category E, some discrepancies are observed between simulated and predicted F α values, as illustrated in Fig. 16. In fact, simulations predict a relatively sharp peak, observed in soil models whose fundamental period -computed according to the formula T 0 = 4H 800 /V S,H -ranges between 0.15 s and 0.4 s. Since the range is included in the reference one adopted to define F α (i.e., 0.07-0.4 s), the resonance of the soil models results in a significant increase of the amplification. Nevertheless, the fitting of the 2021-draft formulation is satisfactory, with an estimation on the safe side for most soil profiles, similarly to what observed in the previous section comparing with the KiK-net records.  (Fig. 17). F β slightly increases with decreasing V S,H . In stiff soil deposits, i.e. site category B, this growth is mainly due to the impedance ratio, similarly to what observed for F α . As for soft models, a reason of the increase of F β at small V S,H could be the strong nonlinearity in the simulated ground response. The soil column, in fact, undergoes a stiffness decay, entailing an increase of the fundamental period of the ground model, which in turn leads to amplify the long-period components. In site category E, both the 2021-draft estimate and the median of the results are characterized by larger F β values for decreasing V S,H and increasing H 800 , corresponding to an elongation of the fundamental period.

Nonlinear effects on site amplification factors
It is well known that NL soil response under high intensity seismic excitation typically shows itself by a reduction of the peak amplification frequency, as a consequence of shear modulus reduction, and by a reduction of the high frequency spectral components, as a consequence of increased damping levels. As shown in the previous section, discrepancies of the 2021-draft and numerical simulated SAFs mostly occur in case of soft/deep soil deposits because of NL effects that, in the numerical simulations, tend to provide smaller amplification levels (or de-amplification) than in the 2021-draft.
Quantification of NL effects on the SAFs is still a controversial issue. As a matter of fact, difficulties to properly identify correction terms on SAFs to account for NL effects mostly stem from the limited amount of worldwide stations presenting both strong and weak motion records with the simultaneous presence of a close-by reference site. Because of that, empirical approaches attempt to check a possible intensity dependence of SAFs by statistical evaluations of ground motion amplitude at generic soil sites with respect to the reference rock conditions (e.g., Perus and Fajfar 2014). In this framework, the NEHRP (2015) factors were supported by the semi-empirical approach introduced by Stewart and Seyhan (2013)  other models such as Lanzano et al. (2018) and Kotha et al. (2020) neglected such intensity dependence. As shown in Fig. 4, the nonlinear parametrization introduced in Table 4 by the coefficients r α and r β , was calibrated in order to approach, albeit within a simplified form that is suitable to a seismic norm, the Akkar et al. (2014) amplification curves at least in the intensity range that is mostly interesting for applications in Europe. It is worth to point out that the nonlinear parametrization in Table 4 is similar to that proposed by Dobry et al. (2000), who used the functional form of the type (V s /V ref ) m , m being dependent on PGA R at the outcropping reference rock. Figure 18 shows the comparison of the formulas in Table 4 for ground categories B to D against the closed-form expressions provided in NEHRP (2015) for the short period amplification factor F a (equivalent to F α ) and for the intermediate period factor F v (equivalent to F β ). A reasonable agreement is obtained, especially for intermediate values of the reference spectral ordinates. However, the spread of SAFs from NEHRP (2015) is larger: more specifically, the SAFs for low seismicity are higher, while those for high seismicity are lower. Indeed, NEHRP (2015) assumes a strong influence of site amplification effects when the intensity of ground motion increases. Considering for example V s,30 = 200 m/s, the NEHRP (2015) F α values decrease by a factor of about 2.5 when increasing S α,RP from 0.2 g to 1.2 g, against a factor of about 1.5 for the 2021-draft. In terms of F β , the decrease is by a factor of about 2 according to NEHRP and of about 1.3 according to 2021-draft, when increasing S β,RP from 0.05 g to 0.5 g.
In alternative to empirical approaches, NL effects on the SAFs are often estimated by 1D NL numerical simulations. Although widely used both in research and in practical applications, as shown in the previous section, the validation of such approaches is still a major issue, which has not, up to now, received satisfactory scientific answers (Pecker et al. 2017). Such validation was addressed in the recent years in the blind prediction test of Turkey Flat (Kwok et al. 2008) and in the Prenolin project (Régnier et al. 2018), the latter having involved 19 worldwide teams with 23 different codes. In the Turkey Flat test, a general bias of simulations towards under-prediction at high frequencies was found and attributed to inaccuracies in prior estimates of near-surface shear-wave velocities. In the Prenolin project, two benchmark sites were available, with large amplitudes of recorded ground motions and an extensive amount of in-situ and laboratory soil properties. A substantial code-to-code variability was found (in the 0.05-0.25 range in log10 scale), that may be considered as a lower-bound of the overall epistemic uncertainty, because all simulations were carried out under the simplifying assumption of vertically propagating polarized SH waves in 1D layered soil. Furthermore, the misfit with records was even more substantial, with increasing values (up to 0.25-0.35 in log10) for increased level of shaking, especially for the deeper soil deposit (depth to the engineering bedrock of about 40 m). It was concluded that, while the prediction of NL site response is relatively easier for shallow soil deposits, provided constitutive models are carefully calibrated and compared with laboratory test in a wide strain range, application of 1D NL approaches to deep deposits not only requires sample measurements at varying depths, a task beyond reach in the vast majority of applications, but the 1D assumption itself may not be valid.
Based on the findings from these validation experiments one can justify the "prudent" attitude in the 2021-draft in the quantification of NL effects: while a moderate reduction of SAFs is assumed by the r α and r β factors in Table 4, introduction of larger reductions may only be justified by a detailed investigation of the local mechanical soil properties that falls largely beyond the accuracy of the soil categorization scheme adopted in a seismic norm and that can be worked out within site-specific response analyses.
As stated previously, the lack of records at variable level of intensity simultaneously obtained on both soft soils and close-by outcropping rock stations prevents a comprehensive validation of NL site amplification factors. However, the Japanese KiK-net provided a substantial amount of useful data for this purpose, being equipped with recording stations installed both at the ground surface and at the borehole level (S and B stations in Fig. 19, respectively), the latter one typically being placed at depth variable from few tens to few hundreds of m. From the statistical analysis of the KiK-net database, Régnier et al. (2013) were able to identify conditions for onset of significant NL effects in the seismic site response. However, only few of the investigated stations were rich in both high-and low-intensity motions to allow for a quantification of sitespecific NL effects in a broad intensity range that could be suitable within a seismic norm. Besides, the amplification levels S/B in the frequency domain are an upper bound of the corresponding levels S/R, because of the destructive interference at location B from reflected motions from the ground surface and from the intermediate layer interfaces. Therefore, the response spectral ratio S/R and its intensity-dependence cannot be directly deduced from S/B.
To provide an indirect estimate of the S/R response spectral ratio and of its intensity dependence, the following procedure was devised: • several KiK-net stations were selected, among the very few ones rich in both high and low-intensity motions; • the soil profile at each station was iteratively modified, starting from the default one provided in the KiK-net database, in order to reproduce, as closely as possible, the observed Fourier spectral ratios S/B and their intensity-dependence: therefore, for each intensity group and each station, a different soil profile was obtained, in terms of V s and soil damping ratio D; • the ground motion at R station was obtained, by linear 1D deconvolution, for each of the available records at S; • the response spectral ratio S/R was computed and plotted as a function of the ground motion amplitude at R.
Referring to Mattevi (2018) for a more detailed introduction of the approach and a broader set of stations, we consider here the most meaningful ones, both in terms of significance of recorded NL effects and of variety of intensity of records, given by the important seismic sequences in Japan of Kumamoto (2016) and Niigata (2004-2007. In Fig. 20, the magnitude-distance distributions for the three selected stations are illustrated. Each station contains at least 10 records with PGA (S station) > 100 gal, that is an almost unique amount of high intensity records.
In Fig. 21 (top row), the S/B average Fourier spectral ratios at each station are shown, grouped according to ranges of PGA (S station). It is clear that, as PGA increases, the ratios at KMMH14 and KMMH16 tend to shrink showing an increasing depletion of high frequencies, while at NIGH11 they remain almost invariant. Based on such experimental results, an iterative set of 1D linear visco-elastic numerical simulations was carried out at each site in order to obtain the best agreement of the numerical with the experimental S/B ratios. The resulting V S profiles are shown in Fig. 21 (medium row), referring the profile S1 to the group of records with the lowest intensity and the profile S4 (for KMMH14) or S5 (for KMMH16) referring to the group with the highest intensity; as expected, V s,30 decreases as PGA increases. A frequency-dependent soil damping inspired to Kausel and Assimaki (2002) was also tuned, to improve the fit with experimental S/B ratios (further details in Mattevi 2018). The bottom row of Fig. 21 shows the comparison of the experimental and numerical S/B ratios for the high intensity group, together with the corresponding S/R ratio. Based on such calibration, it was possible to estimate the ground motion at the ideal outcropping rock station R by deconvolution of recorded ground motion at station S for each recorded earthquake. Figure 22 shows the S/R response spectral ratios at T = 0.2 s and T = 1 s, as a function of PGA R , normalized to the average low-intensity value, in order to highlight the potential variability of such ratio with increasing values of PGA R at the outcropping station R. In the same figure, the red curve shows the normalized decay of the site amplification factors obtained by the r α and r β factors of (1 s) = 0.7 ⋅PGA R were used to refer plots to the same horizontal axis. Furthermore, the corresponding decay curve (black) from the NEHRP is also plotted. V s,30 values are based on the V s profile calibrated to fit the low amplitude S/B spectral ratios (see profiles S1 in Fig. 21).
Nonlinear effects are clearly visible at KMMH16 (V s,30 = 268 m/s), with evidence of decay of the SAFs starting from about PGA R = 0.1 g but limited to short periods, and at KMMH14 (V s,30 = 190 m/s) where the decay is sharper, ranging at short periods (T = 0.2 s) from about 0.9 at around PGA R = 0.05 g down to about 0.6 at around PGA R = 0.25 g. As expected, the decay at intermediate periods (T = 1 s) is smoother, down to about 0.75 at 0.25 g. At NIGH11 (V s,30 = 375 m/s), there is no evidence of significant nonlinear response, as clearly pointed out by the experimental S/B Fourier spectral ratios that remain almost invariable with increasing PGA. The comparison with both the 2021-draft and NEHRP reduction factors shows a reasonably good agreement. NEHRP formulas provide a slightly better approximation at KMMH14, where the V s,30 value is lower, while the 2021-draft ones fit better KMMH16 and NIGH11, for which NEHRP overestimates the decay, especially at short periods. This evidence provides a justification for a conservative approach, as implemented in the 2021-draft.
As recalled previously, very few worldwide recording stations are available to extend this check to further case studies. In Mattevi (2018), results for KiK-net stations KMMH03 and IWTH23 are also shown, not reported here for brevity, confirming the overall good performance of the 2021-draft decaying factors.

Conclusions
In the 2021-draft of Part 1-1 of Eurocode 8 (CEN/TC250/SC8 2021), the site categorization criteria and site amplification factors were modified with respect to EC8-1 (CEN 2004), with the main objective to reduce classification ambiguities (as shown in Fig. 11a) and to provide a continuous formulation of SAFs that could lead to a smooth transition of amplification levels from class to class. This was obtained by the explicit introduction of two proxies (V S,H and H 800 ) to describe site categories, supported by simplified criteria for categorization, presented in the Annex of the 2021-draft, in case of incomplete information.

3
The continuous formulation of short and intermediate period SAFs in the 2021-draft (Table 4) was introduced to approach the median period-dependent values obtained by earthquake records and expressed through the site-dependent coefficients of GMMs. There are several reasons to justify record-based SAFs for seismic norms. First, GMMs are one of the key elements for seismic hazard assessment. Since this is typically evaluated on reference ground through uniform hazard spectra (UHS), the site amplification factors available from GMMs provide the most coherent way to relate the expected ground response at a reference site with that at a generic site. Furthermore, it can be shown (see e.g. Iervolino 2016) that application of GMM-based SAFs provides hazard-consistent spectra at generic sites, at least when NL effects on ground response are limited. This means that the resulting spectra after application of SAFs have the same probability of exceedance as the UHS at reference ground, that is a major advantage from the point of view of a normative application. Second, record-based SAFs implicitly include complex propagation and site amplification effects that may be hard to be quantified in a suitable and easy-to-use format for seismic norms and that 1D numerical simulations do not account for (see e.g. Pilz and Cotton 2019). This is the main reason why the GMM-based SAFs provide generally larger values than predicted by 1D simulations at intermediate and long periods (see e.g. Figure 17c).
However, while the available GMMs can only capture the trend of amplification at generic sites within broad categories, 1D numerical simulations have the role to quantify patterns of amplification in those specific soil conditions not well covered by records, such as velocity inversions, shallow soil profiles with different possible interactions with the frequency content of the input, large soil non-linearities.
Given these premises, we performed a comprehensive check of the SAFs presented in the 2021-draft involving: a. The statistical evaluation of residuals of records from the European (ESM) and Japanese (KiK-net) strong motion datasets with respect to the KOT20 ground motion model, that is in the background of the ESHM20 seismic hazard map produced in support of the revision of the Eurocode 8 , therefore allowing a consistent check in view of the future applications; b. The analysis of migration of available strong motion stations from the EC8-1 to the 2021-draft site categorization, showing some significant changes, involving especially the updated definition of site category E for shallow sites and the introduction of the category F, specific for deep soil sites; particularly, it was shown that the 2021-draft scheme significantly reduces the data dispersion for all site categories with respect to EC8-1, as a natural outcome of the reduced ambiguity of the site categorization criteria; c. The analysis of SAFs from a comprehensive set of 1D numerical simulations involving a stochastically generated set of V S profiles and different intensities of input motions consistent with Italian sites with low-to-high levels of seismicity; d. The check of the r α and r β factors, introduced in Table 4 to account empirically of NL effects on site amplification, by the detailed analysis of amplification at some of the few KiK-net stations with a significant availability of both weak and strong surface and borehole records.
The overall results of this checking activity were satisfactory, showing that the continuous variability of the 2021-draft SAFs with V S,H reproduces reasonably well both the average δS2S A values computed from the ESM and KiK-net datasets with respect to the KOT20 GMM, computed for the reference condition V S,30 = 800 m/s. Agreement with the KiK-net values is noteworthy, because this dataset was not considered in the development of the European dataset and corresponding ground motion models.
The main discrepancies are found for soil category E, for which the SAFs of the 2021draft tend to lie above the average δS2S A values computed only from the KiK-net dataset, because site class E is not sufficiently well represented in the ESM dataset. Such discrepancies, confirmed by the 1D numerical simulations, suggest that it is precisely in such site category where the largest variability of site amplification effects occur, owing to the complex interaction of wave velocity, of soil deposit thickness and of frequency range of the excitation. For this reason, it is difficult to cast such variability in a simplified form suitable for introduction in a norm.
Simulated results are in good agreement with the 2021-draft values for stiff soils. Conversely, for increasing soil deformability, a stronger decay in amplification is found by the numerical simulations, particularly relevant for high seismicity. Although the 2021-draft values incorporate a reduction in the amplification to account for nonlinearity, the decay is not as strong as suggested by numerical results, that typically predict lower amplification at intermediate periods for deep soil deposits and at short periods for shallow soil deposits. On one side, this mismatch may be due to the limited amount of records potentially affected by strong NL response, but, on the other side, numerical results may be biased by limited in-field experimental data to reliably constrain NL parameters for the generic site categories that seismic norms are based on. To throw light on the available in-field data, a site-specific study on NL effects at three KiK-net stations with a significantly large amount of high-intensity recorded ground motions was performed, providing some support to the 2021-draft NL reduction factors, although the very limited number of recording stations allowing such analysis prevents deriving more general conclusions.
To conclude, there are several reasons for supporting the "prudent" attitude of the 2021draft towards the reduction of SAFs owing to NL effects, namely: (i) a proper evaluation of NL response of soils depends on a level of information (e.g., soil strength, MRD curves) that cannot be introduced in the site categorization criteria of a standard; (ii) larger reduction factors may be possible, but they should be supported by site-specific investigations that, as clearly shown by the results of international benchmarks, are crucial to properly constrain the simulation results; (iii) although experimental evidence of NL effects on strong motion records is still quite limited, the 2021-draft NL reduction factors seem to provide a reasonable agreement with records, as shown in Fig. 4 and Fig. 22. numerical study was performed by Mauro Aimar, Andrea Ciancimino and Sebastiano Foti. Each team wrote the corresponding sections of the manuscript and the parts were merged by Roberto Paolucci, who also wrote the descriptive section about the building code. Then, all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Funding Open access funding provided by Politecnico di Torino within the CRUI-CARE Agreement. Partial funding for participants from Politecnico di Milano and Istituto Nazionale di Geofisica e Vulcanologia (RP, GL, FP) was provided by ReLUIS Project WP18 "Research contributions to definition of seismic actions for Design" (2019-21), promoted by the Italian Civil Protection Agency. Partial funding for participants from Politecnico di Torino (MA, AC, SF) was provided by the ReLUIS Project WP16.1 "Site response analysis and liquefaction" (2019-21), promoted by the Italian Civil Protection Agency.
Code availability In the numerical study, ground response simulations were carried out with the SHAKE91 software (available at https:// nisee. berke ley. edu/ elibr ary/ getpkg? id= SHAKE 91) and the Inspector software (available upon request at inspector.igag@gmail.com) was used to help in the selection of the appropriate ground motions. Data processing and figures were done using MATLAB (http:// www. mathw orks. com/ produ cts/ matlab/).

Conflicts of interest
The authors have no conflicts of interest to declare that are relevant to the content of this article. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the funding agencies.
Ethics approval Not applicable.

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