Near-Source Simulation of Strong Ground Motion in Amatrice Downtown Including Site Effects

: On 24 August 2016, a Mw 6.0 earthquake started a damaging seismic sequence in central Italy. The historical center of Amatrice village reached the XI degree (MCS scale) but the high vulnerability alone could not explain the heavy damage. Unfortunately, at the time of the earthquake only AMT station, 200 m away from the downtown, recorded the mainshock, whereas tens of temporary stations were installed afterwards. We propose a method to simulate the ground motion affecting Amatrice, using the FFT amplitude recorded at AMT, which has been modiﬁed by the standard spectral ratio (SSR) computed at 14 seismic stations in downtown. We tested the procedure by comparing simulations and recordings of two later mainshocks (Mw 5.9 and Mw 6.5), underlining advantages and limits of the technique. The strong motion variability of simulations was related to the proximity of the seismic source, accounted for by the ground motion at AMT, and to the peculiar site effects, described by the transfer function at the sites. The largest ampliﬁcation characterized the stations close to the NE hill edge and produced simulated values of intensity measures clearly above one standard deviation of the GMM expected for Italy, up to 1.6 g for PGA.


Introduction
In recent years, the increased number of permanent and temporary seismic stations have allowed researchers to record strong ground motions close to the seismic source, highlighting the large complexity of the waveforms due to the mixture of source, propagation, and site effects [1][2][3]. The ground motion simulation in the near field is then very difficult to perform, especially when the surface geology and morphology are highly heterogeneous. Several advanced simulation techniques have been implemented so far, to account for complex source and site effects, such as full-wavefield simulations (e.g., [4]), empirical Green's function approaches (useful for reproducing the recorded seismograms in a large frequency band without any knowledge of the underground medium; e.g., [5]), hybrid broadband techniques (e.g., [6]), or dynamic models (e.g., [7]). Despite the increasing level of accuracy, the advanced simulation methods have limitations in frequency range and require an accurate knowledge of the propagation medium and rupture process [5,8]. All these uncertainties affect the synthetic earthquake seismograms, whose variability can be induced by the rupture behavior and errors in propagation [9], and limit the reproduction of the experienced ground motion, such as for understanding peculiar effects on the damage distribution in the near-source area. This is the case of Amatrice village (Italy), whose historic center has been destroyed by a Mw 6.0 earthquake that occurred on 24 August 2016. The so-called Amatrice earthquake started a long seismic sequence in central Italy (Figure 1), causing 299 fatalities and about 30,000 to become homeless. It was followed on 26 October by a Mw 5.9 (Visso earthquake) ∼25 km to the north, and on 30 October by a Mw 6.5 (Norcia earthquake), that nucleated in between the source regions of the two quake a few hundred meters away from the historical center (AMT station).
The applied methodology was tested on the 26 October (Mw 5.9) and 30 October (Mw 6.5) 2016 mainshocks, with the aim to highlight strengths and weaknesses of the approach. We then simulated the Mw 6.0 earthquake of 24 August by using the recordings at AMT and the empirical transfer functions at several sites in downtown. The results have been discussed in terms of time series and intensity measures (PGA, PGV, Arias intensity, significant duration).  [35]); blue, red, and yellow dots represent the location of the seismic events used in this study. (B) Schematic geological map of the Amatrice downtown (modified after [27]) generated by GIS (base map source: Esri, user community, geographic information system. Coordinate  [35]); blue, red, and yellow dots represent the location of the seismic events used in this study. (B) Schematic geological map of the Amatrice downtown (modified after [27]) generated by GIS (base map source: Esri, user community, geographic information system. Coordinate system and projection: World Geodetic System 1984-Web Mercator Auxiliary Sphere); black triangle represents the position of seismic stations whose recordings are used for the analysis.
This short review highlights that there is no single interpretation on the explanation of the ground motion experienced in Amatrice, neither on the causes of the variability of heavy damage [25,29,30,34]. Indeed, the ground motion complexity, due to a combination of nearsource and site effects, is difficult to reproduce with standard simulation techniques [21,29].
To overcome the limitations in terms of frequency band, spatial resolution, and knowledge of source and propagation models, we propose a simple methodology for reproducing the ground motion in Amatrice downtown during the 24 August 2016, Mw 6.0 mainshock, in absence of recording stations in the most damaged area. Empirical transfer functions were used to recover the ground motion that could have hit the downtown during the mainshock, through the convolution with the only record available for that earthquake a few hundred meters away from the historical center (AMT station).
The applied methodology was tested on the 26 October (Mw 5.9) and 30 October (Mw 6.5) 2016 mainshocks, with the aim to highlight strengths and weaknesses of the approach. We then simulated the Mw 6.0 earthquake of 24 August by using the recordings at AMT and the empirical transfer functions at several sites in downtown. The results have been discussed in terms of time series and intensity measures (PGA, PGV, Arias intensity, significant duration).

Available Data
In this study, we used seismic signals of 35 earthquakes of the 2016-2017 central Italy seismic sequence (Table 1 and Figure 1A), recorded by 20 seismic stations located in the Amatrice village (Table 2 and Figure 1B).  Table 2. List of the seismic stations installed in the Amatrice municipality (see [25,31] The selected events were characterized by epicentral distances between 7 to 48 km from Amatrice downtown and magnitudes in the range ML ≥ 1.5 and Mw ≤ 6.5, including the three mainshocks of the sequence (Mw 6.0 on 24 August, Mw 5.9 on 26 October, and Mw 6.5 on 30 October 2016). The recordings were downloaded from the European Integrated Data Archive (EIDA, http://eida.ingv.it, accessed on 24 April 2021; [36]) and the Italian Strong Motion Network (RAN, http://ran.protezionecivile.it, accessed on 24 April 2021; [37]) databases. They were selected accounting for a signal-to-noise spectral ratio greater than two or three (depending on the stations) in a wide frequency band, and a clear P and S wave arrivals.
The closest seismic station to Amatrice downtown that recorded the 24 August, Mw 6.0, earthquake was AMT, belonging to the National Accelerometric Network [37] and able to record seismic events with a Ml ≥ 2.5 [17]. AMT was at an epicentral distance of about 8.5 km (Rjb = 1km; [38]) and recorded the maximum horizontal PGA for the event on the EW component (8.5 m/s 2 ; [2]); the difference with the NS component, that recorded half of the PGA value, was probably caused by topographic effects [29], seismic amplification [25], or directivity effects [14].
Another reference station, T1299, was installed by the INGV mobile seismic network about 800 m NW of Amatrice downtown, at the base of the hill [25,39].
Moreover, to record the aftershock sequence and to investigate the local amplification effects for microzonation purposes, a total of 50 temporary seismic stations were installed by INGV in four municipalities, including Amatrice (MZ stations, belonging to the 3A network [40]; see also [25,[31][32][33]). In this study, only six stations distributed along the NW-SE elongated Amatrice terrace were considered: the stations MZ10, MZ12, MZ28 were installed at the top of the sandy conglomeratic deposits of the Amatrice-Sommati units [27], whereas the stations MZ29, MZ30 were located in a NE-SW direction across the Amatrice hill.
According to [25], the station MZ08 was installed a few meters far from the AMT station, to extend the detection capability at the site towards the lower magnitude events. For events before 4 November 2016, we verified a 13 • clockwise orientation difference of AMT with respect to MZ08. However, after that date, the sensor at AMT was changed and the orientation was corrected. The slight difference in the sensor's orientation does not significantly affect the seismic recordings, showing very similar signals (see Figure S1 on Supplementary Materials).
Later on, an array of 24 seismic stations was installed for few hours of activity in the heavily damaged historical center [25]. Out of 12 stations (CS01 to CS12) that worked during 28 June 2017, only four stations (CS01 to CS04) recorded a clear seismic signal of two earthquakes with ML 1.8 and 2.2. Otherwise, 12 stations (CS13 to CS24) worked during the second day of activity, but only seven of them (CS13, CS14, CS16, CS17, CS20, CS23, CS24) recorded four seismic events with a ML between 1.5 to 3.1.

Simulation Method
During the Amatrice earthquake, no station was available in the historical center except for AMT. Later on, tens of temporary stations were installed in Amatrice village and they recorded many earthquakes simultaneously to the reference station.
To overcome the limits of the available data and the network geometry, we implemented a strategy for simulating a given earthquake (e.g., the Amatrice mainshock) at stations (e.g., Table 2) close to the reference one (e.g., AMT) that recorded the event. The strategy is summarized in Figure 2 and explained in the following steps:

1.
The target event is recorded at a reference station, for which we calculated the Fourier transform (FFT phase and amplitude).

2.
Several earthquakes are recorded simultaneously at the reference station and at the other neighboring sites (in this way we can assume similar source effects at all stations), allowing us to calculate the empirical transfer functions respect to the reference site by means of the standard spectral ratio (SSR). 3.
FFT amplitude of the target event recorded at the reference station is multiplied by the SSR of each station; the result is an amplitude spectral content, modified by the contribution of a transfer function.

4.
The modified FFT amplitude spectrum is back-transformed to the time domain through the IFFT, using the FFT phase of the reference station. However, we are aware that the use of the same phase spectrum at all stations results in an unlikely coherence at high frequencies between the signals recorded at different sites, whereas the high frequency incoherency is related to the heterogeneities of the propagation path and it is largely site dependent, but also regionally dependent [41][42][43]. Considering that the spatial coherency decreases with increasing frequency and distance between measuring points [44], the same authors showed that the coherence for high-frequency intensity measures (such as PGA) is still high within a separation distance of one to a few kilometers, which makes the stations configuration in Amatrice reliable for our methodology.
Another possible weakness, due to the use of the same phase spectrum, relates to the inability to capture the increase of the ground-motion duration for 2-3D site effects, which can act to decrease the peak amplitudes.
Finally, the small-magnitude events for the SSR computation do not account for possible nonlinear soil behavior, whose main effect would be a shift of the amplification towards lower frequencies and a decrease of the peak amplitudes [45].  The proposed procedure overcomes the complexity of the source and site effects: the near-source influence is included in the strong motion recording at the reference site, whereas the relative site effects are accounted for by the site transfer functions.
However, we are aware that the use of the same phase spectrum at all stations results in an unlikely coherence at high frequencies between the signals recorded at different sites, whereas the high frequency incoherency is related to the heterogeneities of the propagation path and it is largely site dependent, but also regionally dependent [41][42][43]. Considering that the spatial coherency decreases with increasing frequency and distance between measuring points [44], the same authors showed that the coherence for highfrequency intensity measures (such as PGA) is still high within a separation distance of one to a few kilometers, which makes the stations configuration in Amatrice reliable for our methodology.
Another possible weakness, due to the use of the same phase spectrum, relates to the inability to capture the increase of the ground-motion duration for 2-3D site effects, which can act to decrease the peak amplitudes.
Finally, the small-magnitude events for the SSR computation do not account for possible nonlinear soil behavior, whose main effect would be a shift of the amplification towards lower frequencies and a decrease of the peak amplitudes [45].

Data Processing
The seismic recordings used for the SSR computation are listed in Table 1, together with the stations that recorded each event.
The time series have been band-pass filtered between 0.5 to 20 Hz, corrected by subtracting the mean value and the best-fit line, cut between −0.5 to 30 s with respect to the P-wave arrival (manual picking), and tapered before computing the Fourier spectra. Only the frequencies with a signal-to-noise ratio (s/n) ≥ 3 were used for the recordings at T1299 and MZ stations and (s/n) ≥ 2 for the CS stations, to maximize the information of the few low-magnitude earthquakes available.
As described in [46], the SSR were computed on 18 horizontal components, each rotated every 10 • clockwise, and on the vertical component. After rotation of the time signals, their Fourier spectra were smoothed with Konno-Ohmachi algorithm [47] (b = 40 and fmax = 20 Hz) and divided by the reference ones. The spectral ratios for each station were then averaged (geometrical mean) on all the used seismic events (Figure 3).

Data Processing
The seismic recordings used for the SSR computation are listed in Table 1, together with the stations that recorded each event.
The time series have been band-pass filtered between 0.5 to 20 Hz, corrected by subtracting the mean value and the best-fit line, cut between −0.5 to 30 s with respect to the P-wave arrival (manual picking), and tapered before computing the Fourier spectra. Only the frequencies with a signal-to-noise ratio (s/n) ≥ 3 were used for the recordings at T1299 and MZ stations and (s/n) ≥ 2 for the CS stations, to maximize the information of the few low-magnitude earthquakes available.
As described in [46], the SSR were computed on 18 horizontal components, each rotated every 10° clockwise, and on the vertical component. After rotation of the time signals, their Fourier spectra were smoothed with Konno-Ohmachi algorithm [47] (b = 40 and fmax = 20 Hz) and divided by the reference ones. The spectral ratios for each station were then averaged (geometrical mean) on all the used seismic events (Figure 3).
We computed the SSR of the T1299, MZ, and CS stations with respect to AMT or MZ08 (they are equivalent), using different earthquakes (Table 1 and Figure 1)-for T1299 we used 15 earthquakes with magnitude between 2.6 and 5.4, recorded simultaneously at T1299 and AMT, whereas for the MZ stations we considered 12 events of magnitude range 3.1-5.9, using MZ08 as reference ( Figure S2 in Supplementary Materials).
Regarding the CS stations in downtown, the CS13-14-16-17-20-23-24 recorded four seismic events simultaneously to AMT, but only two, with magnitude 2.6 and 3.1, were recorded by AMT and showed a signal-to-noise ratio greater than 2 in a wide frequency band ( Figure 4 and Figure S3 in Supplementary Materials). Unfortunately, we could not use CS01-02-03-04 stations because they recorded two seismic events only, having too low magnitude to be recorded by the accelerometric station AMT (Table 1).  We computed the SSR of the T1299, MZ, and CS stations with respect to AMT or MZ08 (they are equivalent), using different earthquakes (Table 1 and Figure 1)-for T1299 we used 15 earthquakes with magnitude between 2.6 and 5.4, recorded simultaneously at T1299 and AMT, whereas for the MZ stations we considered 12 events of magnitude range 3.1-5.9, using MZ08 as reference ( Figure S2 in Supplementary Materials).
Regarding the CS stations in downtown, the CS13-14-16-17-20-23-24 recorded four seismic events simultaneously to AMT, but only two, with magnitude 2.6 and 3.1, were recorded by AMT and showed a signal-to-noise ratio greater than 2 in a wide frequency band ( Figure 4 and Figure S3 in Supplementary Materials). Unfortunately, we could not use CS01-02-03-04 stations because they recorded two seismic events only, having too low magnitude to be recorded by the accelerometric station AMT (Table 1).  To quantify the intensity measures of the seismic signals, several parameters have been taken into account: the peak ground acceleration (PGA) and velocity (PGV), providing limited insight to the shaking at high and intermediate frequencies, respectively; the Arias intensity; the significant duration. The Arias intensity (AI) represents the cumulative energy perceived at the site, in a specific time interval during a seismic event, and it is expressed in m/s [48]; we computed it for the total duration of the seismic signal (from t = 0 to t = tmax): The significant duration (SID) uses Equation (1) and measures the time interval (t2-t1), in which Ia is between 5% and 95% of the total.
For the intensity measures on horizontal motion, we used the maximum between the NS and EW components. See [46] for a detailed description of the signal processing procedure. To quantify the intensity measures of the seismic signals, several parameters have been taken into account: the peak ground acceleration (PGA) and velocity (PGV), providing limited insight to the shaking at high and intermediate frequencies, respectively; the Arias intensity; the significant duration. The Arias intensity (AI) represents the cumulative energy perceived at the site, in a specific time interval during a seismic event, and it is expressed in m/s [48]; we computed it for the total duration of the seismic signal (from t = 0 to t = t max ): The significant duration (SID) uses Equation (1) and measures the time interval (t 2 -t 1 ), in which Ia is between 5% and 95% of the total.
For the intensity measures on horizontal motion, we used the maximum between the NS and EW components. See [46] for a detailed description of the signal processing procedure.

Results
The proposed simulation procedure has been first tested by comparing simulated and recorded seismograms for two mainshocks (Mw 5.9 and 6.5 on 26 October and 30 October 2016, respectively) at some stations in Amatrice, using AMT and MZ08 as a reference site. Then, we simulated the expected ground motion in Amatrice downtown during the Mw 6.0 earthquake on 24 August, by means of the recordings at AMT station.

Standard Spectral Ratio (SSR)
The site effects of the Amatrice stations have been already discussed in [25,33], although their reference station was T1299; however, we list here some main findings inferred from their and our analyses (see Figures 3 and 4 and Figure S2 in Supplementary Materials): • T1299 is affected by a slight deamplification with respect to AMT up to 10 Hz; above this value, the amplification increases by up to 2.5 ( Figure 3). and Figure S2 in Supplementary Materials).

Tests
The simulation strategy has been tested to investigate the contribution of the transfer functions and the effect of the reference station phase used in the IFFT ( Figure 5). We evaluated PGV, PGA, Arias intensity (AI) and signal duration (SID) of the simulated and recorded seismic signals of Mw 5.9, 26 October (Test 1 to 3), and Mw 6.5, 30 October (Test 4), at MZ08-10-12-28 stations (Figure 6).
PGA and PGV were also compared with the GMM valid for shallow crustal earthquakes in Italy (ITA18; [49]); the recorded values of the Mw 5.9 earthquake at MZ08 and MZ10 stations were close to the GMM predictions plus 1 standard deviation, while MZ12 and MZ28 showed values larger than those predicted. For the Mw 6.5 earthquake (test 4), the recordings overestimated the GMM values by a larger extent, especially for MZ12 and MZ28, where recorded values almost doubled the prediction plus 1 standard deviation.
To estimate the goodness of fit of the simulations with respect to the observations (GoF), we computed the relative difference as follow: where IM sim is the intensity measure from the simulated time series and IM rec is computed from real records. The obtained values of GoF for the 4 tests can be found on Table S1 (Supplementary Materials). The first test (TEST 1 in Figure 5) was carried out to verify the formal correctness of the procedure, whereas the second test (TEST 2 in Figure 5) aimed to investigate the use of the FFT phase of reference station (MZ08) in the IFFT process. The recorded intensity measures were well reproduced for Test 2 ( Figure 6), except for MZ28 where the simulations underestimated the recorded PGV and PGA (GoF equal to −34% and −26%, respectively). As expected, the use of the reference phase in the IFFT caused a time shift of the simulated signals with respect to the recorded ones, proportional to the distance between the considered station and the reference site ( Figure 7). The first test (TEST 1 in Figure 5) was carried out to verify the formal correctness of the procedure, whereas the second test (TEST 2 in Figure 5) aimed to investigate the use of the FFT phase of reference station (MZ08) in the IFFT process. The recorded intensity measures were well reproduced for Test 2 ( Figure 6), except for MZ28 where the simulations underestimated the recorded PGV and PGA (GoF equal to −34% and −26%, respectively). As expected, the use of the reference phase in the IFFT caused a time shift of the simulated signals with respect to the recorded ones, proportional to the distance between the considered station and the reference site ( Figure 7).
The third test (TEST 3 in Figure 5) was performed to investigate the use of the SSR average on 12 earthquakes (Table 1 and Figure S2, Supplementary Materials). As for Test 2, the fit was very good for all the indicators and it got better for the PGA; MZ28 was again underestimated by an amount of −42% for PGV and −22% for PGA (Figures 6 and 8).  The third test (TEST 3 in Figure 5) was performed to investigate the use of the SSR average on 12 earthquakes (Table 1 and Figure S2, Supplementary Materials). As for Test 2, the fit was very good for all the indicators and it got better for the PGA; MZ28 was again underestimated by an amount of −42% for PGV and −22% for PGA (Figures 6 and 8).  The third test (TEST 3 in Figure 5) was performed to investigate the use of the SSR average on 12 earthquakes (Table 1 and Figure S2, Supplementary Materials). As for Test 2, the fit was very good for all the indicators and it got better for the PGA; MZ28 was again underestimated by an amount of −42% for PGV and −22% for PGA (Figures 6 and 8). The fourth test (TEST 4 in Figure 5) mimicked Test 3 but for the Mw 6.5 earthquake, to investigate the use of small-magnitude SSR in case of possible nonlinear effects for a larger magnitude event. Figure 9 shows the different ingredients of the simulation strategy applied to MZ12: recordings of the reference station MZ08 ( Figure 9A); recorded and simulated accelerations at MZ12, showing a slight amplification of the simulated signal with respect to the recorded one ( Figure 9B); amplitude Fourier spectra of recorded signals at both MZ08 and MZ12, compared with the simulated FFT at MZ12 ( Figure 9C); SSR used to correct the amplitude FFT spectrum of MZ08 ( Figure 9D).
The comparison of the intensity measures ( Figure 6) showed that the PGV at all stations were well fitted, slightly overestimating the recorded values at MZ10 and MZ28 of The fourth test (TEST 4 in Figure 5) mimicked Test 3 but for the Mw 6.5 earthquake, to investigate the use of small-magnitude SSR in case of possible nonlinear effects for a larger magnitude event. Figure 9 shows the different ingredients of the simulation strategy applied to MZ12: recordings of the reference station MZ08 ( Figure 9A); recorded and simulated accelerations at MZ12, showing a slight amplification of the simulated signal with respect to the recorded one ( Figure 9B); amplitude Fourier spectra of recorded signals at both MZ08 and MZ12, compared with the simulated FFT at MZ12 ( Figure 9C); SSR used to correct the amplitude FFT spectrum of MZ08 ( Figure 9D).
The comparison of the intensity measures ( Figure 6) showed that the PGV at all stations were well fitted, slightly overestimating the recorded values at MZ10 and MZ28 of about 10% and 20%, respectively. The PGA values were well reproduced at MZ12, but we overestimated the recorded values of about 47% and 32% at MZ10 and MZ28. The Arias intensity values were well reproduced at MZ10 and MZ12, with a GoF of about ±15%, while at MZ28 we overestimated the recorded values of about 65%. The signal duration (SID) showed an increase of about 35% at MZ12, and ±5% at MZ10 and MZ28.

24 August 2016 Simulation
The Mw 6.0 earthquake on 24 August 2016, was simulated at 14 seismic stations, located mainly in downtown Amatrice, using AMT as a reference site (Table 2, Figure 1). The large number of seismic events recorded by T1299 and MZ and the two events at CS13 to CS24 allowed us to reproduce reliable seismic signals in the 0.5-20 Hz frequency band ( Figure 10): the highest ground shaking was concentrated in the first seconds after the P-wave arrival and the signal duration was similar at all the stations (4 to 6 s in Figure 11D).
Moreover, there was a large variability of the peak values at the different station locations ( Figure 11 for the maximum of horizontal components, Figure S4 for the vertical component), and almost all the stations had PGA and PGV values up to three times the mean plus one standard deviation of the GMM valid for the area (Figure 11). During the Mw 6.0 shaking, the AMT station recorded a PGV value of 0.3 m/s (maximum on the horizontal components), slightly above one standard deviation of the GMM expected value [49]. With the exception of four stations, whose simulated PGV values were similar to AMT (GoF = −15% to −21% at T1299-MZ29-MZ31, GoF = +14% at MZ10), the other stations showed larger values with GoF ranging mostly from 50% to 100%, and with a maximum of 150% or 200% for CS20 and CS16, respectively ( Figure 11A). about 10% and 20%, respectively. The PGA values were well reproduced at MZ12, but we overestimated the recorded values of about 47% and 32% at MZ10 and MZ28. The Arias intensity values were well reproduced at MZ10 and MZ12, with a GoF of about ±15%, while at MZ28 we overestimated the recorded values of about 65%. The signal duration (SID) showed an increase of about 35% at MZ12, and ±5% at MZ10 and MZ28.

24 August 2016 Simulation
The Mw 6.0 earthquake on 24 August 2016, was simulated at 14 seismic stations, located mainly in downtown Amatrice, using AMT as a reference site (Table 2, Figure 1). The large number of seismic events recorded by T1299 and MZ and the two events at CS13 to CS24 allowed us to reproduce reliable seismic signals in the 0.5-20 Hz frequency band ( Figure 10): the highest ground shaking was concentrated in the first seconds after the Pwave arrival and the signal duration was similar at all the stations (4 to 6 s in Figure 11D). Regarding PGA, AMT recorded a value of about 8 m/s 2 , largely above one standard deviation of the GMM ( Figure 11B). Again, T1299-MZ31 had smaller values (GoF about −50%), whereas CS14-CS20-MZ12-MZ28-CS16 exceeded the AMT values by more than 40% (GoF = 100% at CS16). The other stations varied between −20% to 20% of AMT values.
The Arias intensity (AI), shown in Figure 11C, was between 100% and 200% larger than AMT, with two stations reaching 350% and 600% (CS20 and CS16, respectively). T1299 and MZ31 decreased by more than 60%, whereas MZ10 and MZ29 were very similar to AMT. The significant durations (SID, Figure 11D) of the simulated signals were similar to the recorded one at AMT, characterized by a total duration of about 4 s.
Finally, the intensity measures for the vertical component ( Figure S4 in Supplementary Materials) were always larger than AMT with the exception of T1299 and MZ31, and they never exceeded the 100% of AMT for PGA and PGV. Arias intensity values, instead, were larger than AMT by about 100-350%. The significant duration of the simulated signals showed a slight variation, with a maximum increase of about 1.1 s with respect to AMT. Moreover, there was a large variability of the peak values at the different station locations ( Figure 11 for the maximum of horizontal components, Figure S4 for the vertical component), and almost all the stations had PGA and PGV values up to three times the mean plus one standard deviation of the GMM valid for the area (Figure 11). During the Mw 6.0 shaking, the AMT station recorded a PGV value of 0.3 m/s (maximum on the horizontal components), slightly above one standard deviation of the GMM expected value [49]. With the exception of four stations, whose simulated PGV values were similar to AMT (GoF = −15% to −21% at T1299-MZ29-MZ31, GoF = +14% at MZ10), the other stations showed larger values with GoF ranging mostly from 50% to 100%, and with a maximum of 150% or 200% for CS20 and CS16, respectively ( Figure 11A).
The Arias intensity (AI), shown in Figure 11C, was between 100% and 200% larger than AMT, with two stations reaching 350% and 600% (CS20 and CS16, respectively). T1299 and MZ31 decreased by more than 60%, whereas MZ10 and MZ29 were very similar to AMT. The significant durations (SID, Figure 11D) of the simulated signals were similar to the recorded one at AMT, characterized by a total duration of about 4 s.
Finally, the intensity measures for the vertical component ( Figure S4 in Supplementary Materials) were always larger than AMT with the exception of T1299 and MZ31, and they never exceeded the 100% of AMT for PGA and PGV. Arias intensity values, instead,   [49] for Rjb = 1 km [38] and site class B (Vs30 = 670 m/s at AMT).

Discussion
We proposed a procedure that has been tested on recorded earthquakes to evaluate the strengths and the weakness of the implemented method, checking for the contribution of the phase (Test 2) and of the average transfer functions on small earthquakes to reproduce strong shaking (Tests 3-4). The tests showed a fairly good agreement between simulations and recordings, both in terms of time series and intensity measures.

Discussion
We proposed a procedure that has been tested on recorded earthquakes to evaluate the strengths and the weakness of the implemented method, checking for the contribution of the phase (Test 2) and of the average transfer functions on small earthquakes to reproduce strong shaking (Tests 3-4). The tests showed a fairly good agreement between simulations and recordings, both in terms of time series and intensity measures.
However, there are some limits due to the use of the same phase spectrum for all the stations and to the small-magnitude transfer functions. First of all, the reference site phase in the IFFT procedure implied the same P-wave arrival times at all the stations (Figures 7 and 8) and it did not allow us to simulate the variation of the seismic duration due to local site response ( Figure 6D). Secondly, the use of SSR produced signals whose intensity measures were lower than the recordings (Test 2 versus Test 3 in Figure 6) in the case of spectral amplification of the target event larger than the average SSR ( Figure S5 in Supplementary Materials). Instead, when the transfer function is lower than the average SSR, such as for the strongest mainshock (Test 4; Figure 9 and Figure S5 in Supplementary Materials), we overestimated both the total energy AI and the PGV (Figure 6), probably because we were not accounting for possible nonlinear effects of the subsoil during the strong shaking. As a general conclusion, we can assess that our procedure evaluates an upper bound of the expected intensity measures at the target sites.
The proposed methodology has been applied to simulate the 24 August, Mw 6.0 earthquake at sites in the town of Amatrice: seven CS stations located in the historical center and seven farther sites (T1299 and six MZ stations) located on the top and at the base of the Amatrice terrace (Figures 10 and 11). The largest amplification characterized the stations at the proximity of the NE hill edge and tended to decrease in the SW direction toward the central part of downtown. The lack of lateral confinement, due to the morphological escarpment, produced a NE oriented amplification of the low frequencies at the CS stations, on the direction orthogonal to the morphological elongation of the Amatrice terrace. Instead, the amplification at frequencies above 8 Hz were variable in amplitude and preferential orientation, suggesting that the top of the hill is affected by highly variable seismic amplification, due to the presence of low-velocity deposits.
The strong accelerations reproduced at the historical center of Amatrice is then related to the vicinity of the seismic source, accounted for by the ground motion at AMT, and by the local site effects, described by the SSR transfer function at the sites. The presence of layers with different impedance contrasts, together with the steep slopes bounding the Amatrice terrace, is able to generate seismic amplification effects, resulting in the increase of the intensity measures simulated on the top of the hill. The largest motion was simulated at CS16 and CS20, which showed horizontal PGA clearly exceeding the gravitational acceleration (about 16 and 12 m/s 2 , respectively). Similarly, the stations MZ10-12-28-30-31 were characterized by seismic amplification with respect to the reference site: MZ12 and MZ28 showed horizontal PGA of about 14 m/s 2 , whereas MZ10-29-30 had values between 7 to 9 m/s 2 . The two stations at the base of the terrace, MZ31 and T1299, were deamplified with respect to AMT and showed lower intensity measures, highlighting the high variability and the seismic behavior between the base and the top of the Amatrice hill.
Furthermore, the 24 August simulated signals showed a slight variation of the signal duration with respect to AMT (SID in Figure 11). This effect was due to the use of the reference phase spectra in the IFFT procedure, and to the near-source conditions, which concentrated the energy radiation in the direct wave arrivals. Conversely, in case of far-field conditions, as for the low-magnitude earthquakes recorded at CS stations, we verified that the seismic amplification in the Amatrice historical center was associated to an increase of SID up to three times with respect to AMT station (see also [25]).
The PGA and PGV values recorded by AMT, and the simulated ones at the stations, were clearly above one standard deviation of the GMM predictions (Figure 11), suggesting that the Amatrice historical center, especially in the eastern part of downtown, has been subjected to a severe ground shaking larger than the expected average. As already discussed, the proposed procedure did not allow us to simulate the nonlinear effects and it concentrated the energy in a time interval controlled by the reference site, suggesting that we estimated the upper limit of intensity measures in case of a seismic event characterized by similar magnitude and epicentral distance.

Conclusions
The proposed procedure is suitable for simulating strong ground motions of a past earthquakes recorded by at least one reference station close by. In this way, it is possible to overcome the limits of other simulation techniques for which the computation needs the description of the rupture process and propagation model from the source to the site. This is the case of the 24 August 2016, Mw 6.0, earthquake signals which have been reproduced at the historical center of Amatrice, by using the transfer functions with respect to the closest seismic station (AMT) that recorded the event. We simulated the strong ground motion at 11 seismic stations located on the top of the Amatrice hill, out of which seven were installed in the historical center. The largest motion was reproduced at the north-eastern edge of the downtown, where the horizontal PGA values exceeded the gravitational acceleration. The simulated ground motions decreased away from the edges of the hill, suggesting that the topographic effect and the lack of lateral confinement are the predominant factors of the site effect. Two stations located at the base of the terrace showed lower intensity measures with respect to the reference site, highlighting the high variability and the seismic behavior between the base and the top of the Amatrice hill.
The strong ground motion variability during the Mw 6.0 event could help to understand the damage distribution that affected Amatrice downtown (X-XI of MCS; [50]). It is worth noting that the buildings of the historical center belonged to the most vulnerable classes [11,12] and were moderately to completely destroyed by the earthquake. The most damaged buildings were in the eastern part of the town [13], in proximity to the steep slopes of the hill, which was affected by the largest shaking.
Further studies could assess the performance of the simulation strategy for largemagnitude earthquakes having different focal mechanisms and rupture models. Moreover, it should be tested on sites with different geological and morphological conditions (such as soft sediments, fractured rock, or different topographies) and using more than one reference station close by, in order to make a robust evaluation of the simulated intensity measures.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/geosciences11050186/s1, Figure S1: Comparison of the AMT recorded time series of Mw 6.5 earthquake, NS-EW vs. 13 • clockwise rotated components. Figure S2: Comparison between the time series and the FFT amplitude recorded by AMT and MZ08 for the Mw 6.5 earthquake (NS, EW, and vertical components). Figure S3: Standard spectral ratio (SSR) at the seismic stations considered in the study. Figure S4: Recorded and simulated intensity measures on the vertical component of the 24 August 2016, earthquake. Figure S5: comparison between the SSR of MZ10-12-28 averaged on 12 seismic events and the SSR on Mw 5.9 and Mw 6.5 earthquakes individually.
Author Contributions: Conceptualization, G.C.; writing-original draft preparation, A.T. and G.C.; writing-review and editing, A.T. and G.C.; supervision, G.C. All authors have read and agreed to the published version of the manuscript.