Physics-based parametrization of a FAS nonergodic ground motion model for Central Italy

We propose a new fully nonergodic ground motion model for Central Italy, which is one of the most sampled areas in the world after the occurrence of the last seismic sequences of 2009 and 2016–2017. The model predicts 69 ordinates of the Fourier Amplitude Spectrum in the magnitude range 3.2–6.5 and is constrained on a dense set of seismological and geophysical parameters (i.e. stress-drop Δσ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \sigma$$\end{document}, shear-wave velocity in the uppermost 30 m, VS,30 and high-frequency attenuation parameter at source κsource\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\kappa }_{source}$$\end{document} and site κ0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\kappa }_{0}$$\end{document}) made available from a non-parametric generalized inversion technique (GIT). The aim of this work is to capture the underlying physics of ground motion related to different source energy levels, as well as to the crustal and geological structure of the region, thus providing less uncertain predictions. Calibration is performed using a stepwise regression approach which has the advantage of taking a more complex functional form (advanced model) when all physical parameters are known while returning a simpler form (base model) when physical data are missing. As a result, the advanced model reproduces the reference rock motion of the region in case the site additional proxies are set to their average values (VS,30 = 1100 m/s, κ0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\kappa }_{0}$$\end{document}=15 ms). We show that the inclusion of the set of physically-based explanatory variables in the regression has a beneficial effect in constraining the uncertainty, leading to a reduction of the high-frequency variability of about 70% on the between-event and 35% on the site-to-site. This reduction can be viewed as the result of the combination of a more effective physical description through the incorporation of the additional proxies and a calibration embedded in a completely nonergodic framework.


Introduction
Ground motion models (GMM) represent key-elements in applied seismology as they mathematically describe the effect of earthquakes at the site, given a scenario event defined in terms of at least magnitude and distance. They are also basic ingredients in probabilistic and deterministic seismic hazard analyses, also used for loss and risk assessment, as well as in many tools for civil protection planning and engineering applications, like shaking scenarios and Shakemaps (Worden et al. 2018). The main drawback of these models rely in the fact that they represent extreme simplifications of the earthquake process and also because they are typically calibrated on limited data that are available only in high-seismicity areas of the world and then applied as global averages without considering the specific location, thus invoking common ergodic assumption (Anderson and Brune 1999). The use of the ergodic hypothesis typically leads to large uncertainty associated with the ground motion predictions (Stafford 2014), which strongly impacts on the results of hazard analysis and related products, as highlighted in the literature .
For these reasons, in the last years, a strong effort has been made by the research community to reduce ground motion variability and move towards regionalised models at smaller scale. The most acknowledged and advanced strategy to this aim is based on the relaxation of the ergodic assumption; a shift that is nowadays possible thanks to the increasing availability of records consequent to the growing of seismic networks worldwide. A fully nonergodic approach is based on the principle that wherever many records at each station from multiple earthquakes are available, it is possible to identify the systematic contributions of variability into event-, source-, site-and path-effects, through a statistical decomposition technique (Al-Atik et al. 2010). It is estimated that about 60-70% of the aleatory variance in ergodic GMM is due to these systematic effects .
Examples of recent fully nonergodic models, conducted either with continuous or spatial regionalization approaches, can be found in Landwehr et al. (2016), Baltay et al. (2017), Parker et al. (2022), Kuehn et al. (2019), Abrahamson et al. (2019), Sgobba et al. (2019), Lanzano et al. (2021), Sung et al. (2022), Lavrentiadis et al. (2021;, Liu et al. (2022) among the others. All these models are associated with random-terms that act as adjustments of the median prediction, while moving part of the aleatory variability into epistemic uncertainty (Al-Atik et al. 2010;Anderson and Uchiyama 2011). These studies have been typically developed in highly sampled regions such as California or Italy. In this latter country, the occurrence of the last seismic sequences in 2016-2017 in Central Italy, has enhanced a lot the Italian datasets allowing to gather geophysical, geological and seismological information that have driven to a fruitful activity of empirical modeling in this region. In detail, this exceptional set of information has produced many studies in the area that went into detailed parametrization of the regional geology (Chiaraluce et al. 2017;Buttinelli et al. 2021;Di Bucci et al. 2021) and prediction of ground motion at the site (Bindi et al. 2018a, b;Morasca et al. 2019;Lanzano et al. 2020;Felicetta et al. 2018;Sgobba et al. 2021;Lanzano et al. 2021;Castro et al. 2022;Colavitti et al. 2022).
In this work, we want to further exploit the mine of data and seismological parameters in Central Italy, to complement the nonergodic model calibrated by Sgobba et al. (2021) for spectral accelerations (SA), with a new ad-hoc GMM based on Fourier Amplitude Spectra (FAS) by using a set of fully consistent physical information to predict the ground motion parameters and capture more epistemic uncertainty. In doing this, we consider that the aim of setting seismological constraints to ground motion can be achieved more appropriately in terms of FAS compared to response spectral accelerations SA, because FAS are more closely related to the physics, as also highlighted by Bora et al. (2015) and Bindi et al. (2018c). Hence, our study is conducted in terms of a GMM-FAS similarly to Bora et al. (2015) and Bayless and Abrahamson (2019), who proposed a physically-consistent approach for developing FAS predictions easily adjustable to different seismological conditions. Yet, compared to those works, we here investigate the relationships between the systematic terms of uncertainty of the FAS-GMM and the physical parameters related to the source rupture and site response gathered from independent analyses, under a complete nonergodic framework. An advantage linked to the inclusion of more physical constraints in the predictions consists in the opportunity to assess the impact of a boosted parametrization on the aleatory variability reduction. The final purpose is to define a methodological framework for calibrating a model of FAS that, once converted in spectral intensities, could be used to generate non-ergodic and physics-based empirical shaking scenarios (e.g. Sgobba and Pacor 2023;Sgobba et al. 2021;Abrahamson et al. 2019) for civil protection purposes or that could be also exploited to produce expeditive simulations of groundmotion time series.
The paper is outlined as follows: we first develop a nonergodic FAS-GMM on 69 ordinates of the Fourier spectrum including only source and attenuation scaling terms and related to reference rock sites. Then we correlate the repeatable terms of variability with a suite of region-specific spectral parameters (i.e. stress-drop, shear-wave velocity in the top 30 m and high-frequency attenuation kappa) obtained from a non-parametric generalized inversion technique (GIT) developed on the same dataset by Morasca et al. (2023). Finally, we re-parametrise the original model including the additional explanatory variable through a stepwise regression approach and then compare the prediction performance and uncertainty contributions among the two model versions.

Dataset
The collection of records used in this study is composed of both accelerometric and velocimetric three components waveforms of events located in Central Italy and occurred between 2009 and 2018 including the latest major seismic sequences that occurred in the area of study, i.e. the 2009 L'Aquila seismic sequence and the 2016-2017 Amatrice-Visso-Norcia sequence.
This dataset was originally developed in the framework of the seismic microzonation study carried out in Central Italy after the 2016 M w 6.0 Amatrice earthquake (Priolo et al. 2020), but it was also extensively used to study the temporal and spatial variability of the ground motion in the area (Bindi et al. 2018b;Sgobba et al. 2021).
It consists of more than 30,000 waveforms relative to about 450 earthquakes in the magnitude range 3.2-6.5 (local magnitude for M < 4.5 and moment magnitude for M ≥ 4.5) recorded by about 460 stations within 250 km from the epicenters (Fig. 1a).
The magnitude-distance distribution plot (Fig. 1b) shows that the region is highly sampled in the distance and magnitude ranges [0-100] km and [3.2-4.5], respectively. The large amount of records at short distances and for small events is due to the stations belonging to the temporary networks installed during the seismic sequences. This high number of stations and events in such a relatively small area has allowed us to sample a relevant number of ray-paths (Fig. 1c), which is a requirement to apply a fully nonergodic approach.
On this dataset, we analyze the Fourier amplitude spectra (FAS) of S-wave windows band-pass filtered with a variable high-pass corner frequency depending on the signalto-noise ratio. The time windows used to calculate the spectral amplitudes were selected using a distance dependent energy criterion. The FAS are smoothed using the Konno and Ohmachi (1998) algorithm (the smoothing parameter b was set to 40). Details about the data selection and processing are provided by Pacor et al. (2016), Bindi et al. (2017), and Colavitti et al. (2022). Fig. 1 a Map of the study-area with the spatial distribution of events and stations.; b Scatter plot of the dataset records as a function of magnitude and source-to-site distance; c Map showing the source-to-site paths in the study-area

Base model
As the aim of the work is to assess the median and variability of a GMM-FAS including regional seismological parameters, we perform two calibrations: the first one allows us to develop a nonergodic base model including basic dependencies and scaling terms; then a second calibration is developed to include the additional explanatory variables.
The functional form, adopted frequency-by-frequency for the base model is: where fixed-and random-effects can be identified. The fixed part is composed by the term a, which is the offset of equation and the magnitude F M (M) and distance F R (M,R) scalings that are treated adopting standard dependencies, as follows: In particular, we assume a non-linear magnitude scaling through a stepwise linear function ( F M term) and the distance scaling is composed by a geometrical spreading (including a magnitude dependent term) and the anelastic attenuation terms ( F R term).
The explanatory variables are the moment magnitude M w and the source-to-site distance R. In detail, for the strongest events (M > 5.5), R is the Joyner and Boore distance (R JB ), computed from the fault geometries published in the ITalian ACcelerometric Archive (ITACA) database (Russo et al. 2022), whereas for smaller events the epicentral distance R epi is assumed for R. It is worth to note that we do not include the focal depth among the parameters, because the dataset refers to an homogeneous tectonic regime (extensional) characterized by events with roughly constant depths (< 20 km), normally distributed (mean 9.3 km; standard deviation 3.1 km), as also stated by Sgobba et al. (2021), with reference to the spectral acceleration-based model calibrated on the same dataset. Note also that the functional form lacks a site-response scaling term (typically dependent on the shear wave velocity V S,30 ), so that all the unmodelled site effects turn out in the corresponding random-terms S2S s .
In Eqs. (3.2) and (3.3), some parameters are fixed in a first stage by nonlinear regression, i.e. the hinge magnitude M h = 5.0, the reference distance R ref = 1 km and the pseudodepth h = 6 km. The reference magnitude M ref is instead frequency-dependent estimated from a preliminary non-linear regression; it varies between 5.45 (at frequency f ~ 1 Hz) and 3.3 (at frequency f ~ 7.5 Hz).
The sites adopted for model calibration come from the study of Lanzano et al. (2020), who identified a set of 36 recording stations in Italy installed on rock and classified as "reference sites" that are intended as out-cropping rock or stiff soils that show a flat, unamplified response over a frequency range of engineering interest (Lanzano et al. 2022a). These sites have been identified on the basis of a multi-criteria ranking approach including several proxies that influence the site response, i.e., (i) the outcropping geology, (ii) the installation features, (iii) the shear wave velocity V S,30 according to EC8 (CEN, 2003), (iv) the site topography, (v) the horizontal-to-vertical spectral ratios obtained from noise measurements or recordings and (vi) the repeatable site term obtained from residual analysis (δS2S s ) and that can be considered as the empirical amplification function of each station (Priolo et al. 2020;Lanzano et al. 2020;Paolucci et al. 2021).
In the present study, with respect to the initial 36 stations identified according to the above criteria, a further selection is performed on the basis of the high-frequency attenuation κ 0 parameter extracted from the analysis of the Fourier Spectra (named κ 0 -FAS). The κ 0 -FAS are estimated as the semi-logarithmic frequency decay of the S-waves of FAS (Ktenidou et al. 2014), by using the semi-automatic procedure described in Lanzano et al. (2022a). The final set of 6 reference stations is listed in Table 1 and include five stations with κ 0 -FAS < 0.015 s, namely LSS, MNF, SLO, SNO and SDM, plus a sixth station (NRN), which shows a δS2S s trend very similar to the average of the 5 selected sites, although it does not have an estimate of κ 0 -FAS.
The random terms of Eq. (3.1) include some corrective factors related to the systematic source, propagation, and site effects related to event e, source-region r, site s and path p. The random effects are the repeatable terms of variability that are introduced to remove the ergodic assumption, namely: • between-event B e : represents the average deviation away from the median prediction of the GMM for any event; it shows a clear dependence on earthquake stress-drop (Bindi et al. 2017); • location-to-location L2L r : represents the variability term related to the source region due to different releases of radiated energy; indeed, regions with different average stress drops can be identified on smaller scales (Baltay et al. 2019). In this work, this term is estimated via the clustering approach adopted by Sgobba et al. (2021) with reference to a non-ergodic SA-GMM calibrated on the same dataset adopted herein. The clustering is based on the following criteria: (1) seismogenic setting (considering the source zones identified by the database of seismogenic sources, DISS 3.2.1 catalog -DISS Working Group 2018); (2) similarity in the stress drop values associated with the events of each cluster; (3) space-time aggregation based on Reasenberg's algorithm (Reasenberg 1985): the earthquakes are counted within a time window of 10-days bins and assuming a Poisson distribution. Following these criteria, the events are aggregated within polygonal source areas and identified as follows: cluster #1 (area of L'Aquila, southernmost), #2 (area of Amatrice-Norcia in the middle), and #3 (area of Muccia, farther North). In the rest of the region we considered a background seismicity with a mean value of the location term equal to zero. A map of the clusters is reported in Fig. 2a; • site-to-site S2S s : defines the systematic bias of ground motions recorded at each station s with respect to the reference rock motion predicted by the fixed-effect part of the model on the basis of the 6 reference sites listed above; • path-to-path P2P p : defines the variability from one source-to-site path to another due to differences in travel paths across heterogeneous geological layers or main structural discontinuities. The path terms are estimated by dividing the whole region into squared cells (0.2°-spaced) - Fig. 2b that allow to capture the spatial distribution of the attenuation behavior (cell-specific attenuation). The adopted cell size is a compromise value between the available data and the desired spatial resolution given the size of the region, in line with the approach of Dawood and Rodriguez-Marek (2013), Sgobba et al. (2021) and Sung et al (2022). Moreover, we have verified that the observed variability in the paths ranges from 10 to 100 km, but in most cases is larger than 20 km, which is in line with the adopted grid size.
The regression coefficients a, b 1 , b 2 , c 1 , c 2 and c 3 of Eqs. (3.2) and (3.3) are obtained via linear mixed-effects regressions (Stafford 2014;Bates et al. 2015) and are available frequency-by-frequency together with the values of M ref (see "Data availability" section). All the terms in the functional show a p-value (test of the statistical significance of the explanatory variables, Wasserstein and Lazar 2016) lower than 0.05, indicating a significant correlation between the predictor and the response variable.

Additional explanatory variables
The base FAS-model defined before is described by a simple equation dependent on basic explanatory variables mainly including magnitude and source-to-site distance scaling. This means that all the unmodelled effects related to source, site and path are not included in the median model but are centered in the random-terms, which act as adjustments of the FAS predictions. In order to parametrise these terms and move further epistemic uncertainty from the aleatory variability, we adopt the strategy to add constraints by using the physical parameters made available by Morasca et al (2023) from a non-parametric Generalized Inversion Technique (GIT, Andrews 1986;Castro et al. 1990;Oth et al. 2013;Bindi and Kotha 2020). The latter is a data-driven approach useful to resolve source, path and site contributions from the observed FAS that hence has the advantage to derive region-specific spectral parameters for the analyzed events, i.e. seismic moment, corner frequency, stressdrop, source and site kappa values and quality factor.
Starting from the same dataset of almost 460 events used in this study, Morasca et al (2023) selected the FAS corresponding to 283 stations and simultaneously inverted the FAS to isolate the three different factors (source, path and site) by solving an overdetermined linear system in the least-squares sense: where f is the frequency (analyzed in a range 0.5-25 Hz), R es represents the hypocentral distance between the event e and the station s and M is the magnitude associated to the event e.
In the non-parametric approach applied by Morasca et al (2023), no a priori models are assumed to isolate the different terms. This implies the necessity of a second step to extract the source and path parameters from the non-parametric solutions by fitting them a posteriori. To estimate the source parameters for each earthquake, the authors assumed the Brune (1970) model, while non-parametric attenuation functions were fitted assuming a standard model including a distance-dependent bi-linear geometrical spreading and a frequency-dependent attenuation term accounting for source-to-site distance, propagation medium properties, high-frequency attenuation described by the kappa (κ 0 -GIT) parameter, and the regional Quality factor (Q r ). All the physical parameters obtained by Morasca et al. (2023) are available at the link https:// shake. mi. ingv. it/ centr al-italy/ where the moment magnitude, corner frequency, stress-drop, kappa source and seismic moment are available for each analyzed event together with the amplification functions of each station.
In addition to the parameters derived from the GIT of Morasca et al (2023) on Central Italy -hereafter called M-GIT -we also consider the most common site proxy, that is the shear wave velocity in the uppermost 30 m, V S,30 , as discussed in the next section.
In the following, we illustrate how the random effects of the empirical FAS model are investigated and parameterized: we first evaluate the statistical correlation of the randomeffects on the site and seismological parameters get from GIT and, where found relevant, we incorporate them in the regression as additional explanatory variables.
Details on the application of the M-GIT are described in the Morasca et al (2023); we just highlight here that the inversion is performed in such a way that its results are consistent with those from the FAS base model, due to: (i) the correspondence of the original datasets used in both the analyses, though only 283 stations were considered for GIT, each sampling at least 10 events, to obtain more robust inversions; (ii) the adoption of a common (3.4) log 10 FAS es f , R es , M e = log 10 Source e f , M e + log 10 Path es f , R es + log 10 Site s (f ) reference motion: this is an a priori constraint required by GIT to solve the system of linear equations. Namely, the same 6 reference stations introduced in Sect. 3 were adopted both for the calibration of the base empirical FAS model and the GIT analyses. Owing to this consistency, we assume that the random-terms from the empirical FAS model and the outcomes of the M-GIT are inherently related, so we attempt to capture the underlying seismological connections due to the source and geophysical properties in order to parametrise the model variability and adjust the FAS predictions.

Source parameters
Among the random-terms of the empirical FAS model, the between-event δB e and the location-to-location δL2L r are the components related to the source physics due to energy radiation and tectonic setting. Indeed, it is well established that ground motion at highfrequency depends on earthquake stress drop Δσ (e.g., Boore 1983;Bindi et al. 2017;;Baltay et al. 2019), which reflects in the between-event variability. Also the determination of the location-to-location terms δL2L r associated to event clusters  allows us to capture the variations in the Δσ within the Central Italy region, i.e. a given cluster produces earthquakes that are characterized, on average, by higher or lower stress drops. Therefore, also the sum of the event and location terms (i.e., the total source-related random components δL2L r + δB e ) can be considered related to the stress-drop parameter, which is thus assumed a candidate predictor variable.
Both stress-drop and other source parameters are here extracted by model fitting of the non-parametric source spectra based on M-GIT (an ω 2 -model Brune's model is assumed, Brune 1970). Namely, the non-parametric analysis allows the extraction of event-specific values of Δσ and the high-frequency source parameter source , whose distributions are shown in Fig. 3a and b, respectively. This parameter source is defined as the average slope of the source spectra that accounts for the high-frequency attenuation effects close to the fault (e.g. Oth et al. 2011;Bindi et al. 2018c;Bindi et al. 2017) although its correlation with event and location variability is less studied in the literature. The median values of distribution are approximately equal to 0.15 log10 for Δσ distribution (1.5 MPa) and to 20 ms for the source . Both Δσ and source values are reported in the "Data availability" section. We thus compare the inversion-based estimates of Δσ and source with the total sourcerelated random components at each frequency (Fig. 4). Figure 4 shows that for low frequencies (e.g. f = 1 Hz, Fig. 4a), the trend is weak with values of B e + L2L r close to zero throughout the stress drop range (from 0.1 up to 100 MPa), denoting poor or absent correlation. At high-frequency (e.g. f = 24 Hz, Fig. 4b), the data show a trend that increases as the stress-drop values increase up to about 2 MPa and then it flattens out.
An opposite behavior can be noted for B e + L2L r against the source parameter at highfrequency (negative correlation), Fig. 4d, whereas at low-frequency (Fig. 4c) the trend is weaker.
In light of above, it can be stated that stress drop Δ and source values are anti-correlated as also illustrated, in Fig. 5 where we show the frequency-dependent Pearson's correlation coefficient computed between the terms B e + L2L r versus stress drop Δ and versus source .
Regarding the stress drop, the correlation coefficient is negative for frequencies lower than 1 Hz, then rising to positive values around 0.3 for higher frequencies. Concerning the Summarizing, it appears that the correlation of the source-based random terms B e + L2L r with both the parameters is stronger at high-frequency (above about 5 Hz), as expected from physics-based considerations, albeit opposite for Δ and source , and thus we consider these two parameters significant to further parametrise the FAS predictive model.

Site parameters
In this section, we investigate on the site-to-site random terms δS2S s computed with respect to the mean site response of the 6 reference rock sites and whose dependencies can be studied to improve the site description of the final calibration. The statistical dependency is evaluated against geophysical proxies describing the effect of local site response due to deep and shallow geology, not already included as explanatory variables in the base FAS model: (i) the shear wave velocity in the uppermost 30 m, V S,30 (Borcherdt 1994) and (ii) the near-site high-frequency attenuation parameter 0 (Anderson and Hough 1984;Ktenidou et al. 2014). V S,30 is chosen as it provides quantitative and synthetic description of the subsurface structure (first 30 m in depth) that is most accessible to geotechnical investigations. It is also the most common explanatory variable typically included in ground motion models, as well as an easy and low-cost parameter to evaluate site response classification (e.g., Rodriguez-Marek et al. 2001;Abrahamson et al. 2008;Pitilakis et al. 2013).
Herein, the values of V S,30 (m/s) are taken from direct in-situ measurements of the S-wave velocity profile (provided by the Engineering Strong-Motion database, ESM https:// esm-db. eu/; Luzi et al. 2020;Lanzano et al. 2021) and the Italian database (ITACA). For stations where a measurement is not available, an estimate of V S,30 is provided from the empirical correlation with the topographic slope, according to Wald and Allen (2007) relationship, on the basis of slope measurements from high-resolution digital elevation models (DEM) of Italy (Mascandola et al. 2021). An histogram plot of these estimates of V S,30 is  Fig. 6a which denotes a log-normal distribution with a median value of about 900 m/s. Following the taxonomy proposed by Ktenidou et al. (2014), the values of 0 used for the model calibration are obtained from the amplification functions extracted from the M-GIT analysis according to the procedure by Drouet et al. (2011). The average site amplification of the 6 reference stations was assumed flat and equal to one in the M-GIT computation (Morasca et al. 2023), which is equivalent to assuming, on average, 0 =0 for the reference sites. As a result, the 0 estimates from the amplification functions suffer from this assumption and need to be corrected with the mean value of 0 (0.012 s) of the reference rock sites obtained from an independent computation on the observed FAS ( 0 -FAS in Table 1). The corrected values of 0 from GIT ( 0 -GIT in Table 1), in the following simply named 0 , are provided in the "Data availability" section and reported in Table 1 for the 6 reference rock sites. The histogram of 0 is shown in Fig. 6b where it follows an approximately normal distribution with a median value at 0.02 s.
We then examine the dependencies of δS2S s both with the V S,30 and 0 (Fig. 7) for lowand high-frequencies. The site empirical terms δS2S s exhibit a dependence on the V S,30 (anti-correlation) which is more evident at low-frequency (Figs. 7 and 8c) than at highfrequency (Figs. 7 and 8d), as also found by other authors (see Bergamo et al. 2020) who evidenced that this proxy shows larger correlation with site amplification in the low-tointermediate frequency range. On the contrary, a negative correlation between the site terms and 0 is more pronounced at high-frequency (Figs. 7 and 8b) than at low-frequency (Figs. 7 and 8a).

Regression approach
The assumed functional form of the advanced model is consistent with the one proposed by Bora et al. (2015), so that it is formalized as follows:  Table 1.
The distance scaling F R (R, M) remains unchanged compared to the base model, since we did not calibrate a specific relation between the path terms δP2P p and the parameters related to the propagation terms derived from the non-parametric M-GIT analysis. This has meant the inability to introduce in the parametrization an additional proxy describing the effect of anisotropic attenuation due to crustal heterogeneities over the study-area.
Finally, in Eq. (5.1), the terms δB′ e , S2S ′ s , L2L ′ r , P2P ′ p , W ′ 0 are the random terms of the advanced model, homologous to the base model.
A schematic summary of the explanatory parameters for the basic and advanced models is shown in Fig. 9.
In a preliminary stage of the analysis, we consider a first trial model including all the explanatory variables of Eq. (5.1) in a standard one-stage regression process. We apply this approach to explore the potential trade-off effects between the various terms of the model, which are well-known effects in the regression problems. This is a more relevant issue when the number of parameters in a regression increases, as in the present case; indeed, more complex models, while often associated with improved prediction capabilities, are also affected by larger trade-offs among the parameters that make it hard to resolve the individual contributions in the predictive equation.
Results of the one-stage trial regression analysis can be found in the "Data availability" section. In Fig. 10 we report the correlation matrix of the model parameters to investigate potential trade-offs and then assess the most convenient strategy of parametrization.
The matrix of the parameterized model evidences that there is a negligible correlation among the different groups of coefficients, which are associated to different physical effects (i.e. group b n : source-related; group c n : attenuation-related; group d n : site-related), as it can be noted in the example of Fig. 10 both at low-(Fig. 10a) and high-frequency (Fig. 10b). Yet, the source terms are affected by mutual trade-off (negative correlation) particularly among stress-drop (b 3 ), kappa source (b 4 ) and magnitude scaling for M < M h (b 1 ), as expected. A small cross-correlation is obtained between V S,30 (d 1 ) and site 0 (d 2 ): these two terms indeed mainly affect the model amplitudes in different frequency ranges, as shown before. Instead, a very strong correlation exists at high frequency (Fig. 10b) as well-known, between the geometrical spreading (c 1 and c 2 ) and the anelastic attenuation (c 3 ) so they should be always considered as two components of the same attenuation model (Bindi and Kotha 2020).

Stepwise parametric regression
In order to prevent the observed trade-off effects among the coefficients of the different contributions, we adopt a stepwise regression technique similar to the one proposed by Bayless and Abrahamson (2019). In this way, the choice of the predictive variables is carried out with a step-by-step procedure by which at each step, a new variable is considered for addition to the set of base explanatory variables, thus in order to constrain different components of the model using the data relevant to each piece (Bayless and Abrahamson 2018).
The procedure is repeated including the additional dependencies with the seismological and geophysical parameters for a total of 7 steps (Table 2) in order to build-up the base and the final advanced models, as described in the next section. At the end of each step, we also apply a smoothing on the coefficients, which is performed to get smoothed spectra and better constrain the model to more physical behavior . Table 2 shows the regression steps of the base model (i.e. without including the additional parameters) and the advanced model (i.e. including the additional parameters). For the base model (from step #1 to #3), we first constrain the magnitude scaling and the frequency-dependent geometrical spreading using data up to 80 km to isolate the effects of the geometric attenuation from the anelastic one, thus leaving the parameters a 0 , b 1 , b 2 , c 2 and the random terms free in the regression; then the c 2 coefficient is finally estimated and smoothed after the regression (step #1). Once fixed c 2 , in step #2, the magnitude-dependent geometric spreading and the anelastic attenuation are introduced by regressing also the terms c 1 and c 3 , respectively so the coefficients b 1 , b 2 , c 1 and c 3 are finally estimated and smoothed. Note that here the c 3 is forced to assume non-positive values in order to avoid unphysical increasing of ground motion at large distances (R > 80 km). At the end (step #3), the intercept a 0 of the functional and the random terms are estimated along with the corresponding standard deviations.
The same approach is used for the advanced model, (from step #4 to step #7) where we estimate the additional terms d 1 , b 3 , d 2, b 4 , a 1 and random terms by sequential stages. Note that the step model is built in such a way that the advanced model is a specialization of the base model; in fact, for this purpose, the coefficients of the advanced model (b 1 , b 2 , c 1 , c 2 , and c 3 ) are the same of the base model, so the former is calibrated assuming that the scaling with magnitude and distance of the latter are still valid. In addition, the reference values for the site parameters in the calibration (V S,30,ref and κ 0,ref ) are assumed as representative of the 6 reference sites in Table 1. In this way, if any of the additional terms of the advanced model were null, the returned predictions would be close to those provided by the base model. All the model coefficients are reported in the "Data availability" section.
The stepwise regression is also performed via a mixed-effects regression, thus considering the random-terms δB' e , δS2S′ s , δL2L′ r and δP2P′ p introduced in the base model that are zero-mean normally distributed random variables with variances 2 , 2 S2S , 2 L2L , 2 P2P , respectively. These residual components allow us to calculate the repeatable terms

#1
Base R < 80 km a 0 , b 1 , b 2 , c 2 + random effects c 2 #2 Base All a 0 , b 1 , b 2 , c 1 , c 3 (non-positive) + random effects b 1 , b 2 , c 1 , c 3 #3 Base All a 0 + random effects a 0 , standard deviations #4 Advanced All a 1 , d 1 (non-positive) + random effects d 1 #5 Advanced All a 1 , b 3 , d 2 (non-positive) + random effects b 3 , d 2 #6 Advanced All a 1 , b 4 (non-positive) + random effects b 4 #7 Advanced All a 1 + random effects a 1 , standard deviations and associated variability of event, site, source-area, and path and compare them to those obtained from the base model to evaluate the effectiveness of the parameterization. The plot of the raw and smoothed coefficients is reported in Fig. 11. One could note that the intercept value undergoes small variation when passing from the base model to the advanced one -i.e. the difference between the regression coefficient a 0 and a 1 (denoted by δa) is negligible especially below 10 Hz. Conversely, the corresponding offset obtained from the one-step regression (a′ in Fig. 11) shows a sensibly different trend compared to a 0 . This observation confirms that the stepwise approach allows to redistribute the model variability on the random terms and associated uncertainty rather than on the fixed-coefficients and thus it does not produce biased spectral shapes compared to the base model version.
In Fig. 12 (on the top), the coefficients b 1 and b 2 related to the magnitude scaling respectively for small earthquakes ( M w < M h ) and stronger earthquakes ( M w > M h ) tend to mainly affect the low-frequency range, whereas b 3 increases as frequency increases, thus reflecting the greatest relevance of the stress-drop parameter at high-frequencies.
Another aspect concerns the complementary trends of d 1 and d 2 (site effects), which are dominant in different frequency bands (Fig. 12 on the bottom). The former mainly act at low-frequency depending on the shear wave velocities at shallow depth, in line with the findings by Bergamo et al (2020) that show a good correspondence of V S,30 with site amplification across the whole 1.7-6.7 Hz band.
Instead, the coefficient d 2 prevails at high-frequency being dependent on κ 0 . Namely, d 2 is forced to zero up to 5 Hz to limit the contribution of the site attenuation parameter at high-frequency. The same approach is adopted to constrain the b 4 term that leads to linear effect of source only in the high-frequency range starting from about 5 Hz. For similar reasons, the c 3 coefficient of the apparent anelastic term is constrained to be non-positive up to about 3 Hz to avoid unrealistic effects, such as the enhancement of the ground motion from a certain distance onward, which is also a common assumption (e.g. Lanzano et al. 2021). The c 1 and c 2 coefficients (Fig. 12 in the middle) vary in the interval 0 to 0.15 and − 1.5 to − 2, respectively over the entire period range, with a typical trend observed in the geometrical spreading terms of the empirical models. Figure 12 also reveals that the coefficients of the stepwise regression and those related to the one-step calibration are quite similar, and that the only coefficients that show remarkably different trends are c 1 , c 2 and c 3 , confirming that a strong trade-off exists between these terms.
Greater uncertainties are observed for the coefficient b 2 (related to magnitude scaling at high magnitudes) due to poorer constraint to data in this range and d 1 (V S,30 scaling), probably as a consequence of the uncertainty of the V S,30 estimates.
To assess the impact of the additional seismological terms on the variability's reduction, we compare the individual uncertainty contributions of the base model against the advanced one (i.e. the parameterized model of Eq. (5.1)). Figure 13 shows that the between event component ( )- Fig. 13a -is significantly reduced in the advanced model with respect to the base model at high-frequencies, from about 0.25 to 0.08 log10 units for f > 20 Hz (total gain of 68%). Also the standard deviation associated to the location term -L2L Fig. 13b -vanishes at all frequencies except at f > 10 Hz, thus indicating that the variability contribution due to the source region is completely captured in the range 1-10 Hz by the introduction of the stress drop and source in the modelling.
The introduction of the geophysical variables (V S,30 and site kappa) have also a relevant impact on the site-to-site variability S2S as plotted in Fig. 13c, which shows a  Fig. 13d, due to the lack of a pathspecific proxy of crustal properties in the model description. A small reduction of variability is still observed as a result of the mutual trade-off between path and site effects, so that the introduction of the additional site parameters also reflects on the propagation terms. The remaining aleatory variability 0 of the parameterized model (Fig. 13e) is substantially equal to that of the base model, whereas the total sigma ( ) associated with the logarithmic FAS at each frequency, computed from the vectorial combination of all uncertainty contributions = √ 2 0 + 2 + 2 L2L + 2 S2S + 2 P2P , shows a stable trend with frequencies and an appreciable reduction at f > 30-40 Hz. As we can observe in Fig. 13f, this reduction reaches a maximum value of 0.4 log10 units at the highest frequency, due to the fact that the source and site parameters mainly affect the high-frequency range.
We point out that these values of sigma are relatively low compared to the most typical ones published in the literature (Akkar et al. 2014;2014. This aspect Fig. 13 Standard deviations of the random-effects terms of the nonergodic models against frequencies for: a between-event (τ) and location-to-location τ L2L ; b site-to-site S2S and path-to-path P2P ; c aleatory ( 0 ) and total sigma (σ) 1 3 indeed is one of the technical obstacles associated with GMM-FAS models as discussed in Bora et al. (2014). Our findings demonstrate a beneficial effect of the parametrization in constraining the uncertainty, which may be a result of the strategic combination of both an effective physical description of the ground motion through the introduced parameters and the application of a fully nonergodic framework for modelling.

Predictive scenarios
The median estimates of the FAS model (fixed-effects) are here shown to reproduce a set of predictive scenarios with varying magnitude and distance (Fig. 14). The overall trends for different magnitudes are similar to those obtained by Bayless and Abrahamson (2019) and reproduce the decreasing in the corner frequency with increasing magnitude as expected from theoretical source spectra with double-corner frequencies ω −2 model (Brune 1970) (Fig. 14a). Also the trend with distance is well captured by the median model, where is evident a stronger fall-off rate at high-frequency of the curve corresponding to 120 km in the range where c 3 is higher in absolute value, confirming that the trend is controlled by the anelastic attenuation at larger distances (Fig. 14b).
Other FAS scenarios are provided in Fig. 15 where the additional site and seismological parameters are allowed to vary while the other parameters are kept fixed. Variation in the Δσ (Fig. 15a) produces a significant effect both on the median prediction estimates and on the shape of the spectra; i.e., higher stress drop leads to greater FAS in the high-frequency range (according to b 3 increasing with frequency), as well as an increase in the higher corner frequency due to the theoretical scaling laws that link stress drop and corner frequency (Boore 1983). Note that in these examples we fixed all the other parameters to the mean values of each corresponding distribution; this causes the base model to coincide with the advanced model corresponding to the mean value of the stress drop (i.e. the curve corresponding to 1.5 MPa - Fig. 14a). This observation also clarifies that the advanced model leads back to the base one in case all the additional parameters are set to the mean values, as a result of the stepwise regression approach. The increasing source (Fig. 15b) and 0 at site (Fig. 15d) causes the FAS to drop proportionally in the high frequency range (> 5 Hz) as the respective parameter increases, consistently with the trend of b 3 and d 2 respectively, so without any effect at lower frequencies. This dependency of the FAS from the high-frequency attenuation parameter follows a wellknown and expected trend as observed by other authors (Lanzano et al. 2022a among the others).
In contrast, variations of V S,30 (Fig. 15c) show a pronounced effect on FAS at low frequencies and negligible at high frequencies. Even in this case, it can be noted that for values close to the average reference rock site condition in Italy (namely V S,30 = 1100 m/s and 0 =15 ms) and average source properties (Δσ = 1.5 MPa and source =20 ms), the advanced model returns to the base one.

Conclusions
In this work we have proposed a new FAS model for Central Italy constrained to seismological and geophysical parameters with the aim to better capture the underlying physics related to different energetic values of the source and different attenuation scaling due to the crustal properties and geological features of the region. Compared to other approaches, the proposed model takes advantage from the application of a fully nonergodic framework that provides improvements in terms of median predictions and in the constraining of the associated standard deviation.
Our study shows that the incorporation of additional proxies in the regression, such as the stress drop, the V S,30 and the high-frequency source and site kappa parameters, available thanks to the dense collection of data sampled in the area of the seismic sequences occurred in 2016-2017 in Central Italy, allows to capture a significant part of ground motion variability, whose reduction varies to a different extend with the frequency.
Moreover, the use of a stepwise regression technique has not only allowed to prevent the potential trade-offs arising from the inclusion of additional parameters, but it also represents a convenient and flexible strategy of parametrization to obtain the FAS predictions. Hence, the advanced model is built with a modular approach so it can be adopted even when any of the additional parameters is lacking, thus returning predictions close to those provided by the base model. This aspect has the advantage of providing a unique regression model that could be used in different contexts depending on the type and amount of parameters at disposal, so returning in a more complex functional shape with increased accuracy when all the physical parameters are known or, alternatively, with a simpler form when the physical data are missing. In particular, when the model is used in its advanced form for predictive purposes, the source parameters (stress drop and kappa source) could be derived from specific studies conducted on the source area and then used as input parameters with the corresponding variability to be handled as epistemic model uncertainty. Furthermore, it should be noted that, as a result of the regression analysis, the kappa term and the stress drop are highly correlated at lower frequencies, so a simplified version of the advanced model could be considered by including only one of these two terms for the source.
Results have shown that the impact of such a parametrization on the between-event variability, compared to the base model, only including dependency on magnitude and distance, corresponds to a reduction of about 68% at high-frequencies. Also the standard deviation associated to the location term goes to zero at almost all frequencies, indicating that the variability contribution due to the source region is completely captured under 10 Hz by the introduction of Δσ and source in the regression.
Finally, the inclusion of the variables V S,30 and 0 has also a relevant impact on the site-to-site variability, which shows a reduction of about 35% compared to the counterpart contribution of the base model, particularly at high-frequency.
Unchanged values of the remaining aleatory variability indicate that all the effects captured by the advanced parametrization concern only repeatable effects related to the source and site that are quantified as random-terms in the nonergodic base model. Other non-systematic physical phenomena, such as those linked to rupture directivity, slip distribution, radiation pattern etc. are not accounted for here and will deserve further analysis, for instance using the same approach implemented by Colavitti et al (2022) to model directivity on the same target area, therefore to combine those findings with the present parameterization and get a more comprehensive modelling.
With regard to the total variability of the advanced model (0.35-0.45 log10 units), we observe similar or even lower values (especially at high-frequency) with respect to those available in the literature, such as: 0.28-0.5 log10 for the parameterised model of Bora et al. (2015); 0.3-0.52 log10 for Bayless and Abrahamson (2019); 0.35-0.65 log10 units for Kotha et al. (2022); 0.41-0.61 log10 for the model of Sung et al. (2022) for France.
Finally, we suggest that the final FAS model could be coupled with Random Vibration Theory (RVT, Crandall and Mark 1963) methods to convert the median nonergodic predictions into pseudo-spectral accelerations and then used to generate physics-based shaking scenarios, as well as to produce expeditive simulations of ground-motion time series which may help to improve the coverage of recorded data in near-source region and increase understanding of the spatial features of ground shaking within this boundary.
Funding Open access funding provided by Istituto Nazionale di Geofisica e Vulcanologia within the CRUI-CARE Agreement. This research is supported by Istituto Nazionale di Geofisica e Vulcanologia (INGV) in the frame of the project Pianeta Dinamico (Working Earth) -Geosciences for the Understanding of the Dynamics of the Earth and the Consequent Natural Risks (CUP code D53J19000170001), founded by the Italian Ministry of University and Research (MIUR) in the Task S3 -2021 -Seismic attenuation and variability of seismic motion. This study was also partially funded in the framework of INGV and Dipartimento della Protezione Civile (INGV-DPC) agreement B1 2019-2021, with the aim of promoting research activities in the field of seismic hazard in Italy.

Data availability
The seismological parameters and amplification functions derived by the GIT inversion are available in the repository CI-GIT (Morasca et al. 2022). The parametric table of the Fourier Amplitude Spectra ordinates are available in the repository CI-FAS_Flatfile . The datasets of model coefficients developed in the current study are available in the repositories CI-FAS_GMM_one (Lanzano et al. 2022b) and CI-FAS_GMM_step (Lanzano et al. 2022c) for one-step and stepwise regressions, respectively. All the repositories are accessible at the following link: https:// shake. mi. ingv. it/ centr al-italy/.