Spectral decomposition of the Engineering Strong Motion (ESM) flat file: regional attenuation, source scaling and Arias stress drop

We perform a spectral decomposition of the Fourier amplitude spectra disseminated along with the Engineering Strong Motion (ESM) flat file for Europe and Middle East. We apply a non-parametric inversion schema to isolate source, propagation and site effects, introducing a regionalization for the attenuation model into three domains. The obtained propagation and source components of the model are parametrized in terms of geometrical spreading, quality factor, seismic moment, and corner frequency assuming a ω2 source model. The non-parametric spectral attenuation values show a faster decay for earthquakes in Italy than in the other regions. Once described in terms of geometrical spreading and frequency-dependent quality factor, slopes and breakpoint locations of the piece-wise linear model for the geometrical spreading show regional variations, confirming that the non-parametric models capture the effects of crustal heterogeneities and differences in the anelastic attenuation. Since they are derived in the framework of a single inversion, the source spectra of the largest events which have occurred in Europe in the last decades can be directly compared and the scaling of the extracted source parameters evaluated. The Brune stress drop varies over about 2 orders of magnitude (the 5th, 50th and 95th percentiles of the ∆σ distribution are 0.76, 2.94, and 13.07 MPa, respectively), with large events having larger stress drops. In particular, the 5th, 50th and 95th percentiles for M > 5.5 are 2.87, 6.02, and 23.5 MPa, respectively whereas, for M < 5.5, the same percentiles are 0.73, 2.84, and 12.43 MPa. If compared to the residual distributions associated to a ground motion prediction equation previously derived using the same Fourier amplitude spectra, the source parameter and the empirical site amplification effects correlate well with the inter-event and inter-station residuals, respectively. Finally, we calibrated both non-parametric and parametric attenuation models for estimating the stress drop from the ratio between Arias intensity and significant duration. The results confirm that computing the Arias stress drop is a suitable approach for complementing the seismic moment with information controlling the source radiation at high frequencies for rapid response applications.


Introduction
In the context of probabilistic seismic hazard assessment (PSHA), ground motion prediction equations (GMPEs) allow one to compute, for any seismic scenario of interest, the probability of exceeding thresholds for parameters of engineering of interest, such as peak ground acceleration or elastic response spectral accelerations at different periods (Cornell 1968;McGuire 1976). GMPEs are usually derived from regression analysis performed over observed strong motion data, although outcomes from numerical simulations can contribute to improve the explanatory power of the models (Douglas and Edwards 2016). To support the development of GMPEs, parametric tables (also known as 'flat files'), including event and station metadata along with the values of several strong motion parameters of engineering interest, have been disseminated in the framework of several projects. Examples are the NGA-West2 flat file (Ancheta et al. 2014) for the development of GMPEs at global scale (i.e., including recordings from multiple regions distributed worldwide) or the recent Engineering Strong Motion flat file assembled for Europe and Middle East (Lanzano et al. 2018). The predictions from a GMPE describe a normal distribution characterized by the mean value expected for the considered seismic scenario and the standard deviation, which quantifies the spread of the predictions around the mean due to apparent aleatory variability. In the framework of GMPE development and application, several topics are of particular interest such as the study of the aleatory variability and its decomposition into event, station, and path specific components; the characterization of the regional dependence of ground-motion models; the treatment of the epistemic uncertainty through logic tree or back-bone strategies (Douglas 2018). For example, several studies highlighted the correlation existing between the event specific component of the residual distribution (the so called inter-event or between-event residual distributions) and the variability of the stress drop, which is generally estimated from the corner frequency of the source spectra (e.g., Bindi et al. 2007Bindi et al. , 2019aAmeri et al. 2017;Baltay et al. 2017;Oth et al. 2017;Trugman and Shearer 2017). Although the stress drop is a key source parameter controlling the shaking at high frequencies, it is not considered among the explanatory variables of GMPE because this would require stress drop models to be provided for the source models used in PSHA computations whereas constrained stress drop distributions are available only for well monitored regions. Moreover, the possibility to jointly analyze distributions derived for different areas is hampered by the strong dependence of the stress drop values on methodological assumptions.
In this study, we exploit the Fourier amplitude spectra disseminated along with the ESM flat file to develop seismological models for describing source, propagation and site effects using the same data set prepared for the development of GMPEs used to update the PSHA in Europe (Weatherill et al. 2020)., These seismological models can support the GMPE developers in deciding the regionalization, in taking decision about the functional form and in interpreting the residual distributions; furthermore, these models quantify the variability of key seismic parameters, such as the anelastic attenuation and the stress drop, needed for the generation of synthetic ground motion through stochastic simulations (Douglas and Aochi 2008;Edwards and Fäh 2013;Drouet and Cotton 2015;Bora et al. 2017). This study is organized as follows: we first present the dataset used and the methodology applied to decompose the Fourier spectra into source, propagation and site effects; then, we discuss the results in term of non-parametric models, which are later interpreted in terms of parametric models introduced to describe the spectral attenuation and the source spectra, and we present the scaling of the source parameters with the earthquake size. Finally, the recently proposed approach for estimating the stress drop from the Arias intensity (Baltay et al. 2019) is investigated, generalizing the approach used for correcting attenuation and site effects.

Data set
The Engineering Strong Motion flat file (ESM-ff) is a parametric table including information of interest for engineering seismology, relevant to earthquakes occurred in Europe and Middle East since 1970s (Lanzano et al. 2018). ESM-ff mainly includes events with moment magnitude larger than 4 (in the following written using the notation M4 +), although few earthquakes with magnitude in the range M3-4 are also considered (Lanzano et al. 2019). Since moment magnitudes from different catalogs are available in ESM-ff, in the reminder of this article we always refer to the EMEC moment magnitude (Grünthal and Wahlström 2012) which is available for all considered earthquakes. Along with several intensity measures, such as peak parameters, integral parameters, spectral accelerations, Fourier amplitudes, ESM-ff provides station and event metadata needed to derive, for example, ground motion models for seismic hazard assessment. Figure 1 exemplifies few selected ESM-ff characteristics of interest for this study, considering events shallower than 50 km. The left panel shows the magnitude versus distance scatter plot relevant to 20,026 recordings generated by 1372 earthquakes recorded by at least 2 stations; in the central panel, the scaling with distance of Fourier amplitudes computed for earthquakes with M4.5 (± 0.1) and M6 (± 0.1) are compared; in the right panel, the cumulative number of recordings per event and per station are shown. About 20% of stations and 30% of events have one recording whereas 60% of the stations and 50% of the events has at least three recordings. In this study, we analyze the geometrical mean of the two horizontal spectra; in ESM-ff, the FAS are computed considering the whole time histories and we assume that they are mainly representing the contribution of S-waves. A detailed description of ESM-ff is available in Lanzano et al. (2019) and Bindi et al. Fig. 1 Left: moment magnitude versus distance scatter plot of earthquakes shallower than 50 km and recorded by at least two stations. Middle: Fourier amplitude spectra versus distance of recordings for moment magnitude in the range 4.5 ± 0.1 (white dots) and 6 ± 0.1 (black dots). Right: empirical cumulative density function for the number of recordings available either per station (black) or per event (gray) 1 3 (2019b). From ESM-ff, subsets suitable for the spectral analysis developed in this study are extracted, as explained in the following sections.

Spectral decomposition
We apply a spectral decomposition approach to isolate source, attenuation and site amplification contributions to the measured spectral amplitudes. We implement the so-called Generalized Inversion Technique (GIT) (e.g. Castro et al. 1990;Edwards et al. 2008;Oth et al. 2011) following a non-parametric approach where the redundancy in the data set (i.e., the same event recorded at different distances by different stations and each stations recording several events) allows us to estimate, in a least-squares sense, source, propagation and site terms without imposing to them any a priori functional form. Under the assumption of linear system, the GIT approach is based on the following convolution model: where logFAS ij (f) is the logarithm in base 10 of the spectral amplitude at frequency f of event i recorded by a station j located at an hypocentral distance equal to R. We discretize the distance ranges from 1 to 100 km, from 100 to 200 km and from 200 to 300 km into equal spaced bins 5, 10, and 20 km wide, respectively. In Eq. (1), R n are the nodes separating the distance bins and n is selected such that R n ≤ R<R n+1 ; S k (f) and Z k (f) are the source and site spectral amplitudes at frequency f, respectively, and δ ik , δ jk are Kroneker's deltas used to select the correct source and site for the record i, j when constructing the linear system. The spectral attenuation A region (f;R n ) is a non-parametric model (i.e. provided in tabular form) evaluated at nodes R n , being the model region dependent. A linear interpolation is applied between the nodes by setting a n = (R n+1 − R)/∆R, ∆R = (R n+1 − R n ), a n+1 =1 − a n . The linear system in Eq. (1) is solved in least-squares sense (Koenker and Ng 2017). (1) Polygons defining the three regions considered to derive the attenuation models: IT (cyan), GT (green) and RE (red). The epicenter (white circles) and station (blue triangle) locations are shown as well, considering only those with at least two (left) or five (right) recordings Previous studies highlighted significant variations in attenuation models derived across Europe (e.g., Kuehn and Scherbaum 2016;Kotha et al. 2017;Bindi et al. 2019c). Therefore, taking into consideration the distribution of data analyzed in this study, we estimate different attenuation spectral models A region (f;R n ) for three regions: the first one (hereinafter referred to as IT) is including mainly Italy; the second region (named GT) encompasses the Aegean and Anatolia regions; the third region (named RE) includes the rest of Europe, but with data mostly located in the Balkans and in Romania (Fig. 2). The three non-parametric attenuation models are jointly derived by solving Eq. (1). It is worth noting that any model based on the linear combination of three terms as in model (1) is inherently affected by trade-offs among the three different solution blocks, since there are two unresolved degrees of freedom. Therefore, we add two a priori constraints to Eq. (1). Considering that we follow a non parametric approach and we cannot constrain parameters associated to physical models, we require that the attenuation for all regions is 1 at the reference distance R ref = 10 km (i.e. the source spectra are scaled at 10 km) and we set to zero the logarithmic average of the site effects computed for a set of reference stations. It is worth remembering that the trade-off between source and site terms are expected not to affect the attenuation model (Oth et al. 2011) allowing to get the same attenuation models for different choices of the reference site conditions. Finally, in order to smooth the attenuation models, we require that the second derivative of the attenuation curves is small (Castro et al. 1990).
We run the GIT inversion twice; in the first inversion, we retrieve the spectral attenuation with distance whereas, in the second inversion, we estimate the source spectra and the site amplifications. The second step solves the following system of equations: where that is, the attenuation models estimated in the first inversion are used to correct the spectra (Eq. 3) which are, in turn, used to split the source and site contributions (Eq. 2). To remove the trade-off between source and site, in the second inversion we constrain to zero the logarithm of the average site amplification of a subset of stations, where set of reference stations is selected during the analysis. Inversions (1) and (2) are performed by assuming different reference site conditions and by performing different data selection, as detailed in the next paragraphs.

Results: spectral attenuation models
To solve Eq. (1), we select stations and events with at least 5 recordings (Fig. 2). The obtained spectral attenuation curves are shown in Fig. 3 and provided in Tables from ESM1 to ESM3 of the online Resources.
The attenuation in the three regions is characterized by different rates of attenuation and by different dependencies on frequency. For example, the spectral attenuation curves for region RE show in the first 35 km a narrower spread than for the other two regions. Moreover, the attenuation in the three regions show different changes in slope, reflecting different (3) logFAS ij (R;f ) = logFAS ij (R;f ) − a n logA region f ;R n − a n+1 logA region f ;R n+1 contribution of later arrivals controlled by the crustal structure and hypocentral depths. To quantify the differences among the three regions, a parametric model composed by geometrical spreading and anelastic attenuation terms is fit to the non-parametric solutions: where the geometrical spreading G(R) is assumed to be frequency independent and U(f, R) accounts for the anelastic attenuation. Considering the transformation d = log(R/R ref ), the geometrical spreading model is described by a piece-wise linear function with two breakpoints d 1 and d 2 and three different slopes s 1 , s 2 , s 3 , that is: The anelastic attenuation is modeled in term of quality factor Q(f) = Q o f α , that is: where we fix the path-average shear-wave velocity β to 3.6 km/s and R ref to 10 km.
A multiple steps procedure is followed to calibrate G(R) and U(f, R). Since the geometrical spreading is assumed to be frequency independent, we first consider the spectral attenuation at 1 Hz to determine G(R) through a breakpoint regression (Muggeo 2003(Muggeo , 2008. To make more robust the assessment, the average over 4 spectral attenuation curves between 0.9 and 1.15 Hz is used to represent the non-parametric attenuation at 1 Hz. A preliminary fit is performed by considering a geometrical spreading term without breakpoints multiplied by an exponential term given by Eq. (6) evaluated at 1 Hz. The obtained Q 0 is used to correct the data for the anelastic attenuation effects at 1 Hz and then a breakpoint regression is performed over the Q 0 -corrected spectral values, considering the geometrical spreading model as in Eq. (5). Finally, the spectral attenuation curves corrected for the piece-wise linear geometrical spreading are used to Non-parametric attenuation with hypocentral distance for regions RE, IT, and GT (from left to right, respectively). The attenuation curves for frequencies between 0.2 and 30 Hz are shown in gray whereas the curves in color are the spectral attenuation curves averaged over selected frequency intervals as indicated in the legends. The 1/R and 1/R 2 attenuation rates are also shown for comparison assess Q(f) over the frequency range from 1 to 20 Hz. The obtained parameters of models in Eqs. (5) and (6) are listed in Tables 1, 2 and 3. The comparison between non-parametric and parametric models is exemplified in Fig. 4, for three selected frequencies (i.e., f = 2, 8, and 15 Hz). A general good agreement is observed at low frequencies whereas, at high frequencies, the parametric models fail in describing the spectral attenuation for distances above 120 km. For region RE, the parametric model attenuates more than the non-parametric at distances shorter than 35 km. The parametric models are summarized in Fig. 5 in terms of geometrical spreading and frequency-dependent quality factors. Region GT is characterized by the weakest geometrical attenuation with distance but coupled with the strongest anelastic attenuation; regions IT and RE show similar geometric spreading rates, although with breakpoints located at different distances, with stronger anelastic attenuation for IT (see also Tables 1, 2, 3). It is worth remembering that the geometrical spreading and anelastic attenuation models are affected by strong mutual trade-off Left: geometrical spreading attenuation term; right: quality factor Q(f) and they should be always considered as two components of the same attenuation model. Moreover, the standard errors listed in Tables 1, 2, and 3 are conditional on model assumptions that allowed to obtain the non-parametric attenuation curves.
In conclusion, the parametric models are based on strong simplifications and they should be considered as first order approximation of the non-parametric results. In particular, the parametric models derived in this study should not be used for simulating the ground motion at high frequencies (i.e. above 10 Hz) at distances larger than 120 km; in order to capture the high frequency attenuation at longer distances, Fig. 4 suggests that spectral decomposition methods based on parametric approaches should implement more complex parametrizations.

Results: non-parametric source and site spectral models
The source and site spectra are extracted from the FAS corrected for the attenuation models logA region (f, R) (Eqs. 2, 3). Since there are groups of stations and events that are weakly connected when the data set is analyzed as a whole, in order to get stable solutions we split the data set into two subsets, the first covering mainly Italy (hereinafter referred to as DS1) and the second one including the Aegean sea and western Anatolia (DS2). After trial regressions, for DS1 we select stations with at least 4 recordings and events recorded by at least 3 stations whereas, for DS2, we analyze events and stations with at least 5 recordings. The source-station geometries for D1 and D2 are shown in Fig. 6.
In solving Eq.
(2), we remove the trade-off between source and site by constraining to one the average site amplification for a group of stations for both DS1 and DS2. In order to select the set of reference stations, for each data set we run a preliminary inversion by constraining the average amplification over all stations to one. Then, we evaluate the percentiles of the distribution of the spectral amplifications for each station (Fig. 7) and we select those stations satisfying the following conditions: p50 < 0.9, p95 < 1, p05 > 0.1, |p95-p05| < 0.4 where p95, p50, and p05 are the 95th, 50th and the 5th percentiles of the  As expected, the amplifications of stations contributing to the reference site condition (black curves) are almost flat with amplitudes within a factor 2 from their ensemble average. Regarding the source spectra, since they are derived within the same spectral decomposition, it is possible to compare directly the results of the largest events. Since the monitoring networks developed over the last decades, the large earthquakes occurring earlier than 2000 were recorded by a number of analog stations operating on trigger mode smaller than the number of digital recordings available for the more recent events. For example, the number of recordings that we selected from ESM-ff to analyze the M6.4, 1976 Friuli and M6.9, 1980 Irpinia earthquakes are 6 and 4, respectively; for the 1999 Izmit and Düzce earthquakes, the number of analyzed records are 12 and 9, respectively, whereas for the more recent 2009 L'Aquila and 2016 Norcia earthquakes, we used 50 and 175 records, respectively.

Source parameters and source scaling
To extract the source parameters, a source model is fit to the non-parametric source spectra. For the acceleration spectra, we consider the following ω 2 model (Brune 1970(Brune , 1971 where Mo and fc are, respectively, the seismic moment and the corner frequency defining the ω 2 model; R θϕ is the S-wave radiation pattern (fixed to 0.55); F is the free surface amplification factor (fixed to 2); ρ is density (fixed to 2700 kg m −3 ); β is the shearwave velocity (fixed to 3200 m/s considering the average velocity at the events hypocentral depths and the tomographic model derived by Yang et al. 2018). In Eq. (7), the exponential term accounts for deviations of a flat acceleration source spectra at high frequencies as predicted by the ω 2 model. Whereas relative variation of the slope k source among different earthquakes account both for energy radiation effects at high frequency and attenuation Fig. 9 Non-parametric acceleration source spectra obtained for DS1 (left) and DS2 (right) regions (see Fig. 2). The spectra for selected large events are shown in color effects close to the fault, the average slope of the source spectra can be ascribed to the site constraint applied to the linear system since a near surface attenuation effect common to all reference stations not accounted for by the reference site condition is transferred to the source spectra (e.g., Oth et al. 2011;Bindi et al. 2018). In Eq. (7), f k indicates the frequency from which the exponential is applied. Therefore, k source is set to zero for f < f k , while determining f k during the fit. The reference site condition applied to solve Eq.
(2) influences also the average source spectral amplitude, which is in trade-off with the a priori fixed average site amplification. Therefore, in order to get seismic moments compatible with the EMEC moment magnitude, we run a preliminary fit to estimate the average bias of the seismic moments with respect to EMEC. It is worth remembering that within EMEC, complex hierarchical regional-schemes are implemented to convert values from other magnitude scales to Mw when the seismic moment is not available. Therefore, we anchor on average the source spectra to EMEC and then we estimate the seismic moment of each earthquake by performing a spectral fitting. The average bias is computed considering earthquakes in the range M3.8-5.5, to avoid band-limited effects. An average offsets of 0.11 magnitude units obtained for DS1 is applied to both DS1 and D2S source spectra. The offset compensated spectra are fit again, constraining the seismic moment of events with M5.5 + to the EMEC values.
To mitigate the effect of trade-offs existing among different source parameters, the fit is performed through the following steps: (1) the first fit is performed for the displacement spectra in order to emphasize the contribution to the misfit of the low frequencies (seismic moment controlled) with respect to the high frequencies (k source controlled); (2) in the second step, the fit is repeated but integrating the spectra to velocity and constraining the seismic moment to the value provided by the fit in displacement; (3) for the last step, the spectra are integrated to acceleration to determine f k and k source , with both seismic moment and corner frequency fixed to the values obtained in the previous steps. The value of f k is fixed to 0 and 10 Hz in the first two fits whereas a grid search procedure is applied for the fit in acceleration, selecting the value of f k > f c which corresponds to the minimum misfit solution. For earthquakes with EMEC magnitude larger than 5.5, the seismic moment is fixed to the value corresponding to the EMEC moment magnitude and the displacement fit is skipped. For M5.5 +, f k in the velocity fit is fixed to 5 Hz. Figure 10 exemplifies the fitting procedure showing the results obtained for three different earthquakes. Fig. 10 Example of fit for three earthquakes with moment magnitude 3.25, 5, and 6.5. From left to right, the fit is performed in displacement, velocity and acceleration, respectively. Black curves are the source spectra retrieved from the GIT inversion, the red curves are the best fitting models The obtained Mo and fc are used to compute the stress drop ∆σ assuming a circular rupture model with uniform stress drop (Eshelby 1957) where the rupture radius is given by Figure 11 shows the source parameters obtained for the analyzed earthquakes: in the left panel, the moment magnitudes derived from the seismic moment are compared to EMEC whereas the overall scaling of the seismic moment versus corner frequency is shown in the right panel(and the individual source parameters are reported in Table ESM4 of the online Resources). It is worth noting that the assumption of a circular model could lead to an overestimation of the stress drop for the largest events. For example, for the 1999, M7.4 Izmit earthquake we obtain a source radius r = 12.26 km which, using Eq. (8), corresponds to ∆σ = 37.59 MPa. Using the empirical relationships of Thingbaijam et al. (2017) for a strike slip event with M7.4, we obtain average fault length and width equal to L = 124 km and W = 24 km, respectively. Since the length scales controlling Δσ for the actual rupture are expected to be shorter than the length or width of the entire rupture, following Mai and Beroza (2000) we apply a reduction factor of 0.75 and 0.9 for L and W, respectively. Assuming an elliptical slip distribution with semi-axis a = L/2 and b = W/2, the stress drop is given by (e.g., Denolle and Shearer 2016): where η is a function of the Poisson's modulus η (here assumed equal to 0.25), of the aspect ratio a/b and of the complete elliptic integrals of first and second kind (Eshelby 1957, Eq. 5.3;Denolle and Shearer 2016, equation B.8). Using Eq. (10), we get ∆σ = 9.05 MPa whereas considering an elongated rectangular rupture area (Aki 1972 Fig. 11 could be biased toward higher values, it is worth to remember that simulation studies using stochastic approaches implement source models often based on the circular rupture assumption. Baltay et al. (2019) proposed to estimate the stress drop from the Arias intensity Ia. Since Ia is computed from the squared acceleration integrated over the significant duration, it is related to the root mean square acceleration and therefore it can be connected to the stress drop (Hanks 1979;Hanks and McGuire 1981). We use the estimates of Ia and of the 5-95% Arias intensity duration T 90 provided by ESM-ff (see Baltay et al. 2019 for a discussion about the effect of considering different duration estimates). The stress drop derived from Ia (hereinafter referred to as ∆σ arias ) is computed as the average, over the recording stations, of the station-specific estimates obtained after correction for propagation effects, that is where i is the event index, j the station index, < > j indicates the average over the set of stations that recorded the event i, R ij is the rupture (for M5.5 + events) or hypocentral distance. The quantity Y ij measured at station j for event i is given by (Baltay et al. 2019):

Arias stress drop
where the distance dependence of Y ij is carried out by the scaling with distance of the Arias intensity to duration ratio. In Eq. (12), we fix β = 3.2 km/s, ρ = 2800 km/m 3 , R θϕ = 0.55, f max = 25 Hz, whereas Mo is the seismic moment. In this study, we test three different choices for the propagation term P(R ij ) in Eq. (11), that is: Approach 1-Baltay et al. (2019) corrected Y ij for the geometrical spreading term, neglecting both the contribution of anelastic absorption and the site amplifications. They used a geometrical spreading equal to the inverse of the distance, which correspond to P(R) = R 6/5 . The authors showed that for large magnitude events (i.e., M6 +), the distance dependent ratio (Ia/T 90 ) 3/5 scaled similarly to R −1.2 . To support their assumptions also for smaller events, they applied a strict data selection considering M4.5 + ; setting the maximum distance R max to 20 km for magnitudes up to 5.5, R max = 100 km above M6.5, and linearly increasing from 20 to 100 km for M5.5-6.5; using a minimum distance equal to half of the Boore and Thompson (2014) finite fault factor h. We apply the same correction but, (2 ) 2 6 5 1 8.47Mo i 1 5

Ia T90
3 5 ij in order to slightly enlarge the data set, we consider M4 + and R max = 140 km for M7 + , keeping the linear increase of R max from 20 to 140 km for M increasing from 4.5 to 7.
Approach 2-In order to account for the anelastic attenuation and for site effects, and to avoid the application of a strict data selection, we calibrate a non-parametric attenuation model for Y ij by solving the following system of equations: which is similar to (1). In Eq. (13), the unknown are the attenuation table logA arias (R n ;R ref ), where the attenuation is set to 1 at the reference distance R ref = 10 km and R n are the discrete distances from 5 to 300 km, with step 3 km; the non-parametric station corrections Z np for each station j with j = 1, …, N station ; the source parameters Δ̃ k for each earthquake k with k = 1, …, N event . The average of the station corrections Z np is constrained to zero. Given Δ̃ k from (13), we compute the Arias stress drop ∆σ arias as follow: The non-parametric model logA arias (R n ;R ref ) is provided in Table ESM5 of the online Resources.
Approach 3-Similar to approach (2), we calibrate an attenuation model for Y ij but we consider a parametric model based on a functional form which includes a magnitudedependent attenuation term, that is: where Z p are the parametric station corrections for the Arias stress drop (constrained to have a zero average over the stations) and (c 1 , c 2 , c 3 ) are the coefficients defining the parametric attenuation model, being c 3 the coefficient controlling the magnitude dependent geometrical spreading. The obtained coefficients are c 1 = (3.86 ± 0.06) c 2 = (0.0027 ± 0.0001) c 3 = (0.37 ± 0.01). Figure 12 compares the non-parametric attenuation model of Eq. (13) (green curve) with the parametric attenuation model of Eq. (15) evaluated for M5 and M7 (red curves), and with the geometrical spreading term (10/R) 6/5 (dashed curve). Figure 12 also shows the GIT non-parametric spectral attenuation (from Eq. 1) of region IT for 1 and 10 Hz (black curves), and the piece-wise linear geometrical spreading (blue curve). From Fig. 12 we observe that the Arias intensity for moderate events (M5) attenuates as the high frequency spectral amplitudes whereas the Arias attenuation for large events (M7) is closer to the GIT attenuation model obtained for low frequencies; the non-parametric Arias attenuation model (magnitude independent) better describes the high frequency attenuation which is typical of moderate size events. Regarding the geometrical spreading, it underestimates both the non-parametric and parametric attenuation models. Figure 13 compares ∆σ arias obtained considering the three above mentioned approaches with the stress drop ∆σ computed from the corner frequency analysis (Eq. 8). ∆σ and ∆σ arias show a moderate correlation, being the Pearson correlation coefficient (Davis 1986 ik logΔ̃ k + a n logA arias R n ;R ref computed between log(∆σ) and log(∆σ arias ) equal to 0.78, considering the parametric approach of Eq. (15) and M6.5-. The strongest correlation is observed for the parametric model (Eq. 15) (panel a) whereas the non-parametric approach (panel b) generates larger ∆σ arias for the largest events, as expected by the behavior of the attenuation models shown in Fig. 13 (i.e. the non-parametric model for Arias is magnitude independent and better describe the attenuation of moderate events). The Arias stress drop computed correcting for the inverse of the distance (panel c) shows a similar pattern to the parametric model but the applied data selection reduces the number of events significantly. The bottom panels of Fig. 13 show the magnitude dependence of the ratio γ = ∆σ arias /∆σ. Whereas for the parametric model (panel d) γ is on average negative (being the mean ratio equal to 0.96), the ratio computed for the non-parametric model (panel e) shows a positive trend for M5.7and an almost constant positive average ratio equal to 1.06 for larger magnitudes. The scaling with seismic moment of the parametric ∆σ arias is shown in Fig. 14 (top panel). Below M5.5, ∆σ arias values shows a large scatter with an average positive trend with seismic moment; above Mw 5.5, ∆σ arias values are constant, on average (∆σ arias = 6.6 ± 0.2 MPa). Finally, the bottom panel of Fig. 14 shows that the Arias station corrections obtained through the parametric and non-parametric calibrations (Z np and Z p in Eqs. (13) and (15), respectively) are in very good agreement and, for several stations, the station correction to log(∆σ arias ) is large (e.g., the standard deviation of Z p is 0.43).

Fig. 12
Comparison between the spectral attenuation at 1 and 10 Hz (black curves), and the attenuation models for the Arias stress drop (parametric in red, non-parametric in green). The piece-wise geometrical spreading model and the inverse of distance geometrical spreading (see Eq. 5) are also shown

Discussion
In the framework of engineering seismology, it is common practice to compile parametric tables (referred to as flat files) to disseminate information about earthquakes, recording sites, and strong motion parameters of engineering interest. Examples are the flat file for the development of the Next-Generation ground motion models disseminated by PEER (Ancheta et al. 2014), or the flat files compiled for Japan (Dawood et al. 2016) and for Chile (Bastías and Montalva 2016). In this study, we focused on the Engineering Strong Motion flat file (ESM-ff) disseminated for hazard studies in Europe (Lanzano et al. 2018). The flat files are generally compiled for spectral ordinates but ESM-ff is also disseminating Fourier Amplitude Spectra (FAS), allowing us to derive spectral models for source, propagation and site amplifications using the same data used to develop GMPEs for the PSHA update in Europe (Weatherill et al. 2020). We applied a spectral decomposition approach to isolate the three terms without assuming a priori parametric models, letting the data to shape the source and site amplification spectra and to determine the frequency and distance dependent scaling. The advantage is that the solutions are not bounded to specific trends prescribed by the parametric models and alternative parametric models can be tested to describe any non-parametric component obtained after the spectral decomposition without impacting on the others two. Following previous study, such as the spectral decompositions performed for the Vrancea subduction (Oth et al. 2008) or for Japan (Oth et al. 2011), we introduced a regionalization of the attenuation, considering three different domains. The results confirm previous findings obtained for Europe by analyzing either response spectra (Kotha et al. 2016) or high-pass filtered peak displacement for local magnitude calibration (Bindi et al. 2019c): all studies found that the attenuation for the region including Italy is stronger than for the rest of Europe. Moreover, the non-parametric attenuation models capture the modulation of the attenuation with distance and frequency related to the contribution of later arrivals generated by crustal discontinuities, such as the effect of critical reflections from Moho previously highlighted for northern Italy (Bragato and Tento 2005) and northwestern Turkey (Baumbach et al. 2003;Bindi et al. 2006). We described the non-parametric attenuation in terms of parametric models including both geometrical spreading and anelastic attenuation. We allowed the geometrical spreading to change slopes at two breakpoints located at different distances for the different regions, and we moved the frequency dependency to the anelastic term. It is worth remembering that the frequency dependency of Q is still a topic of debate (Mitchell 2010) and frequency-independent but depth-dependent Q models could, in some cases, describe the observations as well (Zollo et al. 2014). The best fit parametric models provide a good first order description of the spectral attenuation that can be used for different tasks such as to develop stochastic simulations. Anyway, the parametric model cannot capture the full range of features shown by the non-parametric curves (e.g., the weak frequency dependence of the attenuation in region RE for distances shorter than 35 km) and tend to over-estimates the attenuation at high frequencies for distances above 120 km.
The spectra corrected for propagation effects were inverted to isolate the source contribution from the local site amplification effects. Since the data set includes stations installed over a wide spectrum of geological conditions, several stations show large amplifications over different frequency ranges. Figure 8 shows that amplifications exceeding 5, or even 10, at both low (e.g. station RTI, Rieti) and high (e.g., station NCR2, Nocera Umbra) frequencies. Considering the cases highlighted in Fig. 8, station GBP (Gubbio Piana, vs30 = 224 m/s) is installed in the middle of the Gubbio basin, a northwest-southeast 20 km long and 6 km wide basin in central Apennines, filled by fluvial lacustrine sediments with a thickness reaching 600 m. In the Gubbio basin, 2D-3D site effects generate large amplifications over a broad frequency range, starting from about 0.3 Hz and affecting both the horizontal and vertical components (Bindi et al. 2009). Station CLF (Colfiorito, vs30 = 145 m/s) is installed in the Colfiorito basin, another intermontane basin located in central Apennines, with an extension of about 3 km and a complex bedrock topography. The amplifications at 1 Hz are connected to edge diffracted surface (Love) waves generated by a 180 m deep depression (Di Giulio et al. 2003). Amplifications at station NCR2 (vs30 = 555 m/s, class E of Eurocode8; CEN 2004) are controlled by a buried wedge of weathered rock (Rovelli et al. 2002) Fig. 15 Top. Comparison between station-specific random effects (i.e., inter-station residuals, black points) and site amplifications (gray lines) obtained through the spectral decomposition for 3 stations (NOR, CLF, NCR), using the results for station BZZ as normalization factor. Bottom. Comparison between inter-event residuals δBe at three frequencies (0.3, 7, 20 Hz) and the Brune stress drop. The inter-station and interevent residuals are taken from Bindi et al. (2019b) whereas station NOR (Norcia Castellina, vs30 = 423 m/s) is installed in the Norcia basin, a 3 by 3 km closed shape basin in central Apennines where amplifications are related to both the presence of deep sediments and to 2D/3D effects Luzi et al. 2018). In all cases, the site amplifications obtained in this study are in good agreement with outcomes from earlier studies that investigated site amplification effects in Central Italy (Luzi et al. 2005). Therefore, the site amplifications obtained through the non-parametric GIT inversion can be used to validate (and interpret) the site specific amplification terms obtained when developing GMPEs. Figure 15 (top panels) compares the GIT site amplifications at three stations (NOR, CLF, and NCR) with the site-specific residuals δS2S associated to the GMPE developed for Fourier amplitudes by Bindi et al. (2019b) to perform a sanity check of ESM-ff. Since the δS2S distribution and the GIT site terms are relative to different reference site conditions, in Fig. 15 both δS2S and the GIT amplifications have been normalized to the result for station BZZ (Bazzano, vs30 = 679 m/s), one of the stations used for defining the GIT reference site condition. The good agreement shown in Fig. 15 confirms that GIT amplifications can be used as empirical site amplifications for site specific (partially non-ergodic) hazard studies, after normalizing both the GMPE median and the empirical site amplification to a common reference level.
The last component of the spectral model corresponds to the source spectra. The possibility to analyze simultaneously the largest events occurred in Europe in the last decades allowed us to get consistent source spectra and, in turn, to compare directly their source parameters. For example, the spectra shown in Fig. 9 for DS1 (Italy) highlight that the M6.4, 1976 Friuli earthquake has the highest spectral content at high frequencies. The source parameters of the Friuli earthquake have been studied by several authors. For example, Cocco and Rovelli (1989) estimated for the Friuli mainshock a (Brune) stress drop of 78 MPa and an average apparent stress (Wyss 1970) of 20 MPa; Castro et al. (1997) estimated a Brune stress drop of 7.14 MPa and an apparent stress drop of 19.3 MPa, and highlighted the complexity of the source spectra which were better described by a double corner frequency model. Rovelli et al. (1991) found an average stress drop of 30 MPa for the aftershocks sequence. The Brune stress drop found in this study for the Friuli mainshock is ∆σ = (37.6 ± 6.5) MPa. The differences among these ∆σ estimates are primarily due to differences in the estimated corner frequencies (f c = 0.47, 0.10, and 0.29 Hz for Cocco and Rovelli 1989;Castro et al. 1997, and this study, respectively) and, to a lesser extent, to different values used for the seismic moments (Mo = 2.9, 4.0 and 5.95 10 18 Nm for Cocco and Rovelli 1989;Castro et al. 1997, and this study, respectively). Despite the large spread in the estimates, the higher stress drop obtained in this study with respect, for example, the Mw 6.5, 2016 Norcia earthquake (∆σ = 6.7 ± 0.8 MPa) is in agreement with the outcomes of previous studies investigating the difference of the reverse fault earthquakes in the Friuli region with respect to the normal faulting earthquakes typical of central Italy (Cocco and Rovelli 1989). For the other large normal faulting earthquakes in central and southern Apennines, we obtained values similar to the value for Norcia (i.e., ∆σ = (6.2 ± 0.9) MPa and ∆σ = (7.3 ± 1.5) MPa, for the 2016, M6.2 Amatrice and the 1980 M6.8 Irpinia earthquakes, respectively) whereas for the 2009, M6.3 L'Aquila earthquake we found a smaller stress drop, i.e. ∆σ = (3.4 ± 0.6) MPa. A factor almost 2 between the stress drop of the Norcia or Amatrice earthquakes and the L'Aquila one is in agreement with the results of a recent regional study focusing on the scaling of source parameters in central Italy, using a large data set of about 1400 earthquakes and analysing both weak and strong motion data . The ∆σ of L'Aquila is similar to the value obtained for the largest shocks of the 1997-1998 Umbria-Marche sequence, which occurred along a segment of the Apennines northern than the area struck by the 2016-2017 sequence (Chiaraluce et al. 2003). For the first (26/09/2016 at 00:32, M5.7), the second (26/09/2016 at 09:40, M6.0) and third (14/10/2016 at 15:23, M5.6) shocks we obtained ∆σ = (3.7 ± 0.5) MPa, (3.6 ± 0.6) MPa, and (3.2 ± 0.5) MPa, respectively. The values for other earthquakes are shown in Table ESM4 of the online Resources. Therefore, ∆σ for normal faulting earthquakes that occurred along the Apennines varies from 3 to 7 MPa, being the largest ones those of the 2016-2017 sequence and of the 1980 Irpinia earthquake.
Regarding DS2, ∆σ of the M6.8, 2014 Aegean earthquake is (15.9 ± 3.2) MPa, whereas ∆σ for the M7.4 Izmit and 7.1 Düzce earthquakes are (37.6 ± 7.2) MPa and (16.0 ± 2.9) MPa, respectively, in agreement with earlier estimates from Edwards and Fäh (2013). Although the absolute values of these estimates could be biased by the usage of a circular rupture model, they suggest a relative high stress drop for the Izmit earthquake and similar stress drops for the 2014 Aegean and 1999 Düzce earthquakes, differently from early suggestions of low stress drop for the 2014 Aegean event (Konca et al. 2018).
The stress drop plays a strong role in driving the level of ground shaking at high frequencies (e.g., Baltay et al. 2019). Considering GMPEs developed for the moment magnitude, which is a static measure of the earthquake size, several studies showed that there is a strong correlation between event-specific residuals computed for high frequency strong motion parameters and the stress drop variability (e.g., Bindi et al. 2007Bindi et al. , 2019aBaltay et al. 2017;Oth et al. 2017;Trugman and Shearer 2017 ). Figure 15 (bottom raw) compares the event specific residuals δBe (inter-event or between event residuals) with ∆σ at three frequencies (0.3, 7, and 20 Hz). The correlation, quantified in terms of Pearson correlation coefficient (cP) (Davis 1986), is strongest at 7 Hz (cP = 0.75), weaker at higher frequencies (cP = 0.56), and no correlation at low frequencies (cP = 0.37). This result agrees with the trend of cP with frequency observed in central Italy (Bindi et al. 2019a) that found a strongest correlation for frequencies close to the corner frequency f c (Eq. 7) whereas the correlation was stronger with k source (Eq. 7) at high frequencies. Therefore, either considering magnitude scales more representative of the earthquake strength and its shaking potential, such as the energy magnitude (Bormann and Di Giacomo 2011), or the rapid assessment of ∆σ could be of relevance for rapid response applications. Recently, Baltay et al. (2019) proposed to determine stable estimate of the stress drop using the Arias intensity (∆σ arias ). In order to extend the distance range of applicability when moderate size events are considered, we derived attenuation models to capture the anelastic attenuation effects and to calibrate station-specific adjustments. The parametric model confirmed the significant magnitude dependency of the attenuation of the Arias intensity to duration ratio. For large magnitudes (M7), the attenuation is closer to the expected 6/5 power of the geometrical spreading term whereas for small magnitudes (M5) the attenuation is much stronger and similar to the spectral attenuation obtained for high frequencies (f = 10 Hz). Contrary to the results obtained by Baltay et al. (2019) for the NGA2 data, we found that ∆σ arias increases with magnitude for M5.5-, whereas it is almost constant for larger events. The calibrated attenuation models confirm that the ratio between Arias intensity and duration for large magnitude events (M7 +) attenuates similarly to the low frequency (< 1 Hz) spectral attenuation whereas, as expected, the effect of anelastic attenuation is stronger for small magnitudes (M5) that attenuate faster, with an attenuation rate close to the spectral attenuation obtained for high frequencies (> 5 Hz). Since a large variability of local site amplification effects are expected at intermediate and high frequencies controlling the Arias intensity, the attenuation models calibrated for computing ∆σ arias include also station corrections in a similar fashion as the station corrections considered when calibrating, for example, a local magnitude scale. Several studies showed that the magnitude station corrections allow to reduce the variance of the distribution of the station magnitudes, leading to more precise event magnitudes (e.g., Baumbach et al. 2003). Also for ∆σ arias , the station corrections turned to be significant, with a standard deviation (for the logarithm of ∆σ arias in base 10) of about 0.4.

Conclusion
In this study, we derived spectral models for source, propagation and site effects for Europe and Middle East by analyzing Fourier amplitude spectra (FAS) of M4 + earthquakes occurred since 1970s. We considered the FAS disseminated along with the flat file extracted from the Engineering Strong Motion data base. We showed that the non-parametric inversion schemes allow to capture regional attributes of the spectral attenuation controlled by heterogeneity at the crustal spatial scale (e.g. changed in the apparent rate of attenuation due to post-critical reflections at strong interfaces) and by differences in the anelastic attenuation. Attenuation models have been derived simultaneously for three different regions, the first two including mainly Italy and the Aegean-western Anatolia regions, and the third the rest of Europe. The attenuation is stronger for Italy, in particular at low-intermediate frequencies.
Once parametrized in terms of geometrical spreading and anelastic attenuation, lower frequency dependent Q values are obtained for the Aegean-western Anatolia region, coupled to a piece-wise linear geometrical spreading showing smaller slopes and the first breakpoint located at shorter distances than for Italy. The non-parametric spectral decomposition provided source spectra for the largest events occurred in Europe in the last decade that is, derived following the same procedure, allowing for the first time to compare them directly. The Brune stress drop varies over about 2 orders of magnitude (the 5th, 50th and 95th percentiles of the ∆σ distribution are 0.76, 2.94, and 13.07 MPa, respectively), with larger events having larger stress drops. In particular, the 5th, 50th and 95th percentiles for M5.5 + are 2.87, 6.02, and 23.5 MPa, respectively; for M5.5-, the same percentiles are 0.73, 2.84, and 12.43 MPa. The comparison of the source spectra confirmed the large spectral content at high frequencies of the M6.4, 1976 Friuli earthquake, occurring over a thrust fault. Once interpreted in terms of Brune model, the correspond a stress drop ∆σ is (37.6 ± 6.5) MPa, larger than the stress drop values obtained for the normal faulting mainshocks occurred in the central and southern sectors of the Apennines (i.e., ∆σ in the range 3.2-7.3 MPa). Finally, following a recent study (Baltay et al. 2019), we considered the estimate of the stress drop from the Arias intensity and we derived attenuation models including station corrections to evaluate the stress drop starting from the Arias intensity and the significant duration. The results confirmed the high potential of this approach for real time and rapid response applications since a rapid evaluation of the stress drop would improve the assessment of the earthquake shaking and damage potential with respect considering the moment magnitude alone.