Yamanaka, H., Dravinski, M., and Kagami, H. (1993). Continuous Measurements of Microtremors on Sediments and Basement in Los Angeles, California. *Bull. Seismological Soc. America* 83 (5), 1595–1609. doi:10.1785/BSSA0830051595

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

**Publisher's Note:** All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of

the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Copyright © 2022 Janusz, Perron, Knellwolf and Fäh. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
<figure><img src="image_54.png" /></figure>

Check for
updates

# Testing Site Amplification Curves in Hybrid Broadband Ground Motion Simulations of M6.0, 24 August 2016 Amatrice Earthquake, Italy

Marta Pischiutta¹*, Aybige Akinci¹, Chiara Felicetta², Francesca Pacor² and Paola Morasca²

¹Istituto Nazionale di Geofisica e Vulcanologia, Rome, Italy, ²Istituto Nazionale di Geofisica e Vulcanologia, Milano, Italy

This research focuses on predicting and assessing earthquake impact due to future scenarios regarding the ground motion seismic hazard by accounting mainly for site effect in the Central Apennines. To this end, we produced synthetic broadband seismograms by adopting a hybrid simulation technique for the M<sub>w</sub>6.0 Amatrice earthquake, Central Italy, on 24 August 2016, accounting for site conditions by means of amplification curves, computed with different approaches. Simulations were validated by comparing with data recorded at 57 strong-motion stations, the majority installed in urban areas. This station sample was selected among stations recording the Amatrice earthquake within an epicentral distance of 150 km and potentially prone to experience site amplification effects because of lying in particular site conditions (sedimentary basins, topographic irregularities, and fault zones). The evaluation of amplification curves best suited to describe local effects is of great importance because many towns and villages in central Italy are built in very different geomorphological conditions, from valleys and sedimentary basins to topographies. In order to well reproduce observed ground motions, we accounted for the site amplification effect by testing various generic and empirical amplification curves such as horizontal-to-vertical spectral ratios (calculated from Fourier spectra using both earthquake, HVSR, and ambient noise, HVNSR, recordings) and those derived from the generalized inversion technique (GIT). The site amplifications emanated from GIT improve the match between observed and simulated data, especially in the case of stations installed in sedimentary basins, where the empirical amplification curve effectively reproduces spectral peaks. On the contrary, the worst performances are for the spectral ratios between components, even compared to the generic site amplification, although the latter ignores the strong bedrock/soil seismic impedance contrasts. At sites on topography, we did not observe any systematic behavior, the use of empirical curves ameliorating the fit only in a small percentage of cases. These results may provide a valuable framework for developing ground motion models for earthquake seismic hazard assessment and risk mitigation, especially in urban areas located in the seismically active central Italy region.

**Keywords:** site amplification effects, Amatrice earthquake, central Italy, ground motion simulations, horizontal-to-vertical spectral ratio (HVSR), generalized inversion technique (GIT)

OPEN ACCESS

**Edited by:**
Nicola Alessandro Pino,
National Institute of Geophysics and
Volcanology (INGV), Italy

**Reviewed by:**
Chuanbin Zhu,
GFZ German Research Centre for
Geosciences, Germany
Yefei Ren,
China Earthquake Administration,
China

*Correspondence:
Marta Pischiutta
marta.pischiutta@ingv.it

**Specialty section:**
This article was submitted to Solid Earth Geophysics, a section of the journal Frontiers in Earth Science

**Received:** 28 February 2022

**Accepted:** 26 April 2022

**Published:** 13 June 2022

**Citation:**
Pischiutta M, Akinci A, Felicetta C, Pacor F and Morasca P (2022) Testing Site Amplification Curves in Hybrid Broadband Ground Motion Simulations of M6.0, 24 August 2016 Amatrice Earthquake, Italy. Front. Earth Sci. 10:886606. doi: 10.3389/feart.2022.886606
# 1 INTRODUCTION

On 24 August 2016, at 03:36 local time, an Mw 6.0 earthquake with a shallow focal depth (8.0 km) occurred close to the Accumoli village and Amarice town in Central Apennines. This event triggered an extended sequence with five Mw ≥ 5 earthquakes, including a strong shock Mw6.5 on 30 October 2016. The first shock caused heavy damage in several villages, mainly in Accumoli, Amatrice, and Arquata del Tronto, with X to XI MCS intensity values (Galli et al., 2016a; Galli et al., 2016b; Quest, 2016; Zanini et al., 2016), and several ancient building collapsed due to the vicinity of the causative fault and the high vulnerability. During the 2016 Central Italy seismic sequence in various municipalities and hamlets, the damage patterns indicated strong evidence of local site effects (Sextos et al., 2018), mainly related to stratigraphic and topographic effects.

It is fundamental to know the features of the ground shaking during an earthquake to support the interventions and actions both in the emergency and the reconstruction phases. For example, shake maps (Wald et al., 1999; Michelini et al., 2008; Licia Faenza et al., 2016), generated in a quasi-real-time, interpolating observed and predicted data represent the distribution of ground-motion parameters following an earthquake. Generating these maps for future events from a given seismic source's selected locations, magnitudes, and rupture mechanisms has important implications for land use planning and seismic risk mitigation of a given area.

There exist a variety of empirical and numerical methods for generating shaking maps from empirical ground motion models to physics-based approaches (Douglas and Aochi, 2008; Goulet et al., 2015) that implement different strategies to include the local site effects. In empirical models, the site-local effects are introduced through site proxies, among which the S-wave velocity in the first 30 m, $V_{s30}$, is the most common. 1D and 3D site effects may directly be introduced in the numerical physics-based approach, but they need detailed knowledge of the site-local geological and geomorphology setting, including geophysical and geotechnical properties. An alternative strategy is to simulate ground motion at rock and then add the amplification curves, empirically or numerically estimated, using 1D simplified models overlaying rigid substrate. The aim of the work is the inclusion of the site response in shaking scenarios calculations using broadband ground motion hybrid modeling for the 24 August 2016 earthquake, to make available a tool useful to reduce seismic hazards and improve risk mitigations in urban areas. The evaluation of amplification curves best suited to describe local effects is of great importance because many towns and villages in Central Italy are built in very different geomorphological conditions, from valleys and sedimentary basins to topographies.

In sedimentary basins, the presence of superficial soft sediments and strong shear-wave velocity and impedance contrasts led to strong amplification of seismic waves, even produced by earthquakes originating at relevant distances (hundreds of meters). First observations date back to the 80s and 90s of the last century all over the world, firstly involving soft soil deposits: in Mexico city during Michoacán earthquakes (1985, M = 8.1 and 7.5, e.g., Sánchez-Sesma et al., 1988); in

Los Angeles basin during Mw 6.7, 1994 Northridge earthquake (e.g., Graves, 1995); in Osaka basin after M 7.2 1995 Kobe earthquake (e.g., Iwata et al., 1996). The physical mechanism at the basis of the phenomenon involves refraction of seismic waves by a velocity contrast between superficial soft sediments and an underlying stiff bedrock and subsequent phase constructive interference causing a resonance effect. Stratigraphic resonance effects are considered in seismic design codes of many countries for seismic risk mitigation (e.g., Eurocode8 in EU, NTC18 in Italy, NEHRP in the United States, NZS1170.5), through the use of scaling factors defined on the basis of the shear-wave velocity profile and the $V_{s30}$ parameter. The Italian seismic design prescribes five classes: A (average $V_{s30}$ over 800 m/s); B ($V_{s30}$ between 360 and 800 m/s); C ($V_{s30}$ between 180 and 360 m/s); D ($V_{s30}$ lower than 180 m/s); E (particular cases). The former represents rock sites that are considered to be unaffected by site amplification, apart from high-frequency effects due to superficial weathering. Nevertheless, many recent studies have highlighted that even at frequencies of engineering interest (0.5–20 Hz) at rock sites, seismic waves can be amplified due to the local properties of the rock (i.e., the presence of pervasive fractures and/or large open cracks in different domains — fault zones, landslides, volcanoes, for example, Pischiutta et al., 2012, 2017; Panzera et al., 2014; Falsaperla et al., 2010; Ben-Zion and Sammis 2003; Lewis and Ben-Zion, 2010, Felicetta et al., 2018; Lanzano et al., 2020 (to cite a few among many). Maximum amplification (with an increase of over 100%) occurs along a site-dependent azimuth at a high angle to the fault strike; this is the reason for calling such effect “directional amplification” (using a term coined by Bonamassa and Vidale, 1991).

Finally, even sites on topography can be affected by seismic amplification. Therefore, design codes account for topographic irregularities considering scaling factors depending on the topographic slope surrounding the studied site. This is particularly important in the framework of seismic hazards and for cultural heritage maintenance and prevention since many historical and ancient settlements in Italy were built on the top of hills, for defense reasons. The topic is complex and has been under debate for the last 5 decades. Seminal papers explained the effect in terms of constructive interference of seismic waves diffracted by the convex shape of topography (“topo-resonant model”, e.g., Géli et al., 1988). However, recent studies underlined that when considering a large number of sites (e.g., Burjánek et al., 2014a; Burjanek et al., 2014b; Pischiutta et al., 2018; Kaiser et al., 2022), this model is often not satisfied because a significant role is played by local velocity distribution and geological setting, as 1) large-scale open cracks (Moore et al., 2011; Burjanek et al., 2012); 2) microcracks in fractured rocks associated to fault activity (Martino et al., 2006; Marzorati et al., 2011; Pischiutta et al., 2012, 2015, 2017); 3) rock instabilities (e.g., Del Gaudio et al., 2019).

In this study, we investigated how to insert in hybrid simulations such amplification effects, observed in different geological and morphological conditions. The modeling was obtained by merging the low-frequency contribution from the kinematic rupture model proposed by Tinti et al. (2016), and the
high-frequency contribution achieved by stochastic simulation performed through the EXSIM code (e.g., Boore, 2009), both at bedrock and the site. In stochastic simulations, we exploit the model parameters validated in previous work (Pischiutta et al., 2021), including regional-specific source scaling, attenuation parameters, and the source complexity. They demonstrated that such a model can adequately explain spectral amplitudes, temporal characteristics of observed seismograms, and detect near-source effects related to the distribution of asperities on the fault plane. However, Pischiutta et al. (2021) also proved that, despite the general good consistency, in some cases simulations were not able to reproduce particular features of the observed acceleration spectrum. They ascribed such discrepancies to the occurrence of site amplification effects that are not accounted for by the use of generic amplification curves obtained through the quarter wavelength technique, due to improper consideration of the site contribution in ground-motion amplification. Also, Boore (2013) has revealed their constraints inferred by the method, which smooths, underestimating the primary resonant peaks provoked by the strong bedrock/soil seismic impedance contrasts. Moreover, several investigations have also suggested the significance of the soil/bedrock impedance contrast, the thickness of soil, and soil belongings in representing the site response in terms of amplitude and frequency content (e.g., Akinci et al., 2021).

Therefore, in this work, we have adopted different amplification curves to include the site contribution, such as the generic site curve; the horizontal-to-vertical spectral ratio (HVSR); the site functions from spectral inversion techniques (GIT). Then, we evaluated their performance by comparing observed and simulated ground motion, calculating the residual and the bias as a function of frequency.

## 2 METHODS

Broadband synthetic motion was generated following a hybrid approach exploited in Akinci et al. (2017), Ojeda et al. (2021), and Pischiutta et al. (2021). Here, the low-frequency (LF) portion of the synthetics (below 1 Hz) was obtained from the rupture model published by Tinti et al. (2016). Conversely, the high-frequency (HF) portion of synthetics (over 1 Hz) was attained by using a stochastic finite-fault simulation model, based on dynamic corner frequency, explained in the following section (Motazedian and Atkinson 2005; Boore 2009). These two frequency portions were merged in the frequency domain at each station following Mai and Beroza (2003). First, the LF component was selected to merge with the HF part using the consistency of the plateau level of acceleration in the Fourier space. Two frequency values are considered, f1 and f2: below f1 the signal is 100% LF, over f2 the signal is 100% HF, and between f1 and f2 the LF and HF spectra should be identical. Considering that LF simulations are reliable up to 1 Hz, we adopted 0.3 and 0.8 Hz for f1 and f2, respectively. However, at stations in deep sedimentary basins where site effects are expected below 1 Hz, we tailored the choice of f1 and f2 values, adopting values of 0.2 and 0.6 Hz, respectively (ex. CLF, GBP, SULA, etc). In this, way we ensure that at these

stations, hybrid synthetics include site amplification effects that are accounted for only by the HF signals.

The HF and LF signals were synchronized using a long- and short-time average (LTA/STA) automatic picking algorithm. To avoid a mismatch in the plateau levels between the HF and LF spectrum, we rotated the two horizontal low-frequency components by increments of 1°. The application of this procedure resulted in hybrid broadband signals related to the horizontal components of ground motion. More details can be found in Akinci et al. (2017).

## 2.1 High-Frequency Stochastic Ground Motion Simulations

In order to simulate the strong ground motion of the Amatrice earthquake we applied the stochastic finite-fault method, and later examined the residual of the ground motions between observed and simulated ground motion parameters both in the time and frequency domain. We follow the approach and parameters already exploited in Pischiutta et al. (2021). The finite-fault simulation employed the EXSIM code, produced by Motazedian and Atkinson (2005) and revised by Boore (2009), which requires as input model parameters a region-specific source model, path, and site contributions.

The total spectrum of the ground motion in this approach consists of earthquake sources, paths, sites, and instruments. These terms can be included in a comprehensive equation in the frequency domain, as follows:

$$A(M_0, r, f) = E(M_0, f) \cdot P(R, f) \cdot G(f) \cdot I(f) \quad (1)$$

where A(M₀, r, f), is the Fourier spectra acceleration, M₀, f and R are the seismic moment, the corner frequency, and the hypocentral distance from the observation point, respectively. The term E(M₀, f) is the earthquake source spectrum, and the term P(R, f) is the path that models the geometric spreading and the anelastic attenuation effects as a function of R and f. The term G(f) is the site effect and I(f) is the instrumental transfer function. The exploited parameters required by the method for the source, the path, and the site terms, are derived from those several current models published in Central Italy. Their adequateness was confirmed in Pischiutta et al. (2021), where the general good consistency found between synthetic and observed ground motion (both in the time and frequency domains) demonstrated that this model can adequately explain spectral amplitudes, temporal characteristics of observed seismograms, and to detect near-source effects related to the distribution of asperities on the fault plane. In order to show the performance of the adopted parameters in reproducing ground-motion estimates for the Amatrice earthquake, in Supplementary Figure S1 we plot the simulated hybrid broadband PGAs, PGVs up to 150 km as a function of RJB for the seismic stations in the “A” site class (blue reverse triangles), and for 961 virtual stations (turquoise dots) distributed in a 4-km grid space (see Pischiutta et al., 2021). Three Italian GMPE models are plotted as well, together with their ±σ standard deviations: Bindi et al., 2011 (cyan); Malagnini et al., 2011 (green); Sgobba et al., 2020 (violet). Simulated PGAs and PGVs lie within the standard deviation of
the three considered GMPEs, suggesting the adequateness of the adopted regional-specific source scaling and attenuation parameters in hybrid simulations in satisfactory reproducing ground-motion estimates.

## 2.1.1 Source

We adopt the same spectral parameters describing the earthquake source, employed and validated in Pischiutta et al. (2021). They are tabulated in Supplementary Table S1. The stress drop parameter $\sigma$, which rules the levels of the acceleration spectrum at high frequencies, was calculated at about 150 bars by Malagnini and Munafò, 2018, and Akinci et al. (2021). We also adopted the kinematic rupture model proposed by Tinti et al. (2016), divided into 0.5 km × 0.5 km sub faults along the strike and dip. According to earlier studies, source parameters, such as geometry (strike 156°, dip 50°), density (2.8 km/m³), and rupture propagation velocity (3.1 km/s), are employed among ordinarily referred values (Supplementary Table S1).

## 2.1.2 Propagation and Attenuation

Seismic wave propagation and seismic attenuation are essential topics, and they are required for the earthquake ground motion estimations in seismic hazard analysis. In our study, we decided to use the most recently described seismic attenuation parameters presented in Malagnini et al. (2011) model. It was obtained from several regressions of 170 weak-motion records belonging to foreshocks and aftershocks of the 2009 L'Aquila seismic sequence and realized through regression analyses of velocity time-series and Fourier spectra from 0.1 to 10 Hz, recorded at distances between 40 and 350 km. It provides the average features of three contributions in the wave propagation: geometrical attenuation, anelastic attenuation, and ground motion duration.

The chosen spectral parameters for seismic wave propagation are reported in Supplementary Table S2. The path spectrum, $P(R, f)$, depends on geometrical spreading, $Z(R, f)$, and quality factor ($Q_s$).

$$P(R, f) = Z(R) \cdot \exp\left(-\frac{\pi f R_{i,j}}{Q_s}\right) \quad (2)$$

Concerning the geometrical spreading coefficient $Z(R, f)$, we adopted a conventional piecewise function as represented by $r^{-1.1}$ at distances smaller than 10 km as a body-wave-like function; within 10 and 40 km, it is defined as $r^{-1}$; within 40 and 100 km, it is defined as $r^{-0.7}$; beyond 100 km distance, it is characterized by $r^{-0.5}$, that is compatible with the surface waves attenuation characteristics in a half-space.

A power-law pattern of $Q_s$ gives the anelastic attenuation:

$$Q_s = Q_0 f'' \quad (3)$$

where $Q_0$ is the value of $Q_s$ at a frequency of 1 Hz, and $\eta$ is the frequency parameter proposed by Aki and Chouet (1975). For the Central Italy region, the quality factor at frequencies > 0.6 Hz is given by

$$Q_s(f) = 140 f^{0.25} \quad (4)$$

Values adopted for frequency < 0.6 Hz are given in Supplementary Table S2. However, at such low frequency ranges the hybrid motion is dominated by the LF contribution due to the merging procedure previously explained.

## 2.1.3 Site Amplification

In the code EXSIM, site amplification is accounted for through the combination of the amplification $A(f)$ and attenuation $D(f)$ contributions as follows:

$$G(f) = A(f) \cdot D(f) \quad (5)$$

$D(f)$ is a diminution operator accounting for deamplification effects from the near-surface:

$$D(f) = \exp(-\pi \kappa_0 f) \quad (6)$$

where an exponential decay marks the kappa parameter ($\kappa_0$), representing the slope of the high-frequency declines of spectra in the stochastic finite-fault method (Anderson and Hough, 1984).

## 2.2 Including the Site Effects

Here, we aimed to involve the site amplifications determined from different approaches to testing the performance and using such generic, empirical, and “specific” site amplification curves at stations selected in conditions potentially prone to experiencing site effects (ex. basin, topography). We used the following site amplification curves, determined throughout commonly affirmed techniques:

1. Generic site curves employed in Pischiutta et al. (2021), and representative of NTC-18 classes A, B, C, and D. According to Eq. 5, they are composed of the product of the wave amplification term $A(f)$ and the diminution term $D(f)$, this latter accounting for high-frequency attenuation.

Many studies in the literature have provided generic amplification curves for the term $A(f)$, in the framework of the National Earthquake Hazard Reduction Program, NEHRP which is the seismic code adopted in the United States (e.g., Boore and Joyner, 1997; Boore, 2003; Boore 2016; Campbell and Boore, 2016). However, considering that the class thresholds for Italian and United States seismic design codes are different (see Table 4 in Pischiutta et al., 2021), in this study we choose to exploit such curves only for NTC-18 classes -B and -C (corresponding to NEHRP -C and -D, respectively). In particular, we adopted two $A(f)$ curves proposed in Boore and Joyner, 1997 without particular frequency peaks and with associated $Vs_{30}$ parameter lying in the middle of the ranges allowed in each class (520 and 255, respectively, for representative curves chosen for -B and -C classes). They are reported in Supplementary Table S3 and graphed in Supplementary Figure S2 (green and red continuous lines). Considering the differences between $Vs_{30}$ thresholds in the Italian NTC-18 and United States seismic codes, for NTC-18 class-A we adopted an $A(f)$ curve generated in Pischiutta et al. (2021) using the quarter wavelength approach (Boore, 2003; Boore, 2005) considering a typical velocity profile for Italian soft rocks (limestones, marls, and flysch), without strong impedance contrasts. It is reported in Supplementary Table S4 and shown
<figure><img src="image_55.png" /><figcaption>**FIGURE 1** | Map showing stations considered in this study up to a distance of 150 km to the epicenter of the Mw6.0 Amatrice earthquake (black star), and lying in particular site conditions, such as sedimentary basins (circle), topographic irregularities (square), fault zones (triangles). They belong to RAN ("IT") and RSN ("IV") networks (Table 1). We also add the fault surface projection of the causative fault proposed in Tinti et al. (2016). Symbol color is related to site classification according to NTC18 based on the V<sub>S30</sub> parameter (blue = class A, green = class B, red = classes C, and pink = D).</figcaption></figure>

in **Supplementary Figure S2** (blue continuous line). Similarly, for the NTC18 class-D (corresponding to NEHRP class-E), we exploited another A(f) curve generated by Pischiutta et al. (2021) considering a velocity profile with Vs<sub>30</sub> in the ranges prescribed by NTC18 (<180 m/s) and without significant impedance and velocity contrasts (see also **Supplementary Table S4** and the orange continuous line in **Supplementary Figure S2**). We remark that generic site amplification curves determined throughout the quarter wavelength procedure are represented by velocity gradients being insensitive to impedance variances connecting the layers (Joyner and Fumal, 1984; Boore and Joyne, 1997; Boore et al., 1994, 2011).

To calculate the diminution term D(f) representing high-frequency attenuation, we used the following values for the kappa parameter (κ<sub>0</sub>):

$$ \kappa_0 = \begin{cases} 0.02 \text{ s} & \text{for } \text{class } -A \\ 0.03 \text{ s} & \text{for } \text{class } -B \\ 0.04 \text{ s} & \text{for } \text{class } -C \\ 0.045 \text{ s} & \text{for } \text{class } -D \end{cases} \qquad (7) $$

The G(f) curves obtained following Eq. 5 are reported in **Supplementary Figure S2** through dotted lines.

2. Generalized inversion technique (GIT, Andrews 1986; Castro et al., 1990). This approach is a reference site
<table><caption>TABLE 1 | Station sample considered in this study.</caption><thead><tr><th>Net code</th><th>Station code</th><th>Station name</th><th>Latitude [°]</th><th>Longitude [°]</th><th>Position</th><th>EC8 code</th></tr></thead><tbody><tr><td>IT</td><td>ANT</td><td>ANTRODOCO</td><td>42,4182</td><td>13,0786</td><td>Fault zones</td><td>A</td></tr><tr><td>IT</td><td>AQG</td><td>L AQUILA V. ATERNO COLLE GRILLI</td><td>42,373474</td><td>13,337026</td><td>Sedimentary basin</td><td>B</td></tr><tr><td>IT</td><td>AQK</td><td>L AQUILA V. ATERNO AQUIL PARK ING.</td><td>42,344967</td><td>13,400949</td><td>Sedimentary basin</td><td>B</td></tr><tr><td>IT</td><td>AQV</td><td>L AQUILA V. ATERNO CENTRO VALLE</td><td>42,377222</td><td>13,343888</td><td>Others</td><td>B</td></tr><tr><td>IT</td><td>ASP</td><td>ASCOLI PICENO</td><td>42,848</td><td>13,6479</td><td>Sedimentary basin</td><td>B</td></tr><tr><td>IT</td><td>ASS</td><td>ASSISI</td><td>43,074982</td><td>12,604141</td><td>Others</td><td>A</td></tr><tr><td>IT</td><td>ATN</td><td>ATINA</td><td>41,620319</td><td>13,801154</td><td>Fault zones</td><td>A</td></tr><tr><td>IT</td><td>AVZ</td><td>AVEZZANO</td><td>42,0274</td><td>13,4259</td><td>Others</td><td>C</td></tr><tr><td>IT</td><td>BSS</td><td>BUSSI</td><td>42,191732</td><td>13,845266</td><td>Fault zones</td><td>A</td></tr><tr><td>IT</td><td>BTT2</td><td>BORGO OTTOMILA - 2 (CELANO)</td><td>41,998333</td><td>13,543056</td><td>Sedimentary basin</td><td>D</td></tr><tr><td>IT</td><td>BVG</td><td>BEVAGNA</td><td>42,932367</td><td>12,611065</td><td>Sedimentary basin</td><td>C</td></tr><tr><td>IT</td><td>CCT</td><td>CITTA DI CASTELLO (TRESTINA)</td><td>43,3683</td><td>12,2346</td><td>Sedimentary basin</td><td>C</td></tr><tr><td>IT</td><td>CLF</td><td>COLFIORITO</td><td>43,036714</td><td>12,920428</td><td>Sedimentary basin</td><td>D</td></tr><tr><td>IT</td><td>CSA</td><td>CASTELNUOVO ASSISI</td><td>43,008015</td><td>12,590602</td><td>Sedimentary basin</td><td>C</td></tr><tr><td>IT</td><td>CTL</td><td>CATTOLICA</td><td>43,955116</td><td>12,735809</td><td>Sedimentary basin</td><td>C</td></tr><tr><td>IT</td><td>CTS</td><td>CITTA DI CASTELLO REGNANO</td><td>43,491987</td><td>12,223396</td><td>Sedimentary basin</td><td>C</td></tr><tr><td>IT</td><td>FBR</td><td>FABRIANO</td><td>43,343601</td><td>12,9119</td><td>Sedimentary basin</td><td>C</td></tr><tr><td>IV</td><td>FEMA</td><td>Monte Fema</td><td>42,9621</td><td>13,04976</td><td>Topography</td><td>A</td></tr><tr><td>IT</td><td>FOC</td><td>FOLIGNO</td><td>43,0263</td><td>12,896506</td><td>Sedimentary basin</td><td>C</td></tr><tr><td>IT</td><td>FOS</td><td>FOLIGNO SEGGIO</td><td>43,01459</td><td>12,83513</td><td>Topography</td><td>B</td></tr><tr><td>IV</td><td>GAG1</td><td>Gagliole</td><td>43,238063</td><td>13,067434</td><td>Topography</td><td>A</td></tr><tr><td>IT</td><td>GBB</td><td>GUBBIO</td><td>43,356972</td><td>12,597252</td><td>Others</td><td>B</td></tr><tr><td>IT</td><td>GBC</td><td>Gubbio</td><td>43,355301</td><td>12,5726</td><td>Sedimentary basin</td><td>C</td></tr><tr><td>IT</td><td>GBP</td><td>GUBBIO PIANA</td><td>43,31381</td><td>12,58949</td><td>Sedimentary basin</td><td>C</td></tr><tr><td>IT</td><td>GRN</td><td>GUARCINO</td><td>41,8134</td><td>13,3169</td><td>Topography</td><td>A</td></tr><tr><td>IT</td><td>LDP</td><td>LAMA DEI PELIGNI</td><td>42,0392</td><td>14,1826</td><td>Topography</td><td>C</td></tr><tr><td>IT</td><td>LSS</td><td>LEONESSA NUOVA</td><td>42,558243</td><td>12,968894</td><td>Others</td><td>A</td></tr><tr><td>IV</td><td>MDAR</td><td>Monte Daria</td><td>43,1927</td><td>13,1427</td><td>Topography</td><td>B</td></tr><tr><td>IV</td><td>MGAB</td><td>Montegabbione</td><td>42,91263</td><td>12,11214</td><td>Topography</td><td>A</td></tr><tr><td>IT</td><td>MMP1</td><td>MOMPEO 1</td><td>42,249229</td><td>12,748319</td><td>Topography</td><td>A</td></tr><tr><td>IV</td><td>MMUR</td><td>Monte Murano</td><td>43,44183</td><td>12,9973</td><td>Topography</td><td>B</td></tr><tr><td>IV</td><td>MNTP</td><td>Montappone</td><td>43,137378</td><td>13,469252</td><td>Topography</td><td>B</td></tr><tr><td>IT</td><td>MTR</td><td>MONTEREALE</td><td>42,524</td><td>13,2448</td><td>Topography</td><td>B</td></tr><tr><td>IV</td><td>MURB</td><td>MONTE URBINO</td><td>43,263</td><td>12,5246</td><td>Topography</td><td>A</td></tr><tr><td>IT</td><td>MVB</td><td>MARSCIANO MONTE VIBIANO</td><td>42,9619</td><td>12,257</td><td>Others</td><td>A</td></tr><tr><td>IT</td><td>NRN</td><td>NARNI</td><td>42,51556</td><td>12,51944</td><td>Topography</td><td>A</td></tr><tr><td>IT</td><td>PGG</td><td>POGGIO PICENZE</td><td>42,322872</td><td>13,539446</td><td>Sedimentary basin</td><td>B</td></tr><tr><td>IV</td><td>PIEI</td><td>PIEIA</td><td>43,53567</td><td>12,535</td><td>Topography</td><td>A</td></tr><tr><td>IT</td><td>PNN</td><td>PENNABILLI</td><td>43,818159</td><td>12,262846</td><td>Topography</td><td>C</td></tr><tr><td>IV</td><td>PP3</td><td>PP3</td><td>43,37783</td><td>13,6095</td><td>Sedimentary basin</td><td>C</td></tr><tr><td>IT</td><td>PSC</td><td>PESCASSEROLI</td><td>41,812042</td><td>13,789196</td><td>Fault zones</td><td>A</td></tr><tr><td>IV</td><td>RM33</td><td>PELLESCRITTA</td><td>42,50898</td><td>13,21452</td><td>Topography</td><td>A</td></tr><tr><td>IT</td><td>SBC</td><td>SUBIACO</td><td>41,9132</td><td>13,1055</td><td>Others</td><td>A</td></tr><tr><td>IT</td><td>SCF</td><td>SCAFA</td><td>42,265117</td><td>13,998489</td><td>Others</td><td>B</td></tr><tr><td>IT</td><td>SNG</td><td>SENIGALLIA</td><td>43,68558</td><td>13,226162</td><td>Others</td><td>C</td></tr><tr><td>IT</td><td>SNS1</td><td>SANSEPOLCRO 2</td><td>43,573502</td><td>12,1312</td><td>Sedimentary basin</td><td>C</td></tr><tr><td>IT</td><td>SPM</td><td>SPOLETO MONTELUCO</td><td>42,72324</td><td>12,751268</td><td>Fault zones</td><td>A</td></tr><tr><td>IT</td><td>SRL</td><td>SIROLO</td><td>43,517905</td><td>13,619388</td><td>Others</td><td>C</td></tr><tr><td>IT</td><td>SUL</td><td>SULMONA</td><td>42,089</td><td>13,934</td><td>Others</td><td>A</td></tr><tr><td>IT</td><td>SULA</td><td>SULMONA AUTOPARCO</td><td>42,0734</td><td>13,9166</td><td>Sedimentary basin</td><td>C</td></tr><tr><td>IT</td><td>SULC</td><td>SULMONA CONSORZIO</td><td>42,068</td><td>13,909</td><td>Sedimentary basin</td><td>C</td></tr><tr><td>IT</td><td>TLN</td><td>TOLENTINO</td><td>43,215904</td><td>13,25838</td><td>Topography</td><td>A</td></tr><tr><td>IT</td><td>TRE</td><td>TREVI</td><td>42,876499</td><td>12,7358</td><td>Sedimentary basin</td><td>C</td></tr><tr><td>IV</td><td>TRE1</td><td>Treia</td><td>43,311198</td><td>13,312848</td><td>Topography</td><td>B</td></tr><tr><td>IT</td><td>TRL</td><td>TERMINILLO</td><td>42,461314</td><td>12,932308</td><td>Others</td><td>B</td></tr><tr><td>IT</td><td>TRN1</td><td>TERNI 2</td><td>42,558201</td><td>12,6461</td><td>Sedimentary basin</td><td>D</td></tr><tr><td>IT</td><td>TVL</td><td>Tivoli</td><td>41,893015</td><td>12,773221</td><td>Sedimentary basin</td><td>B</td></tr></tbody></table>

method such as the standard spectral ratio method (SSR) but, in contrast with the SSR, the GIT curves are the result of a nonparametric inversion scheme applied to a dataset composed of multiple events and stations. In this

study, we used the GIT amplification functions obtained from a dataset of 283 stations and 455 events that occurred in Central Italy (Morasca et al., 2022). The GIT analysis was performed on a frequency range of
<figure><img src="image_56.png" /></figure>

**FIGURE 2** | Exemplification simulation results at station BSS (Bussi), belonging to the IT network. The left-top panel shows amplification curves applied as the term G(f) in the EXSIM code used for stochastic simulations (Eq. 5). The red curve is the generic one that we used for all class B sites defined in NTC-18 based on the V₅₃₀ value (see also **Supplementary Figure S3**). The other tested site-empirical specific curves are: the cyan curve is the HVSR calculated using earthquake weak motion (S waves); the green and turquoise curves are the GIT-derived amplification function by using 455 earthquakes and only the Amatrice earthquake, respectively (Morasca et al., 2022); the blue curve is the HVNSR calculated using ambient noise or earthquake coda waves. Synthetic hybrid broadband horizontal component velocity time histories (left panel) and Fourier velocity amplitude spectra (top-right panel) are compared with the recorded ground motions (black lines). For the latter, we calculated the geometric mean, and mean of the two horizontal components of ground motion, respectively, for time histories and Fourier spectra.

0.5–25 Hz and considered a hypocentral distance range of 10–120 km. The solution of the linear system composed of three terms (source, attenuation, and site contributions) required two prior constraints to remove unresolved degrees of freedom. A first assumption is that for all frequencies the attenuation term is set to unity at the reference distance of 10 km (the smallest in the dataset). The second one is a reference site condition. Considering six reference sites (LSS, MNF, NRN, SNO, SDM, and SLO) located on the rock and carefully selected on the base of Lanzano et al., 2021 analysis, their average amplification is fixed to 1, removing the linear dependence between source and site terms. We also exploit the GIT amplification function obtained only from the Amatrice earthquake (when available). In **Supplementary Figure S3** we provide such GIT curves at the eleven stations that are thoroughly investigated in this paper. At some of them (CLF, FOC, TRE, FEMA, BSS, SULA, GBB) the two GIT

curves show differences both in terms of amplified frequency band and amplitude level. This prompted us to test both of them in stochastic simulations, in a way to better reproduce observed motions.

3. Horizontal-to-vertical (HV) spectral ratio computed on the Fourier Amplitude Spectra (FAS) of strong-motion recordings (HVSR) and ambient noise (HVNSR). Although these methods are usually adopted to estimate the resonance frequencies of the site, we test if they can be also adopted to evaluate the site amplification in specific geomorphological conditions (Molnar et al., 2018; Kawase teal., 2019; Zhu et al., 2020). The HV method may estimate the site amplification function if the vertical amplification is negligible (Lachet and Bard, 1994; Field and Jacob, 1995). The HV curves are available in the Italian Accelerometric Archive v 3.2 (ITACA, Russo, et al., 2022). If the HVNSR curve is missing, we select the HV computed on the coda waves since many studies demonstrated that the curves obtained from noise and coda waves are comparable. Instead, if the HVSR curve is not available in ITACA we use the HV curves
<figure><img src="image_57.png" /></figure>

**FIGURE 3** | Simulation results and records of Amatrice earthquake at Gubbio [panel (A)] and Sulmona sedimentary basins [panel (B)]. For each site, we add: a rough geological map of the area (provided in the ITACA database); the applied generic and empirical site curves at each station (consistently with **Figure 2**); Synthetic velocity Fourier spectra were obtained using the different site curves and recorded velocity Fourier spectra (geometric mean of the two horizontal components).

from an independent study (Priolo et al., 2019) estimated using records of small events (magnitude < 4.5).

Adopted site amplification curves were not concurrently available at all stations. The kappa coefficient $\kappa_0$ was not applied for the simulations for the GIT, HVRS, and HVNSR experimental curves, since implicitly considers the total attenuation effects.

## 2.3 Dataset and Recording Stations

Among stations recording the Amatrice 24 August 2016 earthquake at 150 km from the epicenter, we selected a

subset of 57, installed in different site conditions potentially prone to experience site amplification effects (Figure 1) (Table 1), as

1. Twenty-two (22) stations in sedimentary basins (circles in Figure 1) mostly related to site classes C and D, with only a minor percentage (five stations) in class B. An inspection of the HVSRs and HVNSRs published on the ITACA database resulted in relevant amplification (exceeding a factor of three).

2. Eighteen (18) stations are located on topographic irregularities (squares in Figure 1). Most parts of them lie in site class A (10) and B (6).
A Velocity time histories

<figure><img src="image_58.png" /></figure>

B

<figure><img src="image_59.png" /></figure>

**FIGURE 4** | Velocity time series at Gubbio [panel (A)] and Sulmona [panel (B)] sites. In black, we plot recorded seismograms (arithmetic mean of the two horizontal components). Colored lines depict synthetic hybrid seismograms obtained from the HF and LF contributions (see **Section 2.3**), obtained through the EXSIM code seismic source inversion by Tinti et al. (2016). They were obtained by using empirical site curves in the EXSIM code: GIT (green), HVSR (cyan), and HVNSR (blue).

3. Six (6) stations on class A rock sites and close to fault zones (as identified by a visual inspection of geological maps and information archived in the ITACA database).

4. Eleven (11) and other stations belong to A, B, and C sites (no faults/topography/sedimentary basins) as a further constraint.

Such a station sample is a portion of the one used in Pischiutta et al. (2021) work, where 133 stations at 150 km from the epicenter were selected to validate stochastic simulations through comparison with recorded data in this work.

3 RESULTS

A comparison between observed and simulated ground motion was
led to get insight into simulation reliability and the performance of the
different adopted site amplification curves. **Figure 2**, reports
exemplificative simulation results at station BSS (Bussi) of the IT
network, related to site class A and located close to a tectonized zone.
The left-top panel shows amplification curves used for stochastic

simulations (Eq. 5). The red curve is the generic one that we
used for all class B sites defined in NTC-18 based on the V<sub>s30</sub>
value (see also **Supplementary Figure S2**). According to Eqs 5,
and 6, it is obtained by the product of A(f) and D(f). For
A(f), we adopted as representative of NTC-18 class-B an
amplification curve proposed in Joyner and Boore (1997) and
related to V<sub>s30</sub> of 520 m/s. D(f) was calculated using Eq. 6
considering &kappa;<sub>0</sub> a value of 0.03 s. The other tested site-empirical
specific curves were applied as G(f), considering that
attenuation contribution should be implicitly included:

1) The cyan curve is the HVSR calculated using S-waves of seismic recordings

2) The green and turquoise curves are the GIT-derived amplification for the station by using, respectively, 455 events (Morasca et al., 2022), the empirical function resulting from average overall events for robustness reasons, and only the Amatrice 24 August 2015 earthquake

3) The blue curve is the HVNSR calculated using ambient noise or earthquake coda waves
<figure><img src="image_60.png" /><figcaption>**FIGURE 5** | Bias calculated following Eq. 10, to get an insight into the performance of the generic and empirical site amplification curves exploited in the EXSIM code all over the station sample in sedimentary basins. The top panels show bias versus frequency for both acceleration (A) and velocity (B) spectra. Black curves represent the bias obtained by applying the generic site curve, while color curves are related to the use of GIT (green), HVSR (cyan), and HVNSR (violet) site curves. The standard deviation is plotted as well. The bottom panels show the bias calculated for peak ground acceleration [PGA, panel (B)] and peak ground velocity [PGV, panel (C)].</figcaption></figure>

To validate the effectiveness of our simulations to reproduce observations and test the different site curves, we compared the synthetic velocity time histories (right panel) to recorded horizontal ground motion (arithmetic mean). Velocity Fourier amplitude spectra (left-bottom panel) were also compared with recorded ones (geometric and mean of the two horizontal components of ground motion).

The HVNSR curve and the GIT curve obtained from the Amatrice earthquake consistently show a peculiar characteristic concerning other amplification curves presenting two prominent peaks between 1 and 2 Hz at station BSS. These peaks are easily observed in the simulated time histories, as shown in the left-bottom panel of Figure 2.

As a second step, to quantitatively assess simulations' overall performance we calculated residuals $RF_j(f)$ between observed $F_j$ and simulated $\underline{F_j}$ spectra considering the different generic and empirical curves, related to station j as

$$RF_j(f) = \left( \frac{F_j}{\underline{F_j}} \right) \qquad (8)$$

Similarly, we calculated the residual $R_j$ between observed $Y_j$ and simulated $\underline{Y_j}$ ground-motion parameters (PGA and PGV) at each station j as

$$R_j = \frac{Y_j}{\underline{Y_j}} \qquad (9)$$

We finally computed the bias averaging over the total station sample N:

$$\begin{align}
BIAS RF(f) &= \frac{1}{N} \sum_{j=1}^{N} \ln(R_j(f)) \nonumber \\
BIAS R &= \frac{1}{N} \sum_{j=1}^{N} \ln(R_j) \tag{10}
\end{align}$$

A perfect match between the empirical model and the broadband simulation would have null BIAS values, whereas positive/negative residual shows an underprediction/overprediction of the simulations concerning observed ground motion.
A

<figure><img src="image_61.png" /></figure>

Observed ground motions

Synthetic ground motions obtained by using the following site curves:

B

<figure><img src="image_62.png" /></figure>

Generic HVSR

GIT HVNSR

**FIGURE 6** | Simulation results and records of the Amatrice earthquake at stations GAG1 [panel (A)] and MURB [panel (B)], located on topographic irregularities. For each one, we show a digital elevation model to evince topography shape; the applied generic and empirical site curves at each station (consistently with **Figure 2**); Synthetic velocity Fourier spectra were obtained using the different site curves and recorded velocity Fourier spectra (geometric mean of the two horizontal components).

## 3.1 Stations in Sedimentary Basins

In order to thoroughly investigate how many stratigraphic effects are accounted for/unaccounted for in the earthquake-induced ground motions, we selected twenty-two stations in sedimentary basins.

**Figure 3** details results at two sedimentary basins where Gubbio (panel A) and Sulmona (panel B) towns are settled. Gubbio is a historical town in the Umbria region, central Italy, with more than 30,000 inhabitants and a rich artistic and cultural heritage dating to the Middle-Age. The historical part was built at the lower slope of Ingino hill, but the modern portion expanded towards the alluvial plain. In this area, three stations of the Italian seismic network IT are installed: GBP in the middle of the alluvial basin (site class C); GBC on the basin border, inside Gubbio settlement (site class C); GBB outside the alluvial basin (site class B). Bindi et al. (2009) observed that time series of local earthquakes recorded in the Gubbio plain are characterized by locally generated surface waves, which increase in duration and amplitude with respect to the nearby reference station on a rock (GBB), the spectral energy is distributed over the range 0.4–2 Hz. In particular, they found that the peak ground velocity is amplified by a factor of 5, and the duration is increased by a factor of about 2 where the sedimentary cover is thickest (ca 600 m).

In the top panel A of **Figure 3**, we report a basic geological map of the area. We also add the site curves for each station, consistently with **Figure 2**. While the generic curve (red) which does not account for the specific basin structure, does not show any significant peaks, the GIT curves (green and turquoise) show a relevant amplitude peak at about 5 Hz and between 2 and 10 Hz at stations GBB and GBP, respectively. The use of GIT curves led to better reproduction of the general spectral trend and velocity

recorded signals (**Figure 4A**). At station GBP, in the middle of the basin, the HVSR and HVNSR curves show a peak at about 0.35 Hz, up to a factor of 4 and 7, respectively. In **Figure 4**, we include synthetic velocity time series for Gubbio (panel A) and Sulmona (panel B) cases. Especially by using the HVNSR and GIT curves, we obtained synthetic time histories with amplitudes and spectral content more consistent with recorded signals. At station GBC we applied only the empirical HVSR curve together with the generic one. Bindi et al. (2009) highlighted that in the Gubbio basin the HVSR are strongly affected by amplification on the vertical component, and this is particularly evident close to the basin border. In this study, we found that at station GBC the use of the HVSR curve did not lead to obtaining synthetic signals consistent with observed data.

The sedimentary basin hosting the Sulmona town in the Abruzzi region hosts more than 20,000 inhabitants (**Figure 3B**). Sulmona was founded in Roman times and holds a rich artistic and cultural heritage dating from the Middle-Age. Sulmona rises in the center of the Peligna Valley, which in prehistoric times was occupied by a vast lake. It is located between the Vella torrent and the Gizio rivers, to the west of the Maiella and Morrone mountains, which overlook the city.

In this study, we consider four stations of the Italian seismic network IT: SULA, in the middle of the alluvial basin (site class C); SULC and SULP, on the basin borders (site class C and B, respectively); SUL, outside the alluvial basin (site class A). In **Figure 3**, we also show the geological map of the area, as well as the adopted site curves at each station. While the generic (red) and HVSR (cyan) site curves do not show any significant peaks, the GIT curves at stations SULA, SULC, and SULP show a peak at about 1 Hz, with amplitudes varying from 5 close to the basin borders (SULP and SULP) to over 8 in the middle of the basin
<figure><img src="image_63.png" /><figcaption>**FIGURE 7** | Bias calculated following Eq. 10, to get an insight into the performance of the generic and empirical site amplification curves exploited in the EXSIM code all over the station sample on topography. The top panels show bias versus frequency for both acceleration (A) and velocity (C) spectra. Black curves represent the bias obtained by applying the generic site curve, while color curves are related to the use of GIT (green), HVSR (cyan), and HVNSR (violet) site curves. The standard deviation is plotted as well. The bottom panels show the bias calculated for peak ground acceleration [PGA, panel (B)] and peak ground velocity [PGV, panel (D)].</figcaption></figure>

(SULA). The HVNSR curves show another peak at about 0.45 Hz at SULA and SULC, with amplitudes decreasing from 8 to 6. It is also evident on the HVSR curve even if with lower amplitude levels (ranging from 3 to 4). By applying the GIT site curves we obtained synthetic velocity spectra reproducing the observed spectral trend on data at stations SULA, SULP, and SULC. At this latter station, a better performance was achieved by adopting the HVSR and HVNSR. Finally, at station SUL there are no differences in synthetic spectra obtained by adopting the different site curves (velocity time histories are given in the **Supplementary Figure S4**.

We finally highlight that at stations in the middle of the two investigated sedimentary basins (GBP and SULA), recorded

seismograms show the presence of low-frequency phases, mainly in the coda, due to 3D amplification caused by the impedance contrast between soft sediments and underlying rock formations and by wave reverberation across the sedimentary body in the basin (e.g., Cornou and Bard, 2003). Such empirical 3D effects cannot be thoroughly simulated by employing the stochastic finite-fault technique implemented in this study. However, these effects are not observed at stations installed near the sedimentary basin borders (GBC, SULC, SULP). Finally, stations outside sedimentary basins (GBB and SUL) do not significantly improve by using different amplification curves in fitting the observed time histories. In
A

<figure><img src="image_64.png" /></figure>

Observed ground motions

Synthetic ground motions obtained by using the following site curves:

B

<figure><img src="image_65.png" /></figure>

Generic HVSR

GIT(Average) HVNSR

GIT(Amatrice Eq.) HVNSR

**FIGURE 8** | Simulation results and records of Amatrice earthquake at stations LSS [panel (A)] and SPM [panel (B)], installed on rock and related to site class A. For each one, we show: the geological map and profile furnished in the ITACA database; the applied generic and empirical site curves at each station (consistently with **Figure 2**); Synthetic velocity Fourier spectra were obtained using the different site curves, and recorded velocity Fourier spectra (geometric mean of the two horizontal components).

fact, since they are installed outside the sedimentary basin, they are not affected by any significant site amplification effects.

In **Supplementary Figure S5**, we also provide results from two other study cases in sedimentary basins: Colfiorito and Foligno. They are both located in the Umbria region (central Italy) and experienced an MCS damage intensity of VII during the 1997 seismic sequence (e.g., Camassi et al., 2008), where over six shocks with magnitude larger than 5.0 were produced. In the Colfiorito basin, two stations of the Italian seismic network “IT” are installed: FOC (close to the basin border) and CLF (in the middle of the basin). The former station is inside the small village of Colfiorito, hosting about 5,000 inhabitants and some interesting monuments (i.e., a Church dating back to the fifth century). In **Supplementary Figure S5A**, we report several geological information included in the ITACA database: a geological map and section passing through station FOC and a velocity profile above station CLF. The latter highlights a strong velocity contrast at about 50 m depth, probably responsible for the about 1.0 Hz peak observed on both the GIT and HVNSR curves. However, at this station, empirical site-specific curves in stochastic simulations did not generally result in better reproducing observed data. Conversely, at station FOC the use of the GIT curves (in particular the ones derived by using only the Amatrice earthquake) led to obtaining synthetic spectra similar to the observed trend on recorded spectra at frequencies over 5 Hz.

Foligno municipality has more than 50,000 inhabitants and is located close to Colfiorito (20 km far), in the center of the Umbrian Valley that is crossed by the Topino river. Foligno possesses an important cultural heritage, with many civil and religious edifices dated since the Middle-Ages. The two stations installed in this area (BVG and TRE) show a relevant amplitude peak between 1 and 3 Hz, most prominent on the GIT curves (**Supplementary Figure S5B**). The use of the latter (in particular the ones derived by using only the Amatrice earthquake at station TRE) led to achieving the best consistency between synthetic and recorded spectra.

In order to get an insight into the performance of the generic and empirical site amplification curves exploited in the EXSIM code all over the station sample in sedimentary basins, we calculated the bias of the model as the logarithm (base n) of the ratio of the observed to the simulated following Eq. 10: it is an indication of the difference in the frequency domain between simulated and observed ground motion. In the top panels in **Figure 5**, we plot the bias versus frequency for both acceleration (panel A) and velocity (panel C) spectra. Colored curves represent the bias obtained between observed and simulated spectra by applying the generic site (red), GIT (green), HVSR (cyan), and HVNSR (violet) site curves. The standard deviation is plotted as well. An overall bias reduction is achieved using the GIT site curve, whose performance is better than the generic one at low (0.5–1.5 Hz) and high frequencies (>10 Hz). This is quite evident
in the bias related to velocity spectra (Figure 5B) and to the peak ground velocity (PGV, Figure 5D), where it is lower than 0.5. Moreover, at high frequencies (>10 Hz) there is a tendency in all empirical curves for the bias reduction, that can be ascribed to the implicit inclusion of the attenuation effect in empirical curves. Conversely, in generic curves attenuation was not derived but only hypothesized through the application of the $κ_0$ parameter, arbitrarily assigned on the basis of the site class (Section 2).

## 3.2 Stations on Topography

To investigate the topographic effects on the earthquake-induced ground motions, we have considered 18 stations located on topographic irregularities (middle slope or top), mainly belonging to class A ($V_{s30} > 800$ m/s) and class B ($360 < V_{s30} > 800$ m/s), as prescribed by the Italian NTC18 seismic design code.

**Figure 6** shows two examples, both stations belonging to seismic network IV. The former (panel A) is station GAG1 (B site class), installed in Gagliole (Marche region), a village built in the Middle-Age on the top of a 1,000 m high hill. The historic center retains the original urban structure with the ancient medieval castle, the walls dating back to the 14th century, some ancient churches, and a monastic complex dated back to the 12th century. The GIT site amplification curves show overall higher amplitudes than both the HVSR and generic ones, with a peak at about 3.5 Hz. Its use led to achieving synthetic spectra more consistent with observed ones.

The second example (Figure 6B) is station MURB, installed on the top of Monte Urbino, an 800 m-high uninhabited hill in the Umbria region. While the generic site curve (related to B class) attains amplitude levels lower than 2 (see also Supplementary Figure S2, green-dotted curve), the three- empirical site curves (GIT, HVSR, and HVNSR) show a similar trend, with an amplitude up to eight peaks between 3 and 4 Hz. We, therefore, got synthetic spectra more consistent with observed ones between 2 and 5 Hz. Finally, through the use of the HVNSR curve, synthetic signals better reproduce the observed spectral trend even at high frequencies (>10 Hz).

In the **Supplementary Material**, we show four other study cases. While station FEMA is located on an uninhabited prominent ridge (over 1,500 m high), stations MMP1, PIEI, and TLN are close to small towns and have less prominent topography. All stations site curves show peaks at site-specific frequencies (about 10 Hz at FEMA; 1.5–3 Hz at MMP1; about 3 and 10 Hz at PIEI; about 2 Hz at TLN). This led to generally obtaining synthetic spectra more similar to observations in these frequency bands. Moreover, as for stations in sedimentary basins, at high frequencies (>10 Hz) the use of empirical site curves led to achieving a spectral trend more consistent with observed data. This is particularly evident at station FEMA, where the use of the GIT curve derived from the Amatrice earthquake (turquoise) led to a better fit of the velocity Fourier spectra at 2 Hz, where a prominent peak is observed. Finally, at station MMP1, empirical curves led to overestimating observation, so the generic site curves became more appropriate.

In **Figure 7**, similar to **Figure 5**, we plot the bias following Eq. 10 to get an indication of the difference in the frequency domain between simulated and observed ground motion. At high frequencies (>10 Hz) empirical site curves (HVSR and GIT) are more performant in producing a better fit with data. In fact, while they led to bias values lower than 0.5, the use of the generic curve is associated with bias values increasing from 0.5 to 1 from 10 to 25 Hz. Conversely, at low frequencies, no improvements are observed in the use of empirical site curves. Finally, observed PGA and PGV values are slightly but systematically underestimated and overestimated by the simulations.

## 3.3 Stations on the Rock Site

To study the site amplification effect of the rocks we utilized five stations installed on rock and related to site class A. However the HVSR and HVNSR amplification curves, published on the ITACA database, unexpectedly highlight the presence of peaks over a factor of three at intermediate frequencies (between 1 and 10 Hz), where we consider no amplification.

In **Figure 8A**, we show station LSS as an example of a reference rock site (i.e., no site amplification) installed on Mesozoic limestone lithotypes. All site curves (both generic and empirical) show amplitudes lower than a factor of 1.5, the difference between them being small without sharp peak/s at certain frequencies. Therefore, no differences are observed on synthetic velocity Fourier spectra when using the different site curves. Conversely, SPM is an exemplificative station for site amplification on rock sites (Figure 8B). It is installed on Mesozoic limestone, in an intensely tectonized zone, and at the middle slope of a hill. Empirical site curves show an amplitude exceeding three peaks at about 5 Hz. However, visually speaking, we may say that the GIT curves lead to better reproduction of observed Fourier spectra, particularly at frequencies higher than 5 Hz.

In **Supplementary Figure S7**, we furnish two other examples. The former (panel A) is represented by station ANT, installed in the town of Antrodoco (Latium region) on dolomite lithotypes and in an intensely tectonized area. Using empirical site curves HVSR and HVNSR (showing the amplitude of three peaks at about 2 Hz) led to overestimating observed velocity Fourier spectra, while by using the GIT curve simulated Fourier spectra they are consistent with observed data, especially at high frequencies. Moreover, the use of the GIT curve derived only using the Amatrice earthquake led to better reproduction of observed ground-motion levels. The same findings are shown at station PSC (Figure 8B) installed in Pescasseroli village (Abruzzi region) on calcarenites.

## 4 CONCLUSION

We produced synthetic broadband seismograms using a hybrid simulation technique for the $M_{w6.0}$ Amatrice earthquake, Central Italy, on 24 August 2016, following the previous work of Pischiutta et al. (2021). In the present study, we focused on the site amplification parameters and tried to understand their
impact on ground motions with the aim to improve the hazard assessment for seismic risk reductions, particularly in urban areas in the central Apennines. To do so, the 57 stations selected, mainly located in urban areas, are potentially prone to experience site amplification effects because of lying in particular site conditions (sedimentary basins, topographic irregularities, and fault zones).

In this work, we tested the use of different empirical amplification curves, such as horizontal-to-vertical spectral ratios (calculated using both earthquakes, HVSR, and ambient noise recordings, HVNSR), and site curves derived from the generalized inversion technique (GIT) to improve our simulations. The latter were derived by using 455 earthquakes that occurred in Central Italy (Morasca et al., 2022). We also tested the use of the GIT curve derived only by the Amatrice earthquake. Their performance was linked to the generic curves by comparing the goodness of fit with recorded data. In general, we observed the following:

* In sedimentary basins, the presence of superficial soft sediments and strong shear-wave velocity and impedance contrasts leads to strong amplification due to refraction of seismic waves by an underlying stiff bedrock and subsequent phase constructive interference causing a resonance effect. We found that at the 22 stations in sedimentary basins selected in this study, empirical site curves led to a better fit with the data. They include amplitude peaks at site-characteristic frequencies, which depend on the superficial mean shear-wave velocities and the depth of the velocity contrast. In particular, the GIT site curve was most performant among the empirical adopted curves, achieving the best fit between observed and simulated velocity Fourier spectra. The GIT curve derived only by the Amatrice earthquake is more performant than the one derived as an average among considered earthquake data set in Morasca et al., 2022, suggesting that site correction may be scenario-dependent.

* An overall bias reduction is achieved using the GIT site curve, whose performance is better than the generic one. This behavior is quite evident on velocity spectra both at (i) low frequencies (0.5–1.5 Hz), where amplification effects are generally realized in deep sedimentary basins (hundreds of meters); (ii) high frequencies (>10 Hz), since all empirical curves implicitly include the attenuation effect in empirical curves that is hypothesized on generic curves through the application of $k_0$ parameter, arbitrarily assigned on the basis of the site class. However, due to the limitation intrinsic to the 1D stochastic finite-fault approach implemented in this study, our simulations could not reproduce wave reverberation across the sedimentary body, visible at stations installed in the middle of the basins on time histories as low-frequency coda phases.

* Amplification effects occur on the top of reliefs, due to the constructive interference of seismic waves diffracted by the convex topography. Almost certainly, a significant role is played by the local velocity distribution and subsoil structure. This study selects 18 stations installed on

topography at which empirical site curves show peaks at different frequencies which depend on site specificities. As for stations in sedimentary basins, even at these sites, the use of empirical curves led to achieving a spectral trend more consistent with observed data at high frequencies (>10 Hz). Conversely, at low frequencies, no general improvements are observed in the use of empirical site curves. In fact, each site shows its own peculiar behavior, and a preferential empirical/generic curve leading to achieving synthetic motion more consistent with observed data.

* Rock sites are considered to be unaffected by site amplification. Nevertheless, many recent studies have highlighted that seismic waves can be amplified due to the local properties of the rock (i.e., the presence of pervasive fractures and/or large open cracks). We investigate five stations installed on rock and related to site class A, where empirical site curves highlight the occurrence of amplification effects at intermediate frequencies (between 1 and 10 Hz). Again, we found that the GIT empirical site curves led to obtaining synthetic spectra more consistent with observed ones. This is particularly evident when using the GIT site curve derived only from the Amatrice earthquake suggesting that site correction may be scenario-dependent.

Similar efforts have been shown in different parts of the world. For example, recently Zhu et al. (2022a) have tested and compared the effectiveness of different estimation techniques using a unique benchmark dataset at 1725 K-NET and KiK-net sites, in Japan. Evaluated prediction approaches included: 1) the empirical correction to the horizontal-to-vertical spectral ratio of earthquakes (c-HVSR, see also Kawase et al., 2019 and Zhu et al., 2020); 2) one-dimensional ground response analysis (GRA); 3) the square-root-impedance (SRI) method (also called the quarter wavelength approach). They found that, at the majority of analyzed sites, the empirical correction to HVSR was highly effective in achieving a “good match” in both spectral shape and amplitude. Since this technique has great potential in seismic hazard assessments, even considering that it does not require a velocity model, its use could be evaluated in similar further tests involving hybrid simulations.

We finally stress that, due to the 1D stochastic finite-fault approach implemented in this study, our simulations could not reproduce several scattering and resonance effects. 3D deterministic approaches may rather enhance the ground motion simulations in the sedimentary basins and the presence of topography (Pitarka et al., 2022). Moreover, many studies present clear evidence of rupture directivity in the $M_w$6.0 2016 Amatrice earthquake, Central Italy (e.g., Calderoni et al., 2017; Ren et al., 2017). However, the method used in this study, the stochastic model through the EXSIM code (Motazedian and Atkinson, 2005; Boore, 2009), is a simplistic model and dismisses the effect of rupture directivity on the azimuth dependent variability of ground motions. Although our simulations at low frequencies capture the directivity effect (Tinti et al., 2016; Pischiutta et al., 2021), the goodness of fit between simulated and observed ground motions might be biased for some of the stations only at higher frequencies.
# DATA AVAILABILITY STATEMENT

Publicly available datasets were analyzed in this study. This data can be found at: https://itaca.mi.ingv.it/ItacaNet_32/#/home.

# AUTHOR CONTRIBUTIONS

MP, AA, and FP contributed to the conception and design of the study. MP and CF organized the database. CF and PM furnished the site amplification curves adopted in this study.

# REFERENCES

Aki, K., and Chouet, B. (1975). Origin of Coda Waves: Source, Attenuation, and Scattering Effects. *J. Geophys. Res.* 80, 3322–3342. doi:10.1029/jb080i023p03322

Akinci, A., Aochi, H., Herrero, A., Pischiutta, M., and Karanikas, D. (2017). Physics-Based Broadband Ground-Motion Simulations for Probable Mw≥7.0 Earthquakes in the Marmara Sea Region (Turkey). *Bull. Seismol. Soc. Am.* 107 (3), 1307–1323. doi:10.1785/0120160096

Akinci, A., Cheloni, D., and Dindar, A. A. (2021). The 30 October 2020, M7.0 Samos Island (Eastern Aegean Sea) Earthquake: Effects of Source Rupture, Path and Local-Site Conditions on the Observed and Simulated Ground Motions. *Bull. Earthq. Eng.* 19, 4745–4771. doi:10.1007/s10518-021-01146-5

Anderson, J. G., and Hough, S. E. (1984). A Model for the Shape of the Fourier Amplitude Spectrum of Acceleration at High Frequencies. *Bull. Seismol. Soc. Am.* 74 (5), 1969–1993. doi:10.1785/BSSA0740051969

Andrews, D. J. (1986). Objective Determination of Source Parameters and Similarity of Earthquakes of Different Size. *Earthq. Source Mech.* 37, 259–267. doi:10.1029/GM037p0259

Ben-Zion, Y., and Sammis, C. G. (2003). Characterization of Fault Zones. *Pure Appl. Geophys.* 160, 677–715. doi:10.1007/pl00012554

Bindi, D., Pacor, F., Luzi, L., Puglia, R., Massa, M., Ameri, G., et al. (2011). Ground Motion Prediction Equations Derived from the Italian Strong Motion Database. *Bull. Earthq. Eng.* 9 (6), 1899–1920. doi:10.1007/s10518-011-9313-z

Bindi, D., Parolai, S., Cara, F., Di Giulio, G., Ferretti, G., Luzi, L., et al. (2009). Site Amplifications Observed in the Gubbio Basin, Central Italy: Hints for Lateral Propagation Effects. *Bull. Seismol. Soc. Am.* 99 (2A), 741–760. doi:10.1785/0120080238

Bonamassa, O., and Vidale, J. E. (1991). Directional Site Resonances Observed from Aftershocks of the 18 October Loma Prieta Earthquake. *Bull. Seismol. Soc. Am.* 81 (5), 1945–1957. doi:10.1785/BSSA0810051945

Boore, D. M. (2005). SMSIM—Fortran Programs for Simulating Ground Motions from earthquakes: Version 2.3—A Revision of OFR 96-80-A. Available at: http://www.daveboore.com/smsim.

Boore, D. M. (2009). Comparing Stochastic Point-Source and Finite-Source Ground-Motion Simulations: SMSIM and EXSIM. *Bull. Seismol. Soc. Am.* 99, 3202–3216. doi:10.1785/0120090056

Boore, D. M. (2016). Determining Generic Velocity and Density Models for Crustal Amplification Calculations, with an Update of the Boore and Joyner (1997) Generic Site Amplification for VS(30) = 760 M/s. *Bull. Seismol. Soc. Am.* 106, 316–320. doi:10.1785/0120150229

Boore, D. M., Joyner, W. B., and Fumal, T. E. (1994). Estimation of Response Spectra and Peak Accelerations from Western North American Earthquakes: an Interim Report, Part 2. *U.S. Geol. Surv. Open-File rept.* 94–127, 40. doi:10.3133/of94127

Boore, D. M., and Joyner, W. B. (1997). Site Amplifications for Generic Rock Sites. *Bull. Seismol. Soc. Am.* 87, 327–341. doi:10.1785/bssa0870020327

Boore, D. M. (2003). Simulation of Ground Motion Using the Stochastic Method. *Pure Appl. Geophys.* 160, 635–676. doi:10.1007/pl00012553

MP and AA wrote the first draft of the manuscript. FP wrote sections of the manuscript. All authors contributed to manuscript revision and read and approved the submitted version.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2022.886606/full#supplementary-material

Boore, D. M. (2013). The Uses and Limitations of the Square-Root-Impedance Method for Computing Site Amplification. *Bull. Seismol. Soc. Am.* 103, 2356–2368. doi:10.1785/0120120283

Boore, D. M., Thompson, E. M., and Cadet, H. (2011). Regional Correlations of VS30 and Velocities Averaged over Depths Less Than and Greater Than 30 Meters. *Bull. Seismol. Soc. Am.* 101, 3046–3059. doi:10.1785/0120110071

Burjánek, J., Edwards, B., and Fäh, D. (2014a). Empirical Evidence of Local Seismic Effects At Sites With Pronounced Topography: A Systematic Approach *Geophys. J. Intern.* 197 (1), 608–619. doi:10.1093/gji/ggu014

Burjanek, J., Fäh, D., Pischiutta, M., Rovelli, A., Calderoni, G., and Bard, P.-Y. (2014b). "NERA-JRAI Working Group," in *Site Effects at Sites with Pronounced Topography: Overview & Recommendations* (Research report for EU project NERA), 64.

Burjánek, J., Moore, J. R., Yugsi Molina, F. X., and Fäh, D. (2012). Instrumental Evidence of Normal Mode Rock Slope Vibration. *Geophys. J. Int.* 188 (2), 559–569. doi:10.1111/j.1365-246x.2011.05272.x

Calderoni, G., Rovelli, A., and Di Giovambattista, R. (2017). Rupture Directivity of the Strongest 2016-2017 Central Italy Earthquakes. *J. Geophys. Res. Solid Earth* 122, 9118–9131. doi:10.1002/2017jb014118

Camassi, R., Azzaro, R., and Tertulliani, A. (2008). Macroseismology: the Lessons Learnt from the 1997/1998 Colfiorito Seismic Sequence. *Ann. Geophys.* 51 (2/3). doi:10.4401/ag-4453

Campbell, K. W., and Boore, D. M. (2016). Evaluation of Six NEHRP B/C Crustal Amplification Models Proposed for Use in Western North America. *Bull. Seismol. Soc. Am.* 106, 673–686. doi:10.1785/0120150242

Castro, R. R., Anderson, J. G., and Singh, S. K. (1990). Site Response, Attenuation and Source Spectra of S Waves along the Guerrero, Mexico, Subduction Zone. *Bull. Seismol. Soc. Am.* 80, 1481–1503. doi:10.1785/BSSA08006A1481

Chiara, F., Giovanni, L., Maria, D. A., Rodolfo, P., Lucia, L., and Francesca, P. (2018). Ground Motion Model for Reference Rock Sites in Italy. *Soil Dyn. Earthq. Eng.* 110, 276–283. doi:10.1016/j.soildyn.2018.01.024

Cornou, C., and Bard, P. Y. (2003). Site-to-bedrock over 1D Transfer Function Ratio: An Indicator of the Proportion of Edge-Generated Surface Waves? *Geophys. Res. Lett.* 30, 1453–1457. doi:10.1029/2002gl016593

Del Gaudio, C., De Risi, M. T., Ricci, P., and Verderame, G. M. (2019). Empirical Drift-Fragility Functions and Loss Estimation for Infills Reinforced Concrete Frames under Seismic Loading. *Bull. Earthq. Eng.* 17, 1285–1330. doi:10.1007/s10518-018-0501-y

Douglas, J., and Aochi, H. (2008). A Survey of Techniques for Predicting Earthquake Ground Motions for Engineering Purposes. *Surv. Geophys* 29, 187–220. doi:10.1007/s10712-008-9046-y

Falsaperla, S., Cara, F., Rovelli, A., Neri, M., Behncke, B., and Acocella, V. (2010). Effects of the 1989 Fracture System in the Dynamics of the Upper SE Flank of Etna Revealed by Volcanic Tremor Data: the Missing Link? *J. Geophys. Res.* 115 (B11306). doi:10.1029/2010JB007529

Field, E., and Jacob, K. (1995). A Comparison and Test of Various Site-Response Estimation Techniques, Including Three that Are Not Reference-Site Dependent. *Bull. Seismol. Soc. Am.* 85 (4), 1127–1143. doi:10.1785/BSSA0850041127

Galli, P., Peronace, E., Bramerini, F., Castenetto, S., Naso, G., Cassone, F., et al. (2016b). The MCS Intensity Distribution of the Devastating 24 August
2016 Earthquake in Central Italy (MW 6.2). *Ann. Geophys.* 59. doi:10.4401/ag-7287

Galli, P., Peronace, E., and Tertulliani, A. (2016a). Rapporto sugli effetti macrosismici del terremoto del 24 agosto 2016 di Amatrice in scala MCS, Roma, Rapporto congiunto DPC, CNRIGAG, 15. doi:10.5281/zenodo.161323

Géli, L., Bard, P.-Y., and Jullien, B. (1988). The Effect of Topography on Earthquake Ground Motion: a Review and New Results. *Bull. Seism. Soc. Am.* 78 (1), 42–63. doi:10.1785/BSSA0780010042

Goulet, C. A., Abrahamson, N. A., Somerville, P. G., and Wooddell, K. E. (2015). The SCEC Broadband Platform Validation Exercise: Methodology for Code Validation in the Context of Seismic-Hazard Analyses. *Seismol. Res. Lett.* 86 (1), 17–26. doi:10.1785/0220140104

Graves, R. W. (1995). Preliminary Analysis of Long-Period Basin Response in the Los Angeles Region from the 1994 Northridge Earthquake. *Geophys. Res. Lett.* 22 (2), 101–104. doi:10.1029/94gl02894

Iwata, T., Hatayama, K., Kawase, H., and Irikura, K. (1996). Site Amplification of Ground Motions during Aftershocks of the 1995 Hyogo-Ken Nanbu Earthquake in Severely Damaged Zone. *J. Phys. Earth* 44 (5), 553–561. doi:10.4294/jpe1952.44.553

Joyner, W., and Fumal, T. (1984). Use of Measured Shear-Wave Velocity for Predicting Geologic Site Effects on Strong Ground Motion. *Proc. 8th world conf. Earthq. Eng.* 2, 777–783.

Kaiser, A., Massey, C., Pischiutta, M., Fry, B., and Nicol, A. (2022). *Quantifying Seismic Amplification on Topography in New Zealand and its Relationship to Landslide Occurrence: First Steps under New Zealand's Resilience Challenge Programme*. Bellevue, Washington: Sesmological Society America Annual Meeting.

Kawase, H., Nagashima, F., Nakano, K., and Mori, Y. (2019). Direct Evaluation of S-Wave Amplification Factors from Microtremor H/V Ratios: Double Empirical Corrections to "Nakamura" Method. *Soil Dyn. Earthq. Eng.* 126, 105067. doi:10.1016/j.soildyn.2018.01.049

Lachet, C., and Bard, P.-Y. (1994). Numerical and Theoretical Investigations on the Possibilities and Limitations of Nakamura's Technique. *J. Phys. Earth* 42 (5), 377–397. doi:10.4294/jpe1952.42.377

Lanzano, G., Felicetta, C., Pacor, F., Spallarossa, D., and Traversa, P. (2022). Generic-to-reference Rocks Scaling Factors for the Seismic Ground Motion in Italy. *Bull. Seismol. Soc. Am.* (accepted), 1–24. doi:10.1785/0120210063

Lanzano, G., Felicetta, C., Pacor, F., Spallarossa, D., and Traversa, P. (2020). Methodology to Identify the Reference Rock Sites in Regions of Medium-To-High Seismicity: an Application in Central Italy. *Geophys. J. Int.* 222 (3), 2053–2067. doi:10.1093/gji/ggaa261

Lewis, M. A., and Ben-Zion, Y. (2010). Diversity of Fault Zone Damage and Trapping Structures in the Parkfield Section of the San Andreas Fault from Comprehensive Analysis of Near Fault Seismograms. *Geophys. J. Int.* 183, 1579–1595. doi:10.1111/j.1365-246x.2010.04816.x

Licia Faenza, L., Valentino Lauciani, V., and Alberto Michelini, A. (2016). The ShakeMaps of the Amatrice, M6, Earthquake. *Ann. Geophys.* 59. FAST TRACK 5. doi:10.4401/ag-7238

Mai, P. M., and Beroza, G. C. (2003). A Hybrid Method For Calculating Near Source, Broadband Seismograms: Application To Strong Motion Prediction *Phys. Earth Planet. Inter.* 137 (1/4), 183–199. doi:10.1016/S0031-9201(03)00014-1

Malagnini, L., Akinci, A., Mayeda, K., Munafo', I., Herrmann, R. B., and Mercuri, A. (2011). Characterization of Earthquake-Induced Ground Motion from the L'Aquila Seismic Sequence of 2009, Italy. *Geophys J. Int.* 184 (1), 325–337. doi:10.1111/j.1365-246x.2010.04837.x

Malagnini, L., and Munafo', I. (2018). On the Relationship between $M_L$ and $M_W$ in a Broad Range: An Example from the Apennines, Italy. *Bull. Seism. Soc. Am.* 108 (2), 1018–1024. doi:10.1785/0120170303

Zanini, M. A., Lorenzo Hofer, L., Flora Faleschini, F., Paolo Zampieri, P., Nicola Fabris, N., and Carlo Pellegrino, C. (2016). Preliminary Macroseismic Survey of the 2016 Amatrice Seismic Sequence. *Ann. Geophys.* 59. Fast Track 5. doi:10.4401/ag-7172

Martino, S., Minutolo, A., Paciello, A., Rovelli, A., Mugnozza, G. S., and Verrubbi, V. (2006). Evidence of Amplification Effects in Fault Zone Related to Rock Mass Jointing. *Nat. Hazards* 39, 419–449. doi:10.1007/s11069-006-0001-2

Marzorati, S., Ladina, C., Falcucci, E., Gori, S., Saroli, M., Ameri, G., et al. (2011). Site Effects "on the Rock": the Case of Castelvecchio Subequo (L'Aquila, Central Italy). *Bull. Earthq. Eng.* 9, 841–868. doi:10.1007/s10518-011-9263-5

Michelini, A., Faenza, L., Lauciani, V., and Malagnini, L. (2008). Shakemap Implementation in Italy. *Seismol. Res. Lett.* 79, 688–697. doi:10.1785/gssrl.79.5.688

Moore, J. R., Gischig, V., Burjanek, J., Loew, S., and Fah, D. (2011). Site Effects in Unstable Rock Slopes: Dynamic Behavior of the Randa Instability (Switzerland). *Bull. Seismol. Soc. Am.* 101 (6), 3110–3116. doi:10.1785/0120110127

Molnar, S., Cassidy, J. F., Castellaro, S., Cornou, H., Crow, J. A., Hunter, S., et al. (2018). Application of Microtremor Horizontal-to-Vertical Spectral Ratio (MHVSR) Analysis for Site Characterization: State of the Art *Surv. Geophys.* 39, 613–631. doi:10.1007/s10712-018-9464-4

Morasca, P., D'Amico, M., Lanzano, G., Sgobba, S., Colavitti, L., Pacor, F., et al. (2022). Empirical Correlations between a FAS Non-ergodic Ground Motion Model and a GIT Derived Model for Central Italy. *Geophys. J. Int.* (in prep.).

Motazedian, D., and Atkinson, G. M. (2005). Stochastic Finite-Fault Modeling Based on a Dynamic Corner Frequency. *Bull. Seismol. Soc. Am.* 95, 995–1010. doi:10.1785/0120030207

NEHRP Building Seismic Safety Council (2015). "NEHRP Recommended Seismic Provisions for New Buildings and Other Structures," in *Part 1 Provisions, Part 2 Commentary, FEMA Rept.* 2015 Edition (Washington, D.C: National Institute of Building Sciences), 1, P-1050-1051.

NTC18 (2018). Italian building code 'Norme Tecniche per le Costruzioni - Ministero delle Instrutture e dei Trasporti. Available at: https://www.gazzettaufficiale.it/eli/gu/2018/02/20/42/so/8/sg/pdf. (Accessed December 2020)

NZS (2004). *Standards New Zealand 2004. Structural Design Actions—Part 5 Earthquake Actions—New Zealand*. Standards New Zealand. Private Bag 2439, Wellington 6020. ISBN 1-86975-018-7 (incorporating Amendment 1 (following Canterbury Earthquake Sequence). Available at: https://www.standards.govt.nz/shop/nzs-1170-52004/.

Ojeda, J., Akinci, A., Tinti, E., Arriola, S., and Ruiz, S. (2021). Hybrid Broadband Strong-Motion Simulation to Investigate the Near-Source Characteristics of the M6.5, 30 October 2016 Norcia, Italy Earthquake. *Soil Dyn. Earthq. Eng.* 149, 106866. doi:10.1016/j.soildyn.2021.106866

Panzera, F., Pischiutta, M., Lombardo, G., Monaco, C., and Rovelli, A. (2014). Wavefield Polarization in Fault Zones of the Western Flank of Mt. Etna: Observations and Fracture Orientation Modelling. *Pure Appl. Geophys.* 171, 3083–3097. doi:10.1007/s00024-014-0831-x

Pischiutta, M., Akinci, A., Tinti, E., and Herrero, A. (2020). Broad-band Ground-Motion Simulation of 2016 Amatrice Earthquake, Central Italy. *Cent. Italy. Geophys. J. Int.* 224, 1753–1779. doi:10.1093/gji/ggaa412

Pischiutta, M., Cianfarra, P., Salvini, F., Cara, F., and Vannoli, P. (2018). A Systematic Analysis of Directional Site Effects at Stations of the Italian Seismic Network to Test the Role of Local Topography. *Geophys. J. Int.* 214 (1), 635–650. doi:10.1093/gji/ggy133

Pischiutta, M., Fondriest, M., Demurtas, M., Magnoni, F., Di Toro, G., and Rovelli, A. (2017). Structural Control on the Directional Amplification of Seismic Noise (Campo Imperatore, Central Italy). *Earth Planet. Sci. Lett.* 471, 10–18. doi:10.1016/j.epsl.2017.04.017

Pischiutta, M., Salvini, F., Fletcher, J., Rovelli, A., and Ben-Zion, Y. (2012). Horizontal Polarization of Ground Motion in the Hayward Fault Zone at Fremont, California: Dominant Fault-High-Angle Polarization and Fault-Induced Cracks. *Geophys. J. Int.* 188 (3), 1255–1272. doi:10.1111/j.1365-246x.2011.05319.x

Pischiutta, M., Savage, M. K., Holt, R. A., and Salvini, F. (2015). Fracture-related Wavefield Polarization and Seismic Anisotropy across the Greendale Fault. *J. Geophys. Res. Solid Earth* 120, 7048–7067. doi:10.1002/2014jb011560

Pitarka, A., Akinci, A., De Gori, P., and Buttinelli, M. (2022). Deterministic 3D Ground-Motion Simulations (0-5 Hz) and Surface Topography Effects of the 30 October 2016 Mw 6.5 Norcia, Italy, Earthquake. *Bull. Seismol. Soc. Am.* 112 (1), 262–286. doi:10.1785/0120210133

Priolo, E., Pacor, F., Spallarossa, D., Milana, G., Laurenzano, G., Romano, M. A., et al. (2019). Seismological Analyses of the Seismic Microzonation of
138 Municipalities Damaged by the 2016–2017 Seismic Sequence in Central Italy. *Bull. Earthq. Eng.* 18 (12), 1–41. doi:10.1007/s10518-019-00652-x

Quest, W. G., Raffaele, A., Andrea, T., Filippo, B., Romano, C., Sergio Del, M., et al. (2016). The 24 August 2016 Amatrice Earthquake: Macroseismic Survey in the Damage Area and EMS Intensity Assessment. *Ann. Geophys.* 59. fast track 5. doi:10.4401/ag-7203

Ren, Y., Wang, H., and Wen, R. (2017). Imprint of Rupture Directivity from Ground Motions of the 24 August 2016 M W 6.2 Central Italy Earthquake. *Tectonics* 36, 3178–3191. doi:10.1002/2017tc004673

Russo, E., Felicetta, C., D'Amico, M., Sgobba, S., Lanzano, G., Mascandola, C., et al. (2022). "Italian Accelerometric Archive v3.2 - Istituto Nazionale di Geofisica e Vulcanologia," in *Dipartimento Della Protezione Civile Nazionale*. doi:10.13127/itaca.3.2

Sánchez-Sesma, F., Chávez-Pérez, S., Suàrez, M., Bravo, M. A., and Pérez-Rocha, L. E. (1988). The Mexico Earthquake of September 19, 1985—On the Seismic Response of the Valley of Mexico. *Earthq. Spectra* 4 (3), 569–589. doi:10.1193/1.1585491

Sextos, A., De Risi, R., Pagliaroli, A., Foti, S., Passeri, F., Ausilio, E., et al. (2018). Local Site Effects and Incremental Damage of Buildings during the 2016 Central Italy Earthquake Sequence. *Earthq. Spectra* 34 (4), 1639–1669. doi:10.1193/100317eqs194m

Sgobba, S., Lanzano, G., and Pacor, F. (2020). Empirical Nonergodic Shaking Scenarios Based on Spatial Correlation Models: An Application to Central Italy. *Earthq. Engng Struct. Dyn.* 50, 60–80. doi:10.1002/eqe.3362

Tinti, E., Scognamiglio, L., Michelini, A., and Cocco, M. (2016). Slip Heterogeneity and Directivity of the ML 6.0, 2016, Amatrice Earthquake Estimated with Rapid Finite-Fault Inversion. *Geophys. Res. Lett.* 43, 10745–10752. doi:10.1002/2016gl071263

Wald, D. J., Quitoriano, V., Heaton, T. H., and Kanamori, H. (1999). Relationships between Peak Ground Acceleration, Peak Ground Velocity, and Modified Mercalli Intensity in California. *Earthq. Spectra* 15 (3), 557–564. doi:10.1193/1.1586058

Zhu, C., Cotton, F., Kawase, H., Haendel, A., and Nakano, K. (2022b). How Well Can We Predict Earthquake Site Response So Far? Machine Learning vs. *Physics-Based Modeling*. doi:10.13140/RG.2.2.34933.91365

Zhu, C., Cotton, F., Kawase, H., Haendel, A., Pilz, M., and Nakano, K. (2022a). How Well Can We Predict Earthquake Site Response So Far? Site-specific Approaches. *Earthq. Spectra* 38 (2), 1047–1075. doi:10.1177/87552930211060859

Zhu, C., Pilz, M., and Cotton, F. (2020). Which Is a Better Proxy, Site Period or Depth to Bedrock, in Modelling Linear Site Response in Addition to the Average Shear-Wave Velocity? *Bull. Earthq. Eng.* 18 (3), 797–820. doi:10.1007/s10518-019-00738-6

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

**Publisher's Note:** All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Copyright © 2022 Pischiutta, Akinci, Felicetta, Pacor and Morasca. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
<figure><img src="image_66.png" /></figure>

Check for
updates

# Seismo-Stratigraphic Model for the Urban Area of Milan (Italy) by Ambient-Vibration Monitoring and Implications for Seismic Site Effects Assessment

M. Massa*, S. Lovati, R. Puglia, G. Brunelli, A. Lorenzetti, C. Mascandola, C. Felicetta, F. Pacor and L. Luzi¹

¹National Institute of Geophysics and Volcanology (INGV), Milano, Italy

In this paper, we present the work carried out to characterize the spatial variability of seismic site response related to local soil conditions in the city of Milan and its surroundings, an area with ~3 million inhabitants and a high density of industrial facilities. The area is located at the northwestern end of the Po Plain, a large and deep sedimentary basin in northern Italy. An urban-scale seismo-stratigraphic model is developed based on new passive and active seismic data, supported by the available geological data and stratigraphic information from shallow and deep vertical wells. In particular, 33 single-station and 4 ambient-vibration array measurements are acquired, together with 4 active multichannel analyses of surface waves (MAWS). To estimate the resonant frequencies of the sediments, the horizontal to vertical spectral ratio technique (HVSR) is applied to the ambient-vibration recordings, whereas to determine the Rayleigh-wave dispersion curves from the passive array, the data are analysed using the conventional frequency-wavenumber, the modified spatial autocorrelation and the extended spatial autocorrelation (ESAC) techniques. The array data are used to determine the local shear wave velocity profiles, $V_S$, via joint inversion of the Rayleigh-wave dispersion and ellipticity curves deduced from the HVSR. The results from HVSR show three main bands of amplified frequencies, the first in the range 0.17–0.23 Hz, the second from 0.45 to 0.65 Hz and the third from 3 to 8 Hz. A decreasing trend of the main peaks is observed from the northern to the southern part of the city, allowing us to hypothesize a progressive deepening of the relative regional chrono-stratigraphic unconformities. The passive ambient noise array and MASW highlight the dispersion of the fundamental mode of the Rayleigh-wave in the range 0.4–30 Hz, enabling to obtain detailed Vs. profiles with depth down to about 1.8 km. The seismo-stratigraphic model is used as input for 1D numerical modelling assuming linear soil conditions. The theoretical 1D transfer functions are compared to the HVSR curves evaluated from both ambient noise signals and earthquake waveforms recorded by the IV. MILN station in the last 10 years.

**OPEN ACCESS**

**Edited by:**
Roberto Paolucci,
Politecnico di Milano, Italy

**Reviewed by:**
Alessandro Vuan,
Istituto Nazionale di Oceanografia e di
Geofisica Sperimentale, Italy
Maria Rosaria Gallipoli,
Institute of Methodologies for
Environmental Analysis (IMAA), Italy

**\*Correspondence:**
M. Massa
marco.massa@ingv.it

**Specialty section:**
This article was submitted to Geohazards and Georisks, a section of the journal Frontiers in Earth Science

**Received:** 07 April 2022

**Accepted:** 19 May 2022

**Published:** 14 June 2022

**Citation:**
Massa M, Lovati S, Puglia R, Brunelli G, Lorenzetti A, Mascandola C, Felicetta C, Pacor F and Luzi L (2022) Seismo-Stratigraphic Model for the Urban Area of Milan (Italy) by Ambient-Vibration Monitoring and Implications for Seismic Site Effects Assessment. Front. Earth Sci. 10:915083
doi: 10.3389/feart.2022.915083
**Keywords: milan (Italy), seismic site characterisation, seismostratigraphic model, vs. velocity profile, seismic array, MASW, HVSR (horizontal to vertical spectral ratio), 1D numerical modeling**

## INTRODUCTION

It is widely recognized that a significant proportion of the variability of earthquake ground motion is related to local geological conditions, which can modify the ground-motion amplitude, duration, and frequency. In alluvial basins (or valleys), layers can cause resonance of vertically propagating shear waves at specific frequencies and/or trap surface waves generated locally at the edges of the basin (e.g., Hanks, 1975; Hisada et al., 1993; Bielak et al., 1999; Sato et al., 1999; Joyner, 2000; Kagawa et al., 2004; Somerville et al., 2004; D'Amico et al., 2006; D'Amico et al., 2008; De Ferrari et al., 2010; Ronald Abraham et al., 2015).

In this framework and in particular for urban planning (e.g. Ansal et al., 2009; Ansal et al., 2010, Crespellani 2014, Celikbilek and Sapmaz 2016, Aversa and Crespellani 2016, etc.), the seismic

site effect estimation plays a crucial role in earthquake risk reduction.

The importance of site characterization of urban or metropolitan areas is recognized worldwide, with several studies performed in New York (Stephenson et al., 2009), Taipei (Wang, 2008), Perth (Liang et al., 2009), Seoul (Sun et al., 2014; Kim et al., 2017), Beijing (Liu et al., 2014), Vancouver (Molnar et al., 2020), San José, Costa Rica (Moya et al., 2000), Istanbul (Picozzi et al., 2009; Ansal et al., 2010), Memphis (Schneider et al., 2001), Bucharest (Ritter et al., 2005), among others.

Similar studies in Italy are performed for other cities like Matera (Tragni et al., 2021), Rome (Marcucci et al., 2019), and Firenze (D'Amico et al., 2008), as well as small villages like Mirandola (Tarabusi and Caputo, 2017) and Cavezzo (Lai et al., 2020), both located in the Po plain and strongly struck

<figure><img src="image_67.png" /></figure>

**FIGURE 1** | Panel (A) simplified structural sketch of the Po plain (from Pieri and Groppi, 1981; Bigi et al., 1990). The red square indicates the urban area of Milan. White circles indicate the seismicity from yrs. 1000 to 2022. Dark gray and yellow circles drawn in the red inset indicate the recent (after 1 January 1985, http://terremoti.ingv.it/) and the historical seismicity (yrs. 1000–1985, https://emidius.mi.ingv.it/CPTI15-DBMI15/; Rovida et al., 2020), respectively, occurred in an area of 100 km² around the MMA. The composite seismogenic sources are derived from DISS (Database of Individual Seismogenis Sources, http://diss.ingv.it). Panel (B) Seismic hazard maps of Italy (MPS Working Group, 2004; Stucchi et al., 2011), in terms of peak ground acceleration (PGA) with 10% possibility of exceedance in 50 years (i.e., a return period of 475 years), referred to the rigid soils (i.e., Vs. > 800 m/s; NTC 2018). White circles indicate the epicenters of the earthquakes used for the HVSR analyses (see also **Table 1**).
<figure><img src="image_68.png" /></figure>

**FIGURE 2** | Panel (A) 1:10000 Geological map from the 118 Milan sheet (ISPRA et al., 2016). Po synthem (POI, Upper Pleistocene-Holocene), Cantù synthem (LCN, Upper Pleistocene), Ronchetto delle Rane subsynthem (LCN₄, Upper Pleistocene), Bulgarograsso unit (BXE, Middle-Upper Pleistocene), Minoprio Unit (BMI, Middle-Upper Pleistocene) and Guanzate unit (BEZ, Middle-Upper Pleistocene) are the stratigraphic units characterizing the MMA (ISPRA et al., 2016). The dotted black line indicates the NNW-SSE geological section showed in panel (B). Panel (B) geological section from North Park (NNW) to Vettabbia park (SSE), with thickness of about 200 m. Yellow line represents the Qc3 Unconformity from GeoMol Team, 2015 (or Y-Surface from Scardia et al., 2012), while red line represents the Qc1 Unconformity from GeoMol Team, 2015 (or R-Surface from Muttoni et al., 2003). LS and LI indicate the upper and lower Padano supersynthem, while PD indicates the Padano supersynthem (ISPRA et al., 2016). Panel (C) geological section from North Park (NNW) to Vettabbia park (SSE), with thickness of about 1800 m. MESa (Burdigalian-Tortonian), MESb (Messinian), PL (early-middle Pliocene) and PLMa (Middle lower Pleistocene) indicate the stratigraphic units as recognized by GeoMol project (2015), while PL (base of Pliocene, 5,3 Myrs), GEL (base of Pleistocene, 2,6 Myrs), Qm1 (1,5 Myrs) and Qm2 (1,25 Myrs) represent the related unconformities. Dotted lines indicate the deduced depths on the basis of the VIDEPI project (2009) deep wells.

by the 2012 Emilia seismic sequence (Luzi et al., 2013). Worth mentioning is the intensive studies that have been performed after the 2016–2017 Central Italy seismic sequence (Michele et al., 2016) to obtain detailed microzonation maps in the 138 most damaged municipalities (Hailemikael et al., 2020).

In this paper, we focus on the city of Milan, with an area of about 1.570 km² and a population of more than 3 million (http://

www.istat.it). Milan and its hinterlands represent the area with the highest density of industrial plants, service industries, skyscrapers and high-speed railways in Italy (http://www.istat.it). In particular, Milan is the Italian city that has grown more vertically through the construction of skyscrapers, at present with maximum height of ~230 m and 50 stories (https://en.wikipedia.org/wiki/List_of_tallest_buildings_in_Milan). Recently, a couple
of 5s seismometers, with reserved data, is installed at the floors -3
(i.e., -15 m from surface) and 37 (height ~140 m), respectively, of
the skyscraper Lombardia-building (https://en.wikipedia.org/
wiki/Palazzo_Lombardia), allowing us to preliminary recognize
a natural periods of vibration of ~3.2 s.

The MMA (northern Italy, **Figure 1**, panel A) is located in the northwestern part of the Po Plain, one of the largest alluvial basins worldwide. The Po alluvial basin covers an area of about 50.000 km², with alluvial deposits that can reach a depth of several kilometers and that directly overlie the deep Miocene geological bedrock (Pieri and Groppi, 1981; Regione Lombardia, Eni Divisione Agip, 2002).

The bottom panel of **Figure 1** (panel B) shows the reference Italian seismic hazard map of Italy in terms of peak ground acceleration for rock-site conditions with 10% probability of exceedance in 50 years (MPS Working Group, 2004; Stucchi et al., 2011): despite the low-to-medium seismic hazard of the area, the exposure level of the city, the local geological condition (**Figure 2**) and the proximity with active buried seismogenic structure (**Figure 1**, panel A) make the MMA a medium to high seismic risk zone.

The studied area is, in general, characterized by a low rate of
annual seismicity (**Figure 1**, panel A). In particular, in the last
40 years no local earthquakes with magnitude >4.5 have occurred.
However, the area has often undergone ground motion over the
long-period, namely in the case of distant earthquakes with
higher magnitude, such as the 2012 Emilia, Mw 5.8, seismic
sequence (Luzi et al., 2013).

The importance of the long-period (>1 s) component of seismic ground motion has been recognized worldwide during some strong earthquakes (e.g., 1985, Mw 8.1, Michoacán, Furumura and Kennett, 1998; Furumura and Kennett, 1999, Mw 7.1, Hector Mine, Grazier et al., 2002, Mw 7.9, Denali, Ratchkovski et al., 2003, Mw 8.0, Tokachi-Oki, Miyazaki et al., 2004). Significant damage to long-period structures at large hypocentral distances, as well as prolonged duration of seismic ground motion, has prompted seismologists to focus in understanding the seismic wave propagation and amplification effects in deep sedimentary basins. The Mexico City records of the 1985 Michoacan earthquake lead to far source long-period ground motions being known around the world, with long-period components between 2 and 4 s (Beck and Hall, 1986). Before the 1985 Michoacan earthquake, Japanese seismologists observed far-source long-period ground motion in the Niigata basin, during the 1964 Niigata, 1968 Tokachi-oki, and 1983 Japan Sea earthquakes (Koketsu and Miyake, 2008). Subsequently, the 2003 Tokachi-oki earthquake provided the southern coast of Hokkaido in Japan, with one of the most significant examples (Koketsu et al., 2005; Hatayama, 2008) of far-source long-period ground motion, with dominant amplitudes in the range 7–8 s and displacements of a few meters (Koketsu and Miyake, 2008).

Long-period ground motion amplification is a significant issue
in the Po Plain. In this area, long-period resonance phenomena
have been well-documented by many studies (e.g., Priolo et al.,
2012; Luzi et al., 2013; Martelli and Romani, 2013; Massa and
Augliera, 2013; Milana et al., 2013; Paolucci et al., 2015;
Laurenzano et al., 2017; Mascandola et al., 2017; Tarabusi and

Caputo, 2017). However, most of these are performed in the
epicentral area of the 2012 Emilia seismic sequence (e.g., Luzi
et al., 2013), during which significant ground-motion
amplifications and increased duration are observed. A recent
study by Mascandola et al. (2019) widely mapped 1–3 s long-
period resonance phenomena in the Po Plain, through extensive
ambient-vibration measurements performed in the basin.

This study attempts to investigate the influence of the shallow
to deep (down ~1.8 km) alluvial deposits for the metropolitan
area of Milan, where amplification effects in the frequency range
0.2–8 Hz are observed. Milan has the disadvantage of having few
earthquake recordings due to the relatively low seismicity and
distant reference sites in the Alps or in the Apennine Mountains
(Figure 1, panel A). These features prevent applying the standard
spectral ratio approach (SSR, Bocherdt, 1970) to compute the
empirical soil amplification functions and lead to adopting
techniques based on numerical modelling, which are more
accurate the more detailed the seismo-stratigraphic model is.

Therefore, in this study a detailed seismo-stratigraphic model
obtained from the analysis of existing geological and stratigraphic
data (e.g., Regione Lombardia, Eni Divisione Agip, 2002; Muttoni
et al., 2003; De Franco et al., 2009; ViDEPI project, 2009; Scardia
et al., 2012; GeoMol Team, 2015; ISPRA et al., 2016) and newly
acquired geophysical data is proposed.

Different active and passive seismic measurements are carried
out to characterize the subsoil. The passive analyses represent the
main part of the study, consisting of 33 single-station and 4 array
measurements of ambient noise (Figure 2, panel A). The single-
station measurements allowed defining the resonance frequency
of the site through the Nakamura (1989) horizontal-to-vertical
spectral ratio (HVSR) method. The microtremor arrays,
performed in correspondence of three selected sites (i.e., North
Park, Giuriati sports center, Vettabbia park) located along the
alignment reported in Figure 2 (panel A), enabled assessing vs.
profiles down to depths of about ~1.8 km.

The active measurements consisted of 4 Multichannel
Analyses of Surface Waves (MASW, Park et al., 1999) at
North and Vettabbia public parks (Figure 2, panel A) which
are useful to define the soil categories, according to the Italian
(Nuove Norme Tecniche per le Costruzioni, NTC, 2018) and the
European (Eurocode 8, CEN 2004) seismic codes for buildings.
Soil categories are defined according to the values of the average
seismic shear-wave velocity from the surface to a depth of 30 m
(Vs,30).

The results obtained from ambient noise HVSR are compared
to the HVSR curves considering 17 earthquakes with magnitude
from 3.7 to 5.8 recorded by the IV. MILN station (Figure 2, panel
B) in the last 10 years (Table 1).

The Vs profiles are finally validated through 1D numerical
models performed assuming a linear elastic behaviour of the soil
using the SHAKE91 code (Idriss and Sun, 1993).

Seismicity of the Area and Recorded
Ground Motion

In the last 30 years, thousands of small-energy seismic events
(~200/year, **Figure 1**, panel A, http://terremoti.ingv.it/) have
<table><caption>TABLE 1 | Earthquakes with magnitude ≥3.7 recorded at station IV.MILN in the last 10 years and used for the HVR analyses. The epicentres are depicted in panel B) of Figure 1. In bold are indicated magnitude values > 5.</caption><thead><tr><th>YYYY/MM/DD hh:mm:ss (UTC)</th><th>Lon (')</th><th>Lat (')</th><th>Depth (km)</th><th>Mw (ML)</th><th>Stat</th><th>REPI (km)</th></tr></thead><tbody><tr><td>2012/02/29 15:34:57</td><td>9.0170</td><td>44.8240</td><td>6.50</td><td>3.7</td><td>IV.MILN</td><td>74</td></tr><tr><td>2012/05/20 02:03:50</td><td>11.2300</td><td>44.8900</td><td>6.30</td><td>5.8</td><td>IV.MILN</td><td>169</td></tr><tr><td>2012/05/20 13:18:02</td><td>11.4900</td><td>44.8310</td><td>4.70</td><td>5.1</td><td>IV.MILN</td><td>191</td></tr><tr><td>2012/05/29 07:00:02</td><td>11.0900</td><td>44.8500</td><td>10.20</td><td>5.6</td><td>IV.MILN</td><td>161</td></tr><tr><td>2012/05/29 10:55:57</td><td>11.0080</td><td>44.8880</td><td>6.80</td><td>5.3</td><td>IV.MILN</td><td>153</td></tr><tr><td>2012/06/03 19:20:43</td><td>10.9430</td><td>44.8990</td><td>9.20</td><td>5.1</td><td>IV.MILN</td><td>148</td></tr><tr><td>2018/05/19 16:41:21</td><td>9.684</td><td>44.8160</td><td>28.30</td><td>4.3</td><td>IV.MILN</td><td>81</td></tr><tr><td>2020/04/16 09:42:51</td><td>9.4087</td><td>44.6488</td><td>3.40</td><td>4.2</td><td>IV.MILN</td><td>93</td></tr><tr><td>2020/04/19 09:53:40</td><td>9.2257</td><td>44.9760</td><td>31.80</td><td>3.7</td><td>IV.MILN</td><td>56</td></tr><tr><td>2020/12/17 15:59:22</td><td>9.1282</td><td>45.4918</td><td>56.20</td><td>3.9</td><td>IV.MILN</td><td>8</td></tr><tr><td>2020/12/29 14:36:57</td><td>11.0425</td><td>45.2442</td><td>9.40</td><td>3.9</td><td>IV.MILN</td><td>143</td></tr><tr><td>2021/07/01 11:11:49</td><td>8.3508</td><td>46.6062</td><td>5.00</td><td>3.7</td><td>IV.MILN</td><td>142</td></tr><tr><td>2021/12/18 10:34:47</td><td>9.5860</td><td>45.6752</td><td>26.30</td><td>3.9</td><td>IV.MILN</td><td>35</td></tr><tr><td>2022/02/01 22:43:07</td><td>10.2703</td><td>44.5547</td><td>47.40</td><td>3.7</td><td>IV.MILN</td><td>131</td></tr><tr><td>2022/02/09 20:00:57</td><td>10.7470</td><td>44.7887</td><td>6.40</td><td>4.3</td><td>IV.MILN</td><td>141</td></tr></tbody></table>

occurred in the Po Plain. Except of the 20 May 2012, Mw 5.8, Emilia seismic sequence (Luzi et al., 2013), with epicenters located 170 km southeast of Milan, the western part of the Po Plain area underwent a few moderate events, with magnitude between 4.5 and 5.5, such as the 14th August 2000, Mw 4.8, Asti the 3rd April 2003, Mw 4.8, Novi Ligure, the 14th September 2003, Mw 5.2, Monghidoro, the 24th November 2004, Mw 5.2, Salò and the 23rd December 2008, Mw 5.4, Parma events (http://terremoti.ingv.it). These moderate earthquakes did not cause losses of human life, but had negative economic consequences amounting to hundreds of millions of Euros (e.g., damage to buildings, schools, churches, industrial facilities etc.) highlighting the high degree of exposure of the region.

Focusing on MMA, the most powerful events recorded in the last decades occurred on 17th December 2020 (Mw 3.9 at epicentral distance of 5 km, http://terremoti.ingv.it) and on 18th December 2021 (Mw 3.9 at epicentral distance of 35 km, http://terremoti.ingv.it). In particular, in the past 10 years, the station IV. MILN (http://terremoti.ingv.it/instruments/station/MILN), now installed at INGV Milan department, recorded 17 seismic events in the magnitude (Mw) range 3.7–5.8, with epicentral distances up to 200 km (Table 1). Considering the historical seismicity (Rovida et al., 2020; https://emidius.mi.ingv.it/CPTI15-DBMI15/), the more significant earthquakes that struck the MMA occurred along the buried thrust systems reported in Figure 1 (panel A). The main seismogenic sources of the area (Figure 1, panel A) are derived from the Database of Individual Seismogenesis Sources (DISS, http://diss.ingv.it). Considering 100 km² around MMA, (Figure 1, panel A), the main historical events occurred on 26th November 1369 (Monza, Mw 5.33), on 12th May 1802 (Valle dell’Oglio, Mw 5.60), on 7th April 1786 (Po Plain, Mw 5.22) and on 15th May 1951 (Lodigiano, Mw 5.17). The highest macroseismic intensity (I = 7, https://emidius.mi.ingv.it/CPTI15-DBMI15/) for the area is associated to the 3rd January 1117, Mw 6.52, Verona earthquake, occurred ~140 km, South-East of MMA.

Figure 3 shows the ground motion values recorded by the IV. MILN station in terms of peak ground acceleration (PGA), peak ground velocity (PGV) and response spectral (Sa) ordinates at period (T) of 3.0 s. The highest values of shaking in terms of acceleration is recorded on the 18th December 2021, Mw 3.9, Dalmine event (2021/12/18 10:34:47 UTC, http://terremoti.ingv.it) with maxima horizontal and vertical PGAs of 36 gal and 13 gal, respectively. PGAs between 5 gal and 10 gal, are recorded during the 17th December 2020, Mw 3.9, Pero event (2020/12/17 15:59:22 UTC, http://terremoti.ingv.it) and the 20th and 29th May 2012 Emilia main events (2012/05/20 02:03:52 UTC, Mw 5.8 and 2012/05/29 07:00:03 UTC, Mw 5.6, http://terremoti.ingv.it). This sequence produced seismic motion characterized by the dominance of surface waves (Love and Rayleigh waves) with amplitude higher than S-phases and by the arrival of a quasi-monochromatic wave train centered at 0.2 Hz. This phenomenon, described in Luzi et al. (2013) is more evident in the central and eastern sectors of the Po Plain, even if it is also observed at IV. MILN station in the case of the 20 May 2012, Mw 5.8, Emilia mainshock.

## Geological Setting and Model

The municipality of Milan (Figure 1, panel A) is located in the north-western sector of the Po Plain, a large and deep sedimentary basin between the south-verging thrust system of the Alps and the north-verging thrust system of the Apennines (Pieri and Groppi, 1981; Bigi et al., 1990); Figure 1 (panel A) also shows the expected maximum thickness of the Po plain sediments, represented by light grey contour lines. The Plio-Quaternary deposits of the Po Plain, overlie, in discordance, the geologic bedrock composed of deep-sea turbiditic deposits, which are deformed by compressive tectonic phases that lasted until the Upper Miocene along the piedmont edge of the Southern Alps. Since the Messinian, the tectonic phases of alpine structures stopped and the northern Apennines migrated towards NE. Since that time, the deposition and the geometries of the Po plain have been constrained by rapid uplift and forward
<figure><img src="image_69.png" /><figcaption>**FIGURE 3** | Peak Ground Acceleration (PGA), Velocity (PGV) and response spectral ordinates at period of 3.0 s for all events of magnitude ≥3.5 recorded at IV. MILN station from 1st January 2012. Left and right side indicate the horizontal and the vertical ground motion values, respectively.</figcaption></figure>

movements of Apennine system, and to very long periods of isostatic subsidence of the basin.

From the available geological and geophysical information, a subsoil model with depth down to ~1.8 km is built up along a NNW-SSE section, that crosses the three sites selected for the geophysical surveys (i.e., North Park, Giuriati sports center and Vettabbia park, see **Figure 2**, panel A).

In particular, we made use of:

- the 1:10.000 sheet 118 geological map (ISPRA et al., 2016);

- 5 selected water wells down to 120 m in depth (blue circles, **Figure 2**, panel A) available at https://www.geoportale.regione.lombardia.it/;

- one core drilling down to 100 m in depth (brown circle, **Figure 2**, panel A) performed in the framework of the Carg Project (ISPRA et al., 2016);

- the deep wells of ViDEPI project (2009) surrounding the metropolitan area of Milan (black circles in **Figure 2**, panel A);

- the stratigraphic units from Carg (ISPRA et al., 2016) and GeoMol Team (2015) projects;

- the isobaths published by ISPRA et al. (2016), in the framework of the Carg project (https://www.isprambiente.gov.it/Media/carg/118_MILANO_SOTTO/Foglio.html);

- the North South high-resolution reflection seismic profile acquired in the western area of Milan by De Franco et al. (2009).

The geological map reported in **Figure 2** (panel A) show with different colours all the outcropping units of the 118-Milano geological sheet that belong to the Upper Pleistocene-Holocene continental succession. From a lithological point of view, these
<table><caption>TABLE 2 | Subsoil models at North park, Giuriati sports center and Vettabbia park.</caption><thead><tr><th>Period</th><th>Epoch</th><th>Lithology</th><th>Facies</th><th>Stratigraphic units</th><th>Unconformity</th><th>Hydrogeological units</th><th>North Park (Depth, m)</th><th>Giuriati (Depth, m)</th><th>Vettabbia (Depth, m)</th></tr></thead><tbody><tr><td>Quaternary</td><td>Upper Pleistocene-Holocene</td><td>Coarse gravel</td><td>Continental</td><td>LS</td><td>QC3 (0,45My)</td><td>Aquifer A</td><td>50</td><td>40</td><td>45</td></tr><tr><td>Quaternary</td><td>Upper Pleistocene</td><td>Fine gravel with sand and clay</td><td>Continental</td><td>LI</td><td>QC1 (0,87My)</td><td>Aquifer B</td><td>100</td><td>100</td><td>120</td></tr><tr><td>Quaternary</td><td>Upper Pleistocene</td><td>Mixture of fine gravel, sand, clay</td><td>Continental</td><td>PD</td><td>QM2 (1,25My)</td><td>Aquifer C</td><td>300</td><td>300</td><td>400</td></tr><tr><td>Quaternary</td><td>Middle lower Pleistocene</td><td>Sand (Sabbie di Asti) with % of clay</td><td>Transitional marine</td><td>PD</td><td>QM1 (1,5My)</td><td>Aquifer D</td><td>600</td><td>620</td><td>700</td></tr><tr><td>Quaternary</td><td>Middle lower Pleistocene</td><td>Fine sand (Sabbie di Asti)</td><td>Transitional marine</td><td>PLMa</td><td>GEL (2,6My)</td><td>-</td><td>900</td><td>1000</td><td>1300</td></tr><tr><td>Neogene</td><td>Middle upper Pliocene</td><td>Clays (Argille del Santerno)</td><td>Marine</td><td>PL</td><td>PL (5,3 My)</td><td>-</td><td>1200</td><td>1320</td><td>1800</td></tr><tr><td>Neogene</td><td>Miocene</td><td>Chaotic layer (Sergnano)</td><td>-</td><td>MESb</td><td>ME3</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Neogene</td><td>Miocene</td><td>Marls (Marne di Gallare)</td><td>-</td><td>MESa</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

Dark grey cells indicate the stratigraphic units proposed by the Carg project (ISPRA et al., 2016) locally developed in the framework of the 1:10.000 sheet 118 geological map. In correspondence of the studied area they are described down to the base of group aquifer D (Regione Lombardia, Eni Divisione Agip (2002); for deeper deposits, the stratigraphic units proposed by Geomol project (2015), are considered (dark cells). The main chronostratigraphic unconformities refer to the GeoMol project (2015). The hydrogeological units are described following Regione Lombardia, Eni Divisione Agip (2002). Depths of the chrono-stratigraphic unconformities (GeoMol Team, 2015) reported in the last three columns are referred to the countryside level.

surface units are characterized by coarse and medium gravels, with subordinate levels of medium and coarse sands, and belong to the Supersintema Lombardo Superiore (LS).

The stratigraphy units from the top to the bottom are (Figure 2, panels B,C):

(1) **LS** unit (Supersintema Lombardo Superiore, Upper Pleistocene-Holocene; ISPRA et al., 2016, **Figure 2**, panel B) or group or aquifer-A (Regione Lombardia, Eni Divisione Agip, 2002) is made of clast-supported gravel, medium-coarse gravel and sparse pebbles (Scardia et al., 2012). LS is marked at the base by the Qc3 (Figure 2, panel B) regional chronostratigraphic Unconformity (0,45 My, GeoMol Team, 2015) or Y-surface (Scardia et al., 2012) that defines a regional change in fluvial depositional style (i.e., proximal braidplain according to Scardia et al., 2012), suggesting the occurrence of important erosional processes. All the outcropping units of the 118-Milano geological sheet (Figure 2, panel A) are sub-units of LS. LS unit has an average thickness of 50 m (Figure 2, panel b and Table 2);

(2) **LI** unit (Supersintema Lombardo Inferiore, Upper Pleistocene; ISPRA et al., 2016, **Figure 2**, panel B) or group acquifer-B (Regione Lombardia, Eni Divisione Agip, 2002) is an alternation of medium to coarse grained sand, gravel and subordinate silt, interpreted as a distal braidplain with wandering fluvial channels (Scardia et al., 2012). At the bottom, LI is marked by Qc1 (Figure 2, panel B), a regional chronostratigraphic Unconformity (0,87 My, GeoMol Team, 2015) or R-surface (Muttoni et al., 2003), correlated to the onset of major Pleistocene glaciations in the Alps. The

average thickness of LI unit ranges from about 50 to 75 m moving NNW-SSE (Figure 2, panel B and Table 2);

(3) **PD** unit (Padano Supersynthem, Middle lower Pleistocene; ISPRA et al., 2016, **Figure 2**, panel C) is made of marine coastal and platform deposits (sand, with some clay layers) at the base, that evolve upwards to delta, delta-conoid and lagoon sediments (fine and very fine sand, mud, clay with some organic material, locally gravel) and to completely continental deposits (meandering river depositional system in alluvial plain, Scardia et al., 2012).

In between, the transition from marine to continental environment is defined by the regional 1,25 My unconformity called Qm2 (GeoMol Team, 2015, **Figure 2**, panel C), corresponding to the bottom of group acquifer-C (Regione Lombardia, Eni Divisione Agip, 2002); at the base, PD is marked by the 1,5 My unconformity called Qm1 (GeoMol Team, 2015) that corresponds to the base of group acquifer-D (Regione Lombardia, Eni Divisione Agip, 2002). The average thickness of the PD unit ranges from about 500 to 600 m moving NNW-SSE (Figure 2, panel c and table 2);

(4) **PLMa** unit (Middle lower Pleistocene, Geomol Team, 2015, **Figure 2**, panel C) is made of sand deposits with some clay levels and thin gravel layers saturated with salty water (i.e., "Sabbie di Asti" formation, Regione Lombardia, ENI Divisione Agip 2002). At the bottom, PLMa unit is marked by the 2,6 My GELasian Unconformity (GeoMol Team, 2015). The average thickness of the PLMa unit range from about 300 to 600 m moving NNW-SSE (Figure 2, panel C and Table 2);
(5) PL unit (early-middle Pliocene, Geomol Team, 2015, **Figure 2**, panel C) is made of mud and clay deposits (i.e., “Argille del Santerno” formation, Dondi et al., 1982). At the base, the PL unit is marked by the 5,3 My intra Zanclean Unconformity called PL (“base of Pliocene”, GeoMol Team, 2015). The average thickness of the PL unit ranges from about 300 to 500 m moving NNW-SSE (**Figure 2**, panel C and **Table 2**);

(6) MESa unit (Burdigalian-Tortonian, GeoMol Team, 2015, **Figure 2**, panel C) is characterized by the “Marne di Gallare” formation (Rizzini and Dondi, 1978), that represents the geological bedrock of the area. In between, the formation is broken off at the top by the ME3 intra Messinian Unconformity (Geomol Team, 2015, **Figure 2**, panel C), that marks the transition to the Messinian “Ghiaie del Sergnano” formation (chaotic gravel and sand, cemented levels, marls and sandstone, ViDEPI project 2009), that represents the potential reservoir between the Miocenic marls and the caprock of “Argille del Santerno,” a target of oil exploration.

The subsoil schemes describing the geological cross section of **Figure 2** (panels B, C) are presented in **Table 2**. In particular, the depths of the Qc3, Qc1 and PL unconformities (Geomol project, 2015) are well-constrained by the isobaths published in the framework of the Carg project (https://www.isprambiente.gov.it/Media/carg/118_MILANO_SOTTO/Foglio.html, ISPRA et al., 2016), while the depths of the QM2, QM1, GEL unconformities (Geomol project, 2015) are deduced on the basis of the stratigraphy of ViDEPI project (2009).

## Noise Horizontal-To-Vertical Spectral Ratio

From March to June 2021, 33 single-station ambient noise measurements are performed in the MMA (**Figure 2**, panel A). The sites selected for the measurements cover an area of about 15 km² and are characterized by minimum and maximum inter distances between each measure of ~1 and 4 km, respectively. Although it is not possible to plan a regular grid, the measurements are devised to avoid significant sources of anthropic noise such as industrial facilities, railways etc. Additional data are available from the continuous data stream recorded by the IV. MILN permanent station (**Figure 2**, panel A), equipped with a broadband Nanometrics Trillium 40s seismometer (https://www.nanometrics.ca/) and a Kinemetrics Episensor-FBA accelerometer (https://kinemetrics.com/).

The measurements are performed using a 6-channel Reftek-130 24-bits datalogger (https://reftek.com/) coupled with two Lennartz 5s seismometers (https://www.lennartz.de/en/) with sampling frequency of 100 Hz. Each measurement has a duration spanning from at least 30 min to about 1 h. All measures are performed during working days and in absence of wind. In any case, the couple of seismometers are always covered by a cap and, if possible, partially buried, in order to avoid any disturbances that may adversely affect the results. (e.g., Mucciarelli et al., 2005). The measures are performed in all cases by expert operators and each single point is a priori checked in order to avoid the proximity of underground

services or important skyscrapers. In this regard, a comprehensive description of the good practices for the acquisition and the analysis of ambient noise vibrations are provided by Molnar et al., 2022.

Considering, moreover, the SESAME (SESAME, 2004; Bard, 2008) recommendations for reliable HVSR curves, the first two criteria, based on the availability of adequate data windows length and number of cycles, are always satisfied. Concerning the third criterion, based on the standard deviation of the HVSR amplitude in dependence of the frequency, few measuring points for *pk1* (see black cells in **Supplementary Table S1**) exceed the limit of 3 at frequency <0.5 Hz; at higher frequencies, solely for the measurement 19, the third criterion is not satisfied in case of *pk2* (see black cell in **Supplementary Table S1**).

The data processing removes the mean and the linear trend partitioning the signal in time sub-windows with increasing duration, from 60 to 240 s. A cosine taper function (5%) is applied to better highlight possible low frequency peaks (<0.5 Hz).

As shown in **Figure 4** (panel A), in order to avoid low resolution or aliasing phenomena at low frequency, a minimum window of analyses of 60 s is necessary. After removing possible transient signals identified through visual inspection, the fast Fourier transform (FFT) is calculated for each time window, which is smoothed using the Konno and Ohmachi (1998) algorithm with the b coefficient set to 40 (**Figure 4**, panel B), allowing to distinguish better the single peaks observed in the HVSR curves. Independent of the processing, a further important peculiarity for the evidence of low-frequency peaks in the Po Plain is the atmospheric pressure condition where the measurements are performed. The ambient noise level in the area undergo significant daily and seasonal variations with evident consequences on HVSR results at frequencies <0.5 Hz (e.g., Marzorati and Bindi, 2006). An example is reported in **Figure 4** (panel C) by comparing the FFTs of ambient noise windows recorded in October 2021 at IV. MILN station, during sunny and stormy days, respectively.

To calculate the HVSRs, for each single sub-window the spectral ratio between the root-mean-square average spectra of the horizontal components over the vertical components is calculated (Nakamura, 1989). Finally, the mean and standard deviation are computed.

For completeness, the results obtained for the 33 ambient noise measurements are illustrated in the **Supplementary Material**, in figures from **Supplementary Figure S1–S6**, while **Supplementary Table S1** indicates all frequency peaks (i.e., *pk*) and related uncertainties recognized at urban scale by multiple measurements.

**Figure 5** (panel A) shows the results in terms of mean ± standard deviations (black solid and dotted lines, respectively) of the 33 ambient noise measurements, together with the mean HVSR curves obtained from independent sets of measurements performed in the northern (i.e., North of M002, **Figure 2**, panel A), central and southern (i.e., South of M026, **Figure 2**, panel A) parts of the urban area, respectively, and the mean HVSR curves representative of the seismic passive arrays performed at North Park, Giuriati sports center and Vettabbia park (**Figure 2**, panel A).
<figure><img src="image_70.png" /></figure>

**FIGURE 4** | Panel (A) example of ambient noise HVSRs (Nakamura, 1989) calculated at IV. MILN station during a non-working day (i.e., 8th December 2021) using time sub-windows of increasing length of 20, 60 and 240 s, respectively; Panel (B) example of ambient noise HVSRs (Nakamura, 1989) calculated at IV. MILN station considering FFTs obtained by the Konno-Ohmachi smoothing windows, with *b* equal to 20 and 40; Panel (C) example of Power Spectral Density (PSD) calculated at IV. MILN station from ambient noise recorded at the end of October 2021 during a sunny and a stormy day, respectively.

In general, the *HVSR* curves (**Figure 5**, panel A) show a high degree of coherence, which showed 3 main amplified bands of frequency that are from 0.17 to 0.23 Hz (*pk1*, **Supplementary Table S1**), from 0.48 to 0.66 Hz (*pk2*, **Supplementary Table S1**) and from 3.2 to 7.6 Hz (*pk3*, **Supplementary Table S1**). In general, the *HVSR* curves show a slight decrease in frequency of the three main peaks moving from the Northern to the Southern part of the area.

Considering, the SESAME (2004) criteria for clear HVSR peaks, *pk1* is the resonance frequency that shows a more reliable peak, with 2–4 satisfied criteria for about 70% of the measuring points (**Supplementary Table S1**). *Pk3*, instead, is the resonance peak that shows the highest number of satisfied criteria for a clear HVSR peak (5 satisfied criteria, **Supplementary Table S1**), even if for a lower percentage of measuring points (from 3 to 5 of satisfied criteria for about 40% of measuring points, **Supplementary Table S1**). *Pk2* is the resonance frequency having more uncertainties concerning clear HVSR peaks. In particular, concerning the HVSR amplitude, except *pk1*, which shows values > 2 for the 80% of the measuring points (i.e., the lower bound for significant peaks), at higher frequency more variability is observable. All frequency peaks showing amplitude >2 are indicated in bold in **Supplementary Table S1**.

In general, the *HVSRs* variability in terms of both resonant peaks and amplitude at high frequency (>1 Hz) could be

associated with very local factors (i.e., local soil conditions, anthropic ambient noise levels etc.). For example, **Figure 5** (panel B) shows a comparison between the *HVSR* calculated for the new 33 ambient noise measurements and the *HVSR* calculated from ambient noise recorded by the Trillium 40s installed at IV. MILN station, both during a non-working day (i.e. 8th December 2021) and during the lockdown against the Covid-19 pandemic (i.e. April 2020, Sunday 6th and April 2020, Monday 7th). This example shows the strong influence in an urban area of the anthropic noise on the vertical component of motion and, consequently, on the final results in terms of *HVSR*. The relative Power Spectral Density (PSD) is shown in the **Supplementary Material (Supplementary Figure S7)**.

At low frequency (<1 Hz), the noise level could be affected by variations due to differences in the environmental conditions during the period when the measurement is performed (e.g., Marzorati and Bindi, 2006).

Finally, to detect possible polarization effects of the horizontal components, Fourier spectra are computed by rotating the north–south (N-S) and east-west (E-W) noise components clockwise from the North, from 0° to 180° with step of 10°. For each direction, the *HVSR* is computed considering the ratio between the Fourier spectrum of each rotated horizontal component and the Fourier spectrum of the vertical component. **Figure 5** (panel C) shows the rotated *HVSR* obtained at North
<figure><img src="image_71.png" /></figure>

**FIGURE 5** | Panel (A) HVSRs relative to the 33 new ambient noise measurements (light grey) performed in the Milan urban area from March to June 2021 (in black mean amplification and ±1 standard deviation). Blue, red and green solid lines indicate the mean amplification functions for the north, central and south parts of the target area. Blue, red and green dotted lines indicate the amplification functions for the three target sites of North Park (north of the area), Giuriati (center part) and Vettabbia park (south of the area); Panel (B) comparison in terms of HVSRs between the new ambient noise measurements (light grey; in black mean amplification and ±1 standard deviation) and the results obtained at IV. MILN station using time windows selected both recently (i.e., December 2021) and during the Covid-19 lockdown (i.e., April 2020). Panel (C) directional ambient noise HVSRs obtained for North Park (left), Giuriati (center) and Vettabbia park (right).

park, Giuriati sports center and Vettabbia park, where in all cases just negligible polarization in directions ranging from 40°N to 120°N is detected around 0.2 Hz.

**Figure 6** shows, in the map, the variability of the more relevant frequency peaks detected at high (pk3, panel A) and low (pk1, panel B) frequency, highlighting decreasing values moving from the northern to the southern part of the studied area.

## Earthquake Horizontal-To-Vertical Spectral Ratio

The results obtained from ambient noise HVSRs are compared with those from earthquake HVSRs (e.g., Lermo and Chavez-Garcia, 1993), calculated considering the set of seismic events reported in Table 1. The selected earthquakes, with magnitude (Mw) ranging from 3.7 to 5.8, are recorded by the IV. MILN station and occurred in the Po Plain area and surroundings from the 1st January 2012 until today. The velocimetric data selected for the analyses are downloaded from the EIDA web site (http://www.orfeus-eu.org/data/eida/).

The simple data processing involved the application of a baseline operator to the entire record to obtain the zero-mean signal and subsequently to remove any linear trends, a tapering cosine function (0.05%) applied at the beginning and the end of each selected window and a 4th order Butterworth band-pass

acausal (Boore and Akkar, 2003) filter to remove the high and low-frequency noises; for each analysed event the filter cut-off thresholds are selected by a visual inspection.

Given the long-period effects due to the low frequencies propagation in the Po plain basin during the 2012 Emilia seismic sequence (e.g., Luzi et al., 2013), the earthquakes HVSRs are computed on both the S-waves train and the coda. Time windows from 10 to 20 s are selected for the S-waves starting from their onset after a visual inspection. Time lengths for the surface train waves are determined through a time-frequency analysis on each single earthquake. Finally, coda windows ranging from 15 to 60 s are selected at the end of the S-phase depending on each earthquake's magnitude, distance and ray-path.

An example of propagation effects on train waves crossing along the Po plain alluvial deposits is shown in Figure 7 (bottom panels), where the EW velocimetric waveforms recorded at IV. MILN station during the 20th May 2012, Mw 5.8, Emilia mainshock (Repi 169 km) and during the 18th December 2021, Dalmine, Mw 3.9, event (Repi 35 km, see Table 1) are shown; in the first case, the example shows the generation of surface waves with the consequence of long earthquake ground motion duration and PGV recorded on coda. In the second case, the hypocentre is located at a depth of 31 km and the source to site distance is 35 km, so no surface waves are generated and the PGV is found on the S-phase.
<figure><img src="image_72.png" /></figure>

**FIGURE 6** | Resonance frequency intervals for *pk3* (panel A) and *pk1* (panel B) *HVSR* ambient noise peaks (see **Supplementary Table S1**) superimposed on the geological map of the area (Figure 2, panel A). Measuring points that not satisfied the SESAME (2004) criteria for *clear HVSR peaks* (see black and dark grey cells in **Supplementary Table S1**) are not considered. Top and bottom panels also indicate the isobaths of *Qc3* (base of the upper Padano Supersynthem, see **Figure 2B** and **Table 2**) and *PL* (base of Pliocene, **Figure 2C** and **Table 2**) unconformities as recognized by the Geomol project (2015).

In detail, the *HVSRs* are calculated by splitting the earthquakes in two different sets of data, the first including the 2012 Emilia earthquakes with Mw > 5, epicentral distances ranging from 148 to 191 km (**Table 1**) and homogeneous azimuth of ~290°N (**Figure 1**, panel B); the second including earthquakes occurring mainly at the edge of the basin, with Mw < 4.5, maximum epicentral distance of 142 km (**Table 1**) and sparse azimuths in the range 2°–352°N (**Figure 1**, panel B).

The earthquakes HVSRs are shown in **Figure 7**, where panels A, B show the curves for the higher and lower magnitude set of events, respectively. In particular, top and central panels show the results in terms of mean amplitude calculated from S-phase and coda, respectively. The results obtained considering the ambient noise and the set of earthquakes with magnitude >5 show that the two *HVSR* curves agree at frequencies <0.5 Hz, particularly

considering the coda waves for the higher magnitude set of events.

In case of the available local earthquakes, characterized by short distances and low magnitude, the windows lengths do not always assure the calculation of stable HVSR curves for frequencies lower than ~0.5 Hz (**Figure 7**, right panels). This issue is even more relevant in case of high levels of ambient noise, such as at the IV. MILN station during working days (i.e., -110 < dB < -90 in the range 1–20 Hz).

The most relevant result comparing of earthquake and ambient noise *HVSRs* is that, at frequencies between 1 and 2 Hz, the ambient noise *HVSRs* do not show any amplification peak. In contrast, S-phase and coda depict clear amplification peaks for both set of events. Concerning the amplitude values, earthquakes *HVSRs* calculated both from S-phase and coda show in general higher amplitudes with respect to ambient noise, in particular for frequencies >0.5 Hz.

Concerning the possible influence of the different azimuths of the local earthquakes reported in **Table 1**, just in case of the two events of Mw = 3.9, recorded at IV. MILN station on 2021/12/18 10:34:47 (epicentral distance 35 km, depth 26 km and azimuth 293°N) and on 2020/12/17 15:59:22 (epicentral distance 8 km, depth 56 km, azimuth 83°N) higher amplitudes along the related ray-paths are detected on S-phase in the frequency range 1–8 Hz.

As in the case of ambient noise *HVSRs*, no polarization phenomena due to the recording site are observable from earthquake signals.

## Geophysical Surveys

To characterize the shear wave velocity with depth, active and passive geophysical tests are carried out.

Two Multichannel Analyses of Surface Waves (MASW, Park et al., 1999) are performed at North public park and at Vettabbia public park using in both cases 24 vertical geophones with eigen natural frequency of 4.5 Hz and a vertical impact source (6.8 kg hammer). In both sites, the geophones are set out in a linear array varying the inter geophonic distance from 2.5 to 5 m (total length of 60 and 120 m respectively) and using a minimum offset moving from 5 to 10 m. The linear arrays allow investigating depths ranging from 30 to 35 m. The Rayleigh wave dispersion of the acquired data set is obtained from the stacking of 7 shots. For each shot, a sampling frequency of 2000 Hz and a time of acquisition of 2 s are chosen. To generate the phase velocity spectrum of Rayleigh waves, the phase shift algorithm (Park et al., 1998) is used.

In general, as expected, in the configuration with the greater geophonic inter distance (5 m) an improvement in signal correlation is observed in the low frequency phase velocity spectrum. In particular, at Vettabbia Park, in the configuration with the smaller geophonic inter distance (2.5 m) an improvement in signal correlation is also observed in the high frequency phase velocity spectrum (20–30 Hz). The cause of this observation is probably due to the attenuation properties of the more superficial soil layers that resulted, from local geognostic investigations (https://www.geoportale.regione.lombardia.it/), to be mostly incoherent fine deposits. From the phase velocity spectra obtained with the two configurations (5–2.5 m), the
<figure><img src="image_73.png" /></figure>

**FIGURE 7** | Earthquake HVSRs calculated on the events listed in **Table 1**. Panel (A) shows the results obtained on a sub set of events with magnitude >5 and epicentral distances from 127 to 191 km, while panel (B) indicates the results obtained considering events with magnitude <4.5 and epicentral distances from 8 to 143 km; the vertical orange dotted lines indicate the lower limit for stable HVSR curves obtained by using the available local earthquakes. The epicentres of the considered events are reported in the panel (B) of **Figure 1**. Bottom panels show two examples of near and far field waveforms recorded at IV. MILN station for the 2012/12/18 Dalmine M<sub>w</sub> 3.9 event (right) and for the 2012/05/20 Emilia M<sub>w</sub> 5.8 mainshock (left).

fundamental mode picking covered the frequency range 6–30 Hz at Vettabbia park and the frequency range 5.5–30 Hz at North park (Supplementary Figures S8 and S9 respectively). The Vs.30 calculated from MASW at North and Vettabbia park results equal to 345 m/s and to 265 m/s respectively, allowing classifying the analysed sites in C class (i.e. Vs.30 < 360 m/s) of Italian (NTC, 2018) and European (CEN, 2004) seismic codes.

To define S-wave velocity profiles for larger depths, 4 different microtremor arrays, with increasing apertures and different layouts are installed at North public park (Figure 8, panel A), Giuriati sports center, very close to the IV. MILN station (Figure 8, panel B), and Vettabbia public park (Figure 8, panel C). At Giuriati sports center, two microtremor arrays are carried out using a circular geometry with radii of 12 and 50 m, respectively (Figure 8, panel B). The layout of the two arrays consisted of a central receiver surrounded by 8 stations deployed with an angular distance of 45°. The instrumentation comprised Lennartz-5s high-sensitivity (400 V/m/s) velocimeters coupled with a 24-bit Reftek 130 three-channel digitizer. The sampling frequency is set at 100 Hz, with the recording period for each array of 3 h.

For both North and Vettabbia public parks, one microtremor array is set up using the layouts reported in **Figure 8** (panel A, C) and 14 Lennartz-5s high-sensitivity (400 V/m/s) velocimeters coupled with a 24-bit Reftek 130 three-channel digitizer. The

sampling frequency and recording period are the same as that used at Giuriati. In this case, the two types of geometry enable work with inter distances among sensors in the range 5–600 m.

For each array, the first step of the analysis consists of a visual inspection of the recordings at each station. In particular, in order to identify malfunctioning and to select suitable signal windows for the surface wave analysis, the quality of the recording is evaluated by analysing the signal stationarity in the time domain, the relevant unfiltered Fourier spectra and the HVSR variation over time. During surface wave investigation, it is common practice to verify the reliability of the one-dimensional site structure assumption (Aki, 1957; Okada, 2003). For this reason, we estimated the HVSR at each station of the array and the stability of HVSR throughout the array stations is in the range of frequency (0.2–30 Hz). The analyses use synchronized signal windows, with lengths ranging from 60 to 120 s, extracted from the available recordings, avoiding time windows affected by local disturbances.

The Rayleigh wave dispersion curve relative to the fundamental mode is calculated for each seismic array, taking account of the vertical ground-motion components, using the modified spatial autocorrelation (mSPAC) method (Bettig et al., 2001), the f-k (frequency-wavenumber) beamforming method (Capon 1969; Lacoss et al., 1969) and the ESAC (Extended Spatial Auto Correlation) method (Ohori et al., 2002; Okada 2003).
<figure><img src="image_74.png" /><figcaption>**FIGURE 8** | Panel (A) location of the MASW active surveys and of the 2D passive array performed at North Park; Panel (B) location of the two 2D passive microtremor arrays performed at Giuriati sports center; Panel (C) the same as in panel (A), but for Vettabbia public park; Panel (D) Rayleigh-wave dispersion curves from 2D passive microtremor array and MASW surveys.</figcaption></figure>

Although the *f*-*k* methods are considered stable and robust algorithms, the *m*SPAC and ESAC methods give reliable estimates of the dispersion characteristics within a larger frequency band. In any case, the use of various combinations is always advisable (Bard 2005; Parolai et al., 2006; Puglia et al., 2011; SESAME, 2004).

It is worth mentioning that in this case no higher modes of vibration are detected. It is well known that by recognizing just the fundamental mode of ellipticity, it is generally difficult to correct for the energy of SH and Love waves present in the horizontal components of the ambient vibration wavefield (e.g., Poggi and Fäh, 2010). In this way, the quality of results also depends on the number of the available receivers, on the quality of the sensors installation, on the selected array geometry and, in particular, from the available levels of background ambient noise (Supplementary Figure S7). Moreover, in our case, also the application of the high-resolution frequency-wavenumber analysis and the use of the three-component of motion (e.g., Lai et al., 2020) does not allow to improve the results obtained by using the vertical component.

At Giuriati sports center, the experimental dispersion curve is obtained by using the *f*-*k* and *m*SPAC methods; in particular, the microtremor arrays with radii equal to 12 and 50 m, allow the definition of the dispersion curve in the intervals ~2.5–18 Hz and ~0.4–3.5 Hz, respectively (Figure 8, panel D).

At North and Vettabbia public parks, the experimental dispersion curves are obtained by using the *f*-*k*, *m*SPAC and ESAC methods; in these cases, the analyses on microtremor arrays allow defining the dispersion curves in the intervals ~0.3–3.0 Hz

and ~0.3–3.5 Hz, respectively; at higher frequencies, the dispersion curves are obtained by the interpretation of the MASW surveys in the intervals ~5.5–30 Hz and ~6.0–30 Hz, respectively (Figure 8, panel D).

The inversion procedure to obtain the 1D velocity profile is performed using the nearest neighborhood algorithm, as implemented in the Diver software (Wathelet et al., 2008). To obtain the best results from the available data, at each target site a joint inversion of the experimental dispersion and the ellipticity curves is carried out, since the single station HVSR contains information on the mean S-wave velocity and the thickness of the sedimentary cover (Parolai et al., 2005). The ellipticity curves are evaluated by using the single station analysis (HVSR), selecting of the main peaks or the associated flanks, dividing the amplitude by √2 to reduce the Love-wave and body wave contributions (Fäh et al., 2001; Fäh et al., 2003). At Vettabbia park and Giuriati sports center, the inversion procedure to obtain the Vs. profiles is constrained by using the resonance peaks around 0.17 and 0.19 Hz respectively, while at North park a portion of HVSR curve in the range 0.2–0.8 Hz is used.

## Shear Waves Velocity (V<sub>s</sub>) Profile

The Vs. profiles obtained for the three sites are shown in Figures 9, 10, 11 (panels A). Supplementary Figure S8 shows a detail of the Vs. profiles in the first 400 m of depth. The profiles are calculated using the dispersion curves obtained both from the passive arrays (Figures 9–11, panels B) and the active
<figure><img src="image_75.png" /></figure>

**FIGURE 9** | Panel (A) Vs. shear wave velocity profile obtained from the joint inversion (Rayleigh-wave dispersion curve and ellipticity curve) at North park. The main chrono stratigraphic units and the related unconformities (GeoMol Team, 2015), as reported in **Figure 2** (Panels B, C) and **Table 2**, are also indicated. The central panel shows a detail of the Vs. shear wave velocity profile in the first 400 m of depth. Panel (B) experimental Rayleigh-wave dispersion (black dotted line) relative to the fundamental mode of vibration compared to the theoretical models (grey scale lines). The white line indicates the best-fit. Panel (C) experimental ellipticity curve (black line) obtained from a portion of ambient noise HVSR in the frequency range 0.2–0.8 Hz. Grey scale lines indicate the theoretical models. White line indicates the best fit.

<figure><img src="image_76.png" /></figure>

**FIGURE 10** | The same as for **Figure 9**, but considering the same analyses performed at Giuriati sports center.

multichannel analysis of the surface waves, while the ellipticity contribution is estimated from a selected portion of the ambient noise HVSR curves (Figures 9–11, panels C).

The Vs. profile shows, in general, a similar trend of Vs. with depth and an increasing depth of the main contrasts of impedance towards South.

In particular, all Vs. profiles show a superficial layer with Vs. lower than 250 m/s and a first significant increase of Vs. to ~420 m/s (Vettabbia park), ~440 m/s (Giuriati sports center) and ~500 m/s (North park) at depths ranging from ~10 m (North park) to ~20 m (Vettabbia park).

A second slight increase of Vs. is observable at depth of ~35 m (North park), ~45 m (Giuriati sports center) and ~80 m (Vettabbia park), respectively, where the Vs. approaches the value of ~600 m/s at all investigated sites.

The Vs. values keep constant before a third interface located at depths ranging from ~300 m (North Park and Giuriati sports center) to ~400 m (Vettabbia park); at all sites, the Vs. values reach or exceed 800 m/s, that correspond to the engineering bedrock (i.e., site class A, Vs.30 > 800 m/s) according to the European (CEN, 2004) and the Italian (NTC, 2018) seismic codes.
Vettabbia park

<figure><img src="image_77.png" /><figcaption>FIGURE 11 | The same as for Figure 9, but considering the same analyses performed at Vettabbia public park.</figcaption></figure>

B

<figure><img src="image_78.png" /></figure>

C

<figure><img src="image_79.png" /><figcaption>TABLE 3 | Shear-wave velocities ($V_S$) and layer thickness ($H$) used for 1D modelling at North Park, Giuriati sports center and Vettabbia park, together with the relative theoretical fundamental modes ($f_{0(ID)}$). The $f_{0(ID)}$ values evidenced in bold indicate seismic impedance-contrasts between subsequent layers greater than 1.5. Values for unit weight ($\rho$), mean effective stress ($\sigma'_0$) and damping at low-strains ($D_o$), at each layer (ID) are also reported.</figcaption></figure>

<table><thead><tr><th rowspan="2"></th><th colspan="4">Vettabbia Park</th><th colspan="3">North Park</th><th colspan="3">Giuriati Sport Center</th><th colspan="3">All Sites</th></tr><tr><th>ID</th><th>VS. [m/s]</th><th>H [M]</th><th>f0>ID</th><th>VS. [m/s]</th><th>H [M]</th><th>f0>ID</th><th>VS. [m/s]</th><th>H [M]</th><th>f0>ID</th><th>ρ [t/m3]</th><th>σ'0 [Atm]</th><th>D0 [%]</th></tr></thead><tbody><tr><td rowspan="7">Best-fit</td><td>1</td><td>199</td><td>11</td><td>4.54</td><td>287</td><td>8</td><td>8.96</td><td>249</td><td>9</td><td>6.91</td><td>1.45</td><td>0.2</td><td>1.3</td></tr><tr><td>2</td><td>501</td><td>71</td><td>1.27</td><td>422</td><td>24</td><td>2.95</td><td>455</td><td>42</td><td>1.95</td><td>1.50</td><td>0.8</td><td>0.9</td></tr><tr><td>3</td><td>620</td><td>323</td><td>0.35</td><td>612</td><td>41</td><td>1.66</td><td>576</td><td>59</td><td>1.08</td><td>1.60</td><td>1.4</td><td>0.7</td></tr><tr><td>4</td><td>1080</td><td>346</td><td>0.24</td><td>620</td><td>272</td><td>0.42</td><td>580</td><td>239</td><td>0.39</td><td>1.65</td><td>2.7</td><td>0.6</td></tr><tr><td>5</td><td>1169</td><td>481</td><td>0.17</td><td>1126</td><td>489</td><td>0.25</td><td>773</td><td>572</td><td>0.18</td><td>1.70</td><td>11</td><td>0.4</td></tr><tr><td>6</td><td>1240</td><td>475</td><td>0.14</td><td>1292</td><td>342</td><td>0.20</td><td>1276</td><td>388</td><td>0.15</td><td>2.10</td><td>41</td><td>0.3</td></tr><tr><td>7</td><td>1604</td><td>-</td><td>-</td><td>2051</td><td>-</td><td>-</td><td>1645</td><td>-</td><td>-</td><td>2.50</td><td>56</td><td>0.3</td></tr><tr><td rowspan="4">min</td><td>1–3</td><td colspan="3">Same as above</td><td colspan="3">Same as above</td><td colspan="3">Same as above</td><td colspan="3"></td></tr><tr><td>4</td><td>1000</td><td>425</td><td>0.22</td><td>950</td><td>575</td><td>0.22</td><td>700</td><td>651</td><td>0.16</td><td colspan="3"></td></tr><tr><td>5</td><td>1070</td><td>550</td><td>0.15</td><td>1060</td><td>430</td><td>0.17</td><td>1050</td><td>400</td><td>0.14</td><td colspan="3"></td></tr><tr><td>6</td><td>1170</td><td>420</td><td>0.12</td><td>1700</td><td>-</td><td>-</td><td>1500</td><td>-</td><td>-</td><td colspan="3"></td></tr><tr><td rowspan="4">max</td><td>1–3</td><td colspan="3">Same as above</td><td colspan="3">Same as above</td><td colspan="3">Same as above</td><td colspan="3"></td></tr><tr><td>4</td><td>1160</td><td>285</td><td>0.26</td><td>1300</td><td>425</td><td>0.27</td><td>900</td><td>451</td><td>0.23</td><td colspan="3"></td></tr><tr><td>5</td><td>1250</td><td>465</td><td>0.19</td><td>1450</td><td>280</td><td>0.23</td><td>1500</td><td>400</td><td>0.18</td><td colspan="3"></td></tr><tr><td>6</td><td>1350</td><td>465</td><td>0.15</td><td>2300</td><td>-</td><td>-</td><td>1850</td><td>-</td><td>-</td><td colspan="3"></td></tr><tr><td>7</td><td>1890</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

At larger depths, ranging from ~850 m (North park) to ~1250 m (Vettabbia park), the Vs. values reach ~1150 m/s to ~1250 m/s at all investigated sites. The deepest and observable Vs. interface is located at depths varying from ~1200 m (North park) to ~1700 m (Vettabbia park), although associated with large uncertainties. At these depths, the Vs. profile of the Giuriati sports center (Figure 10, panel A) follows an intermediate trend with respect to the results obtained at North (Figure 9, panel A) and Vettabbia (Figure 11, panel A) parks.

For all investigated sites, both the Vs. values and the depths of the main interfaces appear better constrained down to ~400 m. In

particular, in the uppermost 30 m the Vs. values are constrained by the dispersion curves obtained by both active and passive measurements, while the available stratigraphy from the shallow (https://www.geoportale.regione.lombardia.it/) and deep wells (Carg project, ISPRA et al., 2016; ViDEPI project, 2009), and the available literature (e.g., De Franco et al., 2009; Scardia et al., 2012) allow us to constrain the geological model. In general, the uncertainties of the Vs. models increase with increasing depths, in particular in the southern part of the city (i.e., Vettabbia park), where a limited amount of information is available. The Vs. values and the corresponding thickness of the model showing the
<figure><img src="image_80.png" /><figcaption>**FIGURE 12** | Theoretical transfer functions (TFs) obtained through 1D numerical modelling based on the linear field assumption at North Park (A). Vettabbia park (B) and Giuriati (C,D). Red, yellow and violet lines indicate the TFs obtained considering the best Vs-model and two Vs. profiles (i.e. min, max) that envelop the Vs. values with uncertainties of 10% with respect to the best fit. The vertical grey bars indicate for each site and layer (ID) the values of the theoretical resonance frequencies associated to the fundamental-modes $f_{0,\text{ID}}$, obtained by eq. 1 (see also **Table 3**). All panels show the comparison with the amplitude functions obtained from ambient noise HVSR (blue lines). At Giuriati also earthquakes HVSRs are reported (D).</figcaption></figure>

best fit (white lines in **Figures 9–11**, panels A) are reported in **Table 3**.

## 1D Numerical Modeling

Given the available ground-motion recordings (i.e., maximum PGA = 36 gal), a linear site response model (i.e., within the low-strain field) is preferred to validate the proposed shear wave velocity models (**Figures 9–11**, panels A), by comparing the theoretical transfer functions (TFs) and the fundamental modes with the resonance peaks highlighted by ambient noise and earthquake HVSRs. The 1D linear site response is done by means of the software SHAKE91 (Schnabel, 1972; Idriss and Sun, 1993; Schnabel et al., 1993). In **Table 3**, for each layer (ID) and site, are also reported the values of the theoretical resonance frequencies associated to the fundamental-modes $f_{0,\text{ID}}$, obtained by the equation:

$$f_{0,\text{ID}} = V_{\text{ST,\textit{ID}}} / (4 * H_{\text{T,\textit{ID}}}) \quad (1)$$

where $H_{T,\text{ID}}$ is the sum of the layer thicknesses ($H_{\text{ID}}$), starting from the topmost layer and moving towards, and $V_{\text{ST,\textit{ID}}}$ is the weighted average of shear-wave velocities over $H_{T,\text{ID}}$, obtained by the equation:

$$V_{\text{ST,\textit{ID}}} = H_{T,\textit{ID}} / \sum (H_{\text{ID}}/V_{\text{S,\textit{ID}}}) \quad (2)$$

where $H_{\text{ID}}$ and $V_{\text{S,\textit{ID}}}$ are the thickness and the shear-wave velocity, respectively.

**Table 3** also indicates the values for unit weight ($\rho$), useful for the impedance-contrast calculation, together with an estimation of the mean effective stress ($\sigma'_{0}$) at each layer (ID), to be used in the equation by Darendeli (2001) for the low-strain damping ($D_{0}$) calculation. In particular, the parameters of the Darendeli equations are fixed to: plasticity index: 0.0%; over-consolidation ratio: 1.0; cycles of loading: 10; characteristic frequency: 1 Hz. The mean effective stress ($\sigma'_{0}$) is estimated from unit weight ($\rho$) considering a coefficient of lateral earth pressure equal to 0.50.

It is well known that the choice of the seismic input for a low-deformation field is theoretically not relevant. Then, as the seismic input for modeling, a generic strong-motion waveform is selected, just characterized by enough energy to overcome the numerical limits of the software in the frequency band of interest (i.e., 0.1–10 Hz). In fact, the fundamental modes estimated by eq. (1) are within the same frequency band of about 0.1–10 Hz, for the numerical models in question.

Finally, in order to also account for the variations of both $V_{\text{S}}$ and $H$ in the modeling, **Table 3** reports the two further $V_{\text{S}}$ models min and max. These models represent the minimum and maximum limit of the envelope that includes all $V_{\text{S}}$ profiles (i.e., black lines in panels a of **Figures 9–11**) obtained with
uncertainties of 10% with respect to the best-fit model (i.e., white line in panels A of **Figures 9–11**).

The theoretical transfer functions are computed using the best model, together with the min and max models and are plotted in **Figure 12**, applying the Konno-Ohmachi smoothing (1998) operator with $b = 40$. The TFs are compared with the ambient noise *HVSRs*; except in the case of Giuriati sports centre, which is sited very close to IV. MILN (~300 m), where the *HVSR* from earthquakes is considered (**Figure 12D**).

For all the analysed sites, the numerical model matches the resonance frequencies recognized by the experimental results, in particular for *pk1* and *pk3*. On the contrary, more uncertainties are present in the range 0.5–3 Hz, where many peaks and valleys are present at different periods, as a consequence of the one-dimensional modeling adopted, which is composed by 7 homogeneous layers (**Table 3**). For all cases, a significant discrepancy in terms of peak amplitudes is evident by comparing the obtained TFs and the ambient noise *HVSR* curves.

In particular, both at North and Vettabbia parks the resonance frequencies obtained from the noise *HVSRs* at ~ 0.2 Hz better match the resonant frequency of the TFs derived from the min model (**Table 3**). At Giuriati site, the mean earthquake *HVSR* agrees with TFs in terms of frequency and amplitudes, with except of the peak at ~ 0.2 Hz, where the TFs slightly underestimate the experimental *HVSRs* obtained both from ambient noise and earthquakes.

## DISCUSSIONS AND CONCLUSION

In this paper, we present the work carried out to characterize the spatial variability of local soil conditions in the city of Milan and its surroundings and propose a seismo-stratigraphic model oriented to seismic site amplification. The model is derived from a combination of passive and active geophysical tests supported by a dense grid of geological and stratigraphic information from shallow and deep vertical wells and it is cross-checked with 1D numerical modeling.

The work consists of performing different geophysical tests; both active (i.e., MASW) and passive measurements are carried out to investigate both shallow and deep deposits. The passive analyses comprise 33 single-station ambient noise measurements and 4 microtremor arrays at 3 sites (i.e., North Park, Giuriati, Vettabbia) selected along a NNW-SSE alignment. The passive arrays are performed using different geometries and increasing aperture to explore the frequency range from 0.35 to 30 Hz. The geophysical analysis is supported by a subsoil model built along the NNW-SSE section after the collection and re-interpretation of all available geological and stratigraphic information and the documentation related to surficial and deep wells for exploration. Therefore, the resulting Vs. profiles, down to ~1.8 km in depth, can be associated with the geophysical interfaces and stratigraphic discontinuities.

The transfer functions obtained from numerical models and cross-checked by experimental curves (*HVSR* from ambient noise and earthquakes) identify a first resonant frequency at about 0.2 Hz, namely the expression of a deep contrast of impedance,

that we associate with the PL transition (**Figures 10–12**) between the Miocene sedimentary rocks (i.e., Gallare marls) and the overlying Pliocene deposits (i.e., Santerno clays).

The second contrast of impedance, relevant for site effects, is found at a depth between 300 and 400 m (**Figures 10–12**), at the transition between marine and continental sedimentation. In this case, a chronostratigraphic unconformity, Qm2 (**Figures 10–12**), is functional to geophysics, although the corresponding frequency is not evident from ambient noise *HVSRs*.

Shallower layers with low shear wave velocity (<250 m/s, **Figures 10–12**) result in high frequency peaks (between 3 and 8 Hz) that are evident in the ambient noise *HVSRs*.

A general trend, from all measurements, is the decrease in the frequency of the three main identified peaks, interpreted as the increase in depth of the deposits towards the centre of the Po plain, to the south.

Our hypotheses on the Vs. model should be verified due to the uncertainties and trade-off typical of the passive techniques (e.g., Lai et al., 2020). As an example, De Franco et al. (2009) interpret the subsoil of Milan by analysing a very deep active reflection profile, performed in the western part of the city and, in that case, except for superficial impedance contrast, the Vs. gradually increase down to ~900 m without marked changes.

Another source of uncertainty can be the amplification peak around 0.2 Hz, which is interpreted as the coincidence with the microseismic long-period double-frequency peak by Marzorati and Bindi (2006) or Ferretti et al. (2013). Other authors, e.g., Paolucci et al. (2015), highlight a correlation between the amplitude values of the low-frequency peak and the sea-wave height, not observed for the high-frequency peaks. This issue could be resolved performing a generalised inversion of seismic data (e.g., Parolai et al., 2000; Ameri et al., 2011).

Moreover, the 1D approach which, in our opinion, is acceptable for the target area because surface stratigraphy is typically characterized by horizontal layers of alluvial deposits, does not account for the generation of surface waves, and the response of soil deposits is assumed to be predominantly related to *SH* waves propagating vertically from the underlying bedrock. Numerical modelling performed in the Po Plain highlighted that 2D and 3D effects can be significant because of the complex buried morphology of the basin (e.g., Vuan et al., 2011; Paolucci et al., 2015; Klin et al., 2019; Mascandola et al., 2021). In particular, due to the very large dimensions of the Po basin few attempts are available in the literature to objectively quantify 2D or 3D morphological effects on wave field propagation, in particular at low frequency, which makes this still an open issue.

This study will be useful for future site response assessments, numerical modelling of seismic-wave propagation, dynamic ground response analyses, and site-specific seismic hazard evaluation, contributing to future urban planning and risk mitigation (Eurocode 8, CEN, 2004; NTC, 2018).

## DATA AVAILABILITY STATEMENT

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
## AUTHOR CONTRIBUTIONS

This study started from an original idea of LL, MM and SL. MM, SL, GB, AL and CF contributed to the data acquisition of the passive seismic arrays and MAWS. All the Authors participated to the ambient noise single measurements performed in the urban area of Milan. MM analysed ambient noise and earthquakes data through the HVSR method. MM, SL and RP worked on the ambient noise arrays to estimate the surface wave dispersion curves. MM and RP worked on the shear wave velocity profiles. RP performed the 1D numerical modeling. SL interpreted the geological and stratigraphic data and worked about the subsoil model.

## FUNDING

National Institute of Geophysics and Volcanology (INGV).

## REFERENCES

Aki, K. (1957). Space and Time Spectra of Stationary Stochastic Waves, with Special Reference to Microtremors. *Bull. Earthq. Res. Inst.* 35, 415-457.

Ameri, G., Oth, A., Pilz, M., Bindi, D., Parolai, S., Luzi, L., et al. (2011). Separation of Source and Site Effects by Generalized Inversion Technique Using the Aftershock Recordings of the 2009 L'Aquila Earthquake. *Bull. Earthq. Eng.* 9, 717-739. doi:10.1007/s10518-011-9248-4

Ansal, A., Tönük, G., and Kurtuluc, A. (2009). "Microzonation for Urban Planning," *Earthquakes and Tsunamis, Geotechnical, Geological, and Earthquake Engineering*. Editor A. T. Tankut (Dordrecht: Springer), Vol 11, 133-152. doi:10.1007/978-90-481-2399-5_9

Ansal, A., Kurtuluş, A., and Tönük, G. (2010). Seismic Microzonation and Earthquake Damage Scenarios for Urban Areas. *Soil Dyn. Earthq. Eng.* 30, 1319-1328. doi:10.1016/j.soildyn.2010.06.004

Aversa, S., and Crespellani, T. (2016). Seismic Microzonation: an Essential Tool for Urban Planning in Seismic Areas. *Upl. J. Urban Plann Landsc. Environ. Des.* 1 (1), 121-152. doi:10.6092/2531-9906/5035

Bard, P.-Y. (2008). Foreword. *Bull. Earthq. Eng.* 6 (1), 1-2. doi:10.1007/s10518-008-9059-4

Bard, P. Y. (2005). *Site Effects Assessment Using Ambient Excitations (SESAME). WP 6: Derivation of Dispersion Curves*. Grenoble: Technical report.

Beck, J. L., and Hall, J. F. (1986). Factors Contributing to the Catastrophe in Mexico City during the Earthquake of September 19, 1985. *Geophys. Res. Lett.* 13, 593-596. doi:10.1029/gl013i006p00593

Bettig, B., Bard, P. Y., Scherbaum, F., Riepl, J., Cotton, F., Cornou, C., et al. (2001). Analysis of Dense Array Measurements Using the Modified Spatial Auto-Correlation Method (SPAC). Application to Grenoble Area. *Bol. Geofis. Teor. Appl.* 42 (3-4), 281-304.

Bielak, J., Xu, J., and Ghattas, O. (1999). Earthquake Ground Motion and Structural Response in Alluvial Valleys. *J. Geotechnical Geoenvironmental Eng.* 125 (5), 413-423. doi:10.1061/(asce)1090-0241(1999)125:5(413)

Bigi, G., Costantino, D., Parotto, M., Sartori, R., and Scandone, P. (1990). *Structural Model of Italy. Firenze, Società Elaborazioni Cartografiche (S.EL.CA.), Consiglio Nazionale Delle Ricerche Progetto Finalizzato Geodinamica, Scala 1:500.000, 9 Fogli.*

Bocherdt, R. D. (1970). Effects of Local Geology on Ground Motion Near San Francisco Bay. *Bull. Seismol. Soc. Am.* 60, 29-61. doi:10.1785/BSSA0600010029

Boore, D. M., and Akkar, S. (2003). Effect of Causal and Acausal Filters on Elastic and Inelastic Response Spectra. *Earthq. Engng. Struct. Dyn.* 32, 1729-1748. doi:10.1002/eqe.299

Capon, J. (1969). High-resolution Frequency-Wavenumber Spectrum Analysis. *Proc. IEEE* 57, 1408-1418. doi:10.1109/proc.1969.7278

Celikbilek, A., and Sapmaz, G. (2016). Risk Management and Microzonation in Urban Planning: an Analysis for Istanbul. *Disaster Sci. Eng.* 2 (2), 59-66.

Crespellani, T. (2014). Seismic Microzoning in Italy: a Brief History and Recent Experiences. *Ing. Sism* 31 (2), 3-31.

## ACKNOWLEDGMENTS

The authors thank the colleagues of INGV-Milano for the cooperation during the geophysical surveys, in particular Antonucci A., Manganello P., Mirenna S., Pessina V., Rizzo A. and Rovida A. Particular thanks to the Associate Editor, Roberto Paolucci, and to the reviewers who provided helpful reviews that improved this paper.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2022.915083/full#supplementary-material

D'Amico, V., Picozzi, M., Baliva, F., Albarello, D., Menichetti, M., Bozzano, F., et al. (2006). "Test Sites in Europe for the Evaluation of Ground Motion Amplification: Site Response of the Gubbio Basin (Central Italy) Using Geological Data and Seismic Noise Measurements," in Proceeding of the 1st European Conference on Earthquake Engineering and Seismology, Geneva, Switzerland (Geneva: Abstract Book), 301.

D'Amico, V., Picozzi, M., Baliva, F., and Albarello, D. (2008). Ambient Noise Measurements for Preliminary Site-Effects Characterization in the Urban Area of Florence, Italy. *Bull. Seismol. Soc. Am.* 98 (3), 1373-1388. doi:10.1785/0120070231

Darendeli, M. B. (2001). *Development of a New Family of Normalized Modulus Reduction and Material Damping Curves*. (PhD thesis). Austin, Texas: University of Texas. Available at: https://repositories.lib.utexas.edu/bitstream/handle/2152/10396/darendelimb016.pdf March, accessed 2022).

De Ferrari, R., Ferretti, G., Barani, S., and Spallarossa, D. (2010). Investigating on the 1920 Garfagnana Earthquake (Mw=6.5): Evidences of Site Effects in Villa Collemandina (Tuscany, Italy). *Soil Dyn. Earthq. Eng.* 30, 1417-1429. doi:10.1016/j.soildyn.2010.07.004

De Franco, R., Biella, G., Caielli, G., Berra, F., Guglielmin, M., Lozej, A., et al. (2009). Overview of High Resolution Seismic Prospecting in Pre-Alpine and Alpine Basins. *Quat. Int.* 204 (1-2), 65-75. doi:10.1016/j.quaint.2009.02.011

Dondi, L., Mostardini, F., and Rizzini, A. (1982). "Evoluzione sedimentaria e paleogeografia nella pianura padana," in *Guida alla geologia del margine appenninico-padano, Guida Geol. Reg. S.G.I.* Editors G. Cremonini and F. Ricci Lucchi (Bologna (Italy), 47-58.

Eurocode 8, CEN (2004). Eurocode 8: Design of Structures for Earthquake Resistance. P1: General Rules, Seismic Actions and Rules for Buildings. *Draft* 6, 522. Doc CEN/TC250/SC8/N335.

Fäh, D., Kind, F., and Giardini, D. (2001). A Theoretical Investigation of Average H/V Ratios. *Geophys. J. Int.* 145 (2), 535-549. doi:10.1046/j.0956-540x.2001.01406.x

Fäh, D., Kind, F., and Giardini, D. (2003). Inversion of Local S-Wave Velocity Structures from Average H/V Ratios, and Their Use for the Estimation of Site-Effects. *J. Seismol.* 7 (4), 449-467. doi:10.1023/b:jose.0000005712.86058.42

Ferretti, G., Zunino, A., Scafidi, D., Barani, S., and Spallarossa, D. (2013). On Microseisms Recorded Near the Ligurian Coast (Italy) and Their Relationship with Sea Wave Height. *Geophys. J. Int.* 194, 524-533. doi:10.1093/gji/ggt114

Furumura, T., and Kennett, B. L. N. (1998). On the Nature of Regional Seismic Phases-III. The Influence of Crustal Heterogeneity on the Wavefield for Subduction Earthquakes: the 1985 Michoacan and 1995 Copala, Guerrero, Mexico Earthquakes. *Geophys. J. Int.* 135 (3), 1060-1084. doi:10.1046/j.1365-246x.1998.00698.x

GeoMol Team(2015). in *GeoMol – Assessing Subsurface Potentials of the Alpine Foreland Basins for Sustainable Planning and Use of Natural Resources – Project Report*. Lf. U. Augsburg, 188.

Grazier, V., Shakal, A., Scrivner, C., Hauksson, E., Polet, J., and Jones, L. (2002). TriNet Strong-Motion Data from the M 7.1 Hector Mine, California, Earthquake of 16 October 1999. *Bull. Seismol. Soc. Am.* 92, 1525-1542. doi:10.1785/0120000925

Hailemikael, S., Amoroso, S., and Gaudiosi, I. (2020). Guest editorial: Seismic Microzonation of Central Italy Following the 2016-2017 Seismic Sequence. *Bull. Earth. Eng.* 18, 5415-5422. doi:10.1007/s10518-020-00929-6
Hanks, T. C. (1975). Strong Ground Motion of the San Fernando, California, Earthquake: Ground Displacements. *Bull. Seismol. Soc. Am.* 65 (1), 193–225. doi:10.1785/bssa0650010193

Hatayama, K. (2008). Lessons from the 2003 Tokachi-Oki, Japan, Earthquake for Prediction of Long-Period Strong Ground Motions and Sloshing Damage to Oil Storage Tanks. *Journal of Seismol* 12, 255–263. doi:10.1007/s10950-007-9066-y

Hisada, Y., Aki, K., and Teng, T. L. (1993). 3-D Simulations of Surface Wave Propagation in the Kanto Sedimentary Basin, Japan Part 2: Application of the Surface Wave BEM. *Bull. Seismol. Soc. Am.* 83 (6), 1700–1720. doi:10.1785/BSSA0830061700

Idriss, I. M., and Sun, J. I. (1993). *User's Manual for SHAKE91: A Computer Program for Conducting Equivalent Linear Seismic Response Analyses of Horizontally Layered Soil Deposits*. Davis, CA, USA: Center for Geotechnical Modeling, Department of Civil and Environmental Engineering, University of California.

ISPRAGeologico d'Italia, S., and Lombardia, R. (2016). in *Carta Geologica d'italia alla scala 1:50.000, Foglio 118 Milano (con note illustrative)*. Editors V. A cura di Francani, A. Piccin, M. Credali, F. Berra, D. Battaglia, P. Gattinoni, et al.

Joyner, W. B. (2000). Strong Motion from Surface Waves in Deep Sedimentary Basins. *Bull. Seismol. Soc. Am.* 90 (6B), S95–S112. doi:10.1785/0120000505

Kagawa, T., Zhao, B., Miyakoshi, K., and Irikura, K. (2004). Modeling of 3D Basin Structures for Seismic Wave Simulations Based on Available Information on the Target Area: Case Study of the Osaka Basin, Japan. *Bull. Seismol. Soc. Am.* 94 (4), 1353–1368. doi:10.1785/012003165

Kim, H.-S., Sun, C.-G., and Cho, H.-I. (2017). Geospatial Big Data-Based Geostatistical Zonation of Seismic Site Effects in Seoul Metropolitan Area. *Int. J. Geo-Information* 6 (6), 174. doi:10.3390/ijgi6060174

Klin, P., Laurenzano, G., Romano, M. A., Priolo, E., and Martelli, L. (2019). ER3D: A Structural and Geophysical 3-D Model of Central Emilia-Romagna (Northern Italy) for Numerical Simulation of Earthquake Ground Motion. *Solid earth*. 10 (3), 931–949. doi:10.5194/se-10-931-2019

Koketsu, K., Hatayama, K., Furumura, T., Ikegami, Y., and Akiyama, S. (2005). Damaging Long-Period Ground Motions from the 2003 Mw 8.3 Tokachi-Oki, Japan Earthquake. *Seismol. Res. Lett.* 76, 67–73. doi:10.1785/gsslrl.76.1.67

Koketsu, K., and Miyake, H. (2008). A Seismological Overview of Long-Period Ground Motion. *J. Seismol.* 12 (2), 133–143. doi:10.1007/s10950-007-9080-0

Konno, K., and Ohmachi, T. (1998). Ground-motion Characteristics Estimated from Spectral Ratio between Horizontal and Vertical Components of Microtremor. *Bull. Seismol. Soc. Am.* 88, 228–241. doi:10.1785/bssa0880010228

Lacoss, R. T., Kelly, E. J., and Toksöz, M. N. (1969). Estimation of Seismic Noise Structure Using Arrays. *Geophysics* 34, 21–38. doi:10.1190/1.1439995

Lai, C. G., Poggi, V., Famà, A., Zuccolo, E., Bozzoni, F., Meisina, C., et al. (2020). An Inter-disciplinary and Multi-Scale Approach to Assess the Spatial Variability of Ground Motion for Seismic Microzonation: the Case Study of Cavezzo Municipality in Northern Italy. *Eng. Geol.* 274, 105722. doi:10.1016/j.enggeo.2020.105722

Laurenzano, G., Priolo, E., Mucciarelli, M., Martelli, L., and Romanelli, M. (2017). Site Response Estimation at Mirandola by Virtual Reference Station. *Bull. Earthq. Eng.* 15 (6), 2393–2409. doi:10.1007/s10518-016-0037-y

Lermo, J., and Chávez-García, F. J. (1993). Site Effect Evaluation Using Spectral Ratios with Only One Station. *Bull. Seismol. Soc. Am.* 83 (5), 1574–1594. doi:10.1785/bssa0830051574

Liang, J. Z., Hao, H., Wang, Y., and Bi, K. M. (2009). Design Earthquake Ground Motion Prediction for Perth Metropolitan Area with Microtremor Measurements for Site Characterization. *J. Earthq. Eng.* 13 (7), 997–1028. doi:10.1080/13632460802687710

Liu, L., Chen, Q.-f., Wang, W., and Rohrbach, E. (2014). Ambient Noise as the New Source for Urban Engineering Seismology and Earthquake Engineering: a Case Study from Beijing Metropolitan Area. *Earthq. Sci.* 27 (1), 89–100. doi:10.1007/s11589-013-0052-x

Luzi, L., Pacor, F., Ameri, G., Puglia, R., Burrato, P., Massa, M., et al. (2013). Overview on the Strong-Motion Data Recorded during the May-June 2012 Emilia Seismic Sequence. *Seismol. Res. Lett.* 84 (4), 629–644. doi:10.1785/0220120154

Maddalena Michele, M., Raffaele Di Stefano, R., Lauro Chiaraluce, L., Marco Cattaneo, M., Pasquale De Gori, P., Giancarlo Monachesi, G., et al. (2016). The Amatrice 2016 Seismic Sequence: a Preliminary Look at the Mainshock and Aftershocks Distribution. *Ann. Geophys. Fast track* 59, 59. doi:10.4401/ag-7227

Marcucci, S., Milana, G., Hailemikael, S., Carlucci, G., Cara, F., Di Giulio, G., et al. (2019). The Deep Bedrock in Rome, Italy: a New Constraint Based on Passive

Seismic Data Analysis. *Pure Appl. Geophys.* 176 (6), 2395–2410. doi:10.1007/s00024-019-02130-6

Martelli, L., and Romani, M. (2013). *Microzonazione Sismica e analisi della condizione limite per l'emergenza delle aree epicentrali dei terremoti della pianura emiliana di maggio-giugno 2012, relazione illustrativa, Servizio geologico, sismico e dei suoli Regione Emilia Romagna (in Italian).*

Marzorati, S., and Bindi, D. (2006). Ambient Noise Levels in North Central Italy. *G-Cube* 7 (9), 1–14. doi:10.1029/2006gc001256

Mascandola, C., Massa, M., Barani, S., Lovati, S., and Santulin, M. (2017). Long-Period Amplification in Deep Alluvial Basins and Consequences for Site-Specific Probabilistic Seismic-Hazard Analysis: An Example from the Po Plain (Northern Italy). *Bull. Seismol. Soc. Am.* 107 (2), 770–786. doi:10.1785/0120160166

Mascandola, C., Massa, M., Barani, S., Albarello, D., Lovati, S., Martelli, L., et al. (2019). Mapping the Seismic Bedrock of the Po Plain (Italy) through Ambient-Vibration Monitoring. *Bull. Seismol. Soc. Am.* 109 (1), 164–177. doi:10.1785/0120180193

Mascandola, C., Barani, S., Massa, M., and Albarello, D. (2021). New Insights on Long-Period (>1s) Seismic Amplification Effects in Deep Sedimentary Basins: the Case of the Po Plain (Northern Italy). *Bull. Seism. Soc. Am.* 111 (4), 2071–2086. doi:10.1785/0120200315

Massa, M., and Augliera, P. (2013). Teleseisms as Estimators of Experimental Long-Period Site Amplification: Application to the Po Plain (Italy) for the 2011 Mw 9.0 Tohoku-Oki (Japan) Earthquake. *Bull. Seismol. Soc. Am.* 103 (5), 2541–2556. doi:10.1785/0120120164

Milana, G., Bordoni, P., Cara, F., Di Giulio, G., Hailemikael, S., and Rovelli, A. (2013). 1D Velocity Structure of the Po River Plain (Northern Italy) Assessed by Combining Strong Motion and Ambient Noise Data. *Bull. Earthq. Eng.* 12 (5), 2195–2209. doi:10.1007/s10518-013-9483-y

NTC (2018). *Aggiornamento delle Norme Tecniche per le Costruzioni. Part 3.2.2: Categorie di sottosuolo e condizioni topografiche*, 42. Rome: Gazzetta Ufficiale n. del 20 febbraio 2018 (in Italian).

Miyazaki, S. I., Segall, P., Fukuda, J., and Kato, T. (2004). Space Time Distribution of Afterslip Following the 2003 Tokachi-oki Earthquake: Implications for Variations in Fault Zone Frictional Properties. *Geophys. Res. Lett.* 31 (6), 410. doi:10.1029/2003gl019410

Molnar, S., Assaf, J., Sirohey, A., and Adhikari, S. R. (2020). Overview of Local Site Effects and Seismic Microzonation Mapping in Metropolitan Vancouver, British Columbia, Canada. *Eng. Geol.* 270, 105568. doi:10.1016/j.enggeo.2020.105568

Molnar, S., Sirohey, A., Assaf, J., Bard, P. Y., Castellaro, S., Cornou, C., et al. (2022). A Review of the Microtremor Horizontal-To-Vertical Spectral Ratio (MHVSR) Method. *J. Seismol* 16, 1–33. doi:10.1007/s10950-021-10062-9

Moya, A., Schmidt, V., Segura, C., Boschini, I., and Atakan, K. (2000). Empirical Evaluation of Site Effects in the Metropolitan Area of San José, Costa Rica. *Soil Dyn. Earthq. Eng.* 20 (1-4), 177–185. doi:10.1016/s0267-7261(00)00049-x

MPS Working Group (2004). *Redazione della mappa di pericolosità sismica prevista dall'Ordinanza PCM del 20 marzo 2003 Rapporto Conclusivo per il Dipartimento della Protezione Civile.*

Mucciarelli, M., Gallipoli, M. R., Di Giacomo, D., Di Nota, F., and Nino, E. (2005). The Influence of Wind on Measurements of Seismic Noise. *Geophys. J. Int.* 161 (2), 303–308. doi:10.1111/j.1365-246x.2004.02561.x

Muttoni, G., Carcano, C., Garzanti, E., Ghielmi, M., Piccin, A., Pini, R., et al. (2003). Onset of Major Pleistocene Glaciations in the Alps. *Geol* 31 (11), 989–992. doi:10.1130/g19445.1

Nakamura, Y. (1989). A Method for Dynamic Characteristics Estimation of Subsurface Using Microtremor or the Ground Surface. *QR Railw. Tech. Res. Inst.* 30 (1), 25–33.

Ohori, M., Nobata, A., and Wakamatsu, K. (2002). A Comparison of ESAC and FK Methods of Estimating Phase Velocity Using Arbitrarily Shaped Microtremor Arrays. *Bull. Seismol. Soc. Am.* 92, 2323–2332. doi:10.1785/0119980109

Okada, H. (2003). *The Microtremor Survey Method*. Houston, TX, United States: SEG library.

Paolucci, E., Albarello, D., D'Amico, S., Lunedei, E., Martelli, L., Mucciarelli, M., et al. (2015). A Large Scale Ambient Vibration Survey in the Area Damaged by May-June 2012 Seismic Sequence in Emilia Romagna, Italy. *Bull. Earthq. Eng.* 13 (11), 3187–3206. doi:10.1007/s10518-015-9767-5
Park, C. B., Miller, R. D., and Xia, J. (1998). "Imaging Dispersion Curves of Surface Waves on Multi-Channel Record," in *SEG Technical Program Expanded Abstracts* 1998 (Houston, TX, United States: Society of Exploration Geophysicists), 1377-1380.

Park, C. B., Miller, R. D., and Xia, J. (1999). Multichannel Analysis of Surface Waves. *Geophysics* 64, 800–808. doi:10.1190/1.1444590

Parolai, S., Bindi, D., and Augliera, P. (2000). Application of the Generalized Inversion Technique (GIT) to a Microzonation Study: Numerical Simulations and Comparison with Different Site-Estimation Techniques. *Bull. Seismol. Soc. Am.* 90 (2), 286–297. doi:10.1785/0119990041

Parolai, S., Picozzi, M., Richwalski, S. M., and Milkereit, C. (2005). Joint Inversion of Phase Velocity Dispersion and H/V Ratio Curves from Seismic Noise Recordings Using a Genetic Algorithm, Considering Higher Modes. *Geophys. Res. Lett.* 32, L01303. doi:10.1029/2004GL021115

Parolai, S., Richwalski, S. M., Milkereit, C., and Fäh, D. (2006). S-wave Velocity Profiles for Earthquake Engineering Purposes for the Cologne Area (Germany). *Bull. Earthq. Eng.* 4, 65–94. doi:10.1007/s10518-005-5758-2

Picozzi, M., Strollo, A., Parolai, S., Durukal, E., Özel, O., Karabulut, S., et al. (2009). Site Characterization by Seismic Noise in Istanbul, Turkey. *Soil Dyn. Earthq. Eng.* 29 (3), 469–482. doi:10.1016/j.soildyn.2008.05.007

Pieri, M., and Groppi, F. (1981). "Subsurface Geological Structure of the Po Plain, Italy," in *Progetto Finalizzato Geodinamica/Sottoprogetto "Modello Strutturale"* (Rome: Consiglio Nazionale delle Ricerche Publ. N°), 414.

Poggi, V., and Fäh, D. (2010). Estimating Rayleigh Wave Particle Motion from Three-Component Array Analysis of Ambient Vibrations. *Geophys. J. Int.* 180 (Issue 1), 251–267. doi:10.1111/j.1365-246x.2009.04402.x

Puglia, R., Tokeshi, K., Picozzi, M., D'Alema, E., Parolai, S., and Foti, S. (2011). Interpretation of Microtremor 2D Array Data Using Rayleigh and Love Waves: the Case Study of Bevagna (Central Italy). *Near Surf. Geophys.* 9 (6), 529–540. doi:10.3997/1873-0604.2011031

Priolo, E., Romanelli, M., Barnaba, C., Mucciarelli, M., Laurenzano, G., Dall'Olio, L., et al. (2012). The Ferrara Thrust Earthquakes of May-June 2012: Preliminary Site Response Analysis at the Sites of the OGS Temporary Network. *Ann. Geophy.* 55 (4), 591–597. doi:10.4401/ag-6172

Ratchkovski, N. A., Hansen, R. A., Stachnik, J. C., Cox, T., Fox, O., Rao, L., et al. (2003). Aftershock Sequence of the Mw 7.9 Denali Fault, Alaska, Earthquake of 3 November 2002 from Regional Seismic Network Data. *Seismol. Res. Lett.* 74 (6), 743–752. doi:10.1785/gssrl.74.6.743

Regione Lombardia, Eni Divisione Agip (2002). *Geologia degli acquiferi Padani della Regione Lombardia, a cura di Cipriano Carcano e Andrea Piccin*. Firenze: S.EL.CA. (in Italian).

Ritter, J. R. R., Balan, S. F., Bonjer, K.-P., Diehl, T., Forbriger, T., Marmureanu, G., et al. (2005). Broadband Urban Seismology in the Bucharest Metropolitan Area. *Seismol. Res. Lett.* 76 (5), 574–580. doi:10.1785/gssrl.76.5.574

Rizzini, A., and Dondi, L. (1978). Erosional Surface of Messinian Age in the Subsurface of the Lombardian Plain (Italy). *Mar. Geol.* 27, 303–325. doi:10.1016/0025-3227(78)90037-3

Ronald Abraham, J., Lai, C. G., and Papageorgiou, A. (2015). Basin-effects Observed during the 2012 Emilia Earthquake Sequence in Northern Italy. *Soil Dyn. Earthq. Eng.* 78, 230–242. doi:10.1016/j.soildyn.2015.08.007

Rovida, A., Locati, M., Camassi, R., Lolli, B., and Gasperini, P. (2020). The Italian Earthquake Catalogue CPT115. *Bull. Earthq. Eng.* 18 (7), 2953–2984. doi:10.1007/s10518-020-00818-y

Sambridge, M. (1999). Geophysical Inversion with a Neighbourhood Algorithm-I. Searching a Parameter Space. *Geophys. J. Int.* 138, 479–494. doi:10.1046/j.1365-246x.1999.00876.x

Sato, T., Graves, R. W., and Somerville, P. G. (1999). Three-dimensional Finite-Difference Simulations of Long-Period Strong Motions in the Tokyo Metropolitan Area during the 1990 Odawara Earthquake (MJ 5.1) and the Great 1923 Kanto Earthquake (MS 8.2) in Japan. *Bull. Seismol. Soc. Am.* 89 (3), 579–607. doi:10.1785/bssa0890030579

Scardia, G., De Franco, R., Muttoni, G., Rogledi, S., Caielli, G., Carcano, C., et al. (2012). Stratigraphic Evidence of a Middle Pleistocene Climate-driven Flexural Uplift in the Alps. *Tectonics* 31 (6), 3108. doi:10.1029/2012tc003108

Schnabel, P. B., Lysmer, J., and Seed, H. B. (1993). *SHAKE-91: Equivalent Linear Seismic Response Analysis of Horizontally Layered Soil Deposits*. Berkeley, CA, USA: The Earthquake Engineering Online Archive NISEE E-Library.

Schnabel, P. B. (1972). *SHAKE a Computer Program for Earthquake Response Analysis of Horizontally Layered Sites*. Berkeley, CA, USA: EERC Report; University of California.

Schneider, J. A., Mayne, P. W., and Rix, G. J. (2001). Geotechnical Site Characterization in the Greater Memphis Area Using Cone Penetration Tests. *Eng. Geol.* 62 (1-3), 169–184. doi:10.1016/s0013-7952(01)00060-6

SESAME (2004). Guidelines for the Implementation of H/V Spectral Ratio Technique on Ambient Vibration Measurements, Processing and Interpretation. Available at: http://sesame-fp5.obs.ujf-grenoble.fr/Delivrables/Del-D23-HV_User_Guidelines.pdf Last accessed March, 2022.

Somerville, P. G., Collins, N. F., Graves, R. W., and Pitarka, A. (2004). "An Engineering Ground Motion Motion Model for Basin Generated Surface Waves," in Proc. 13th World Conference on Earthquake Engineering, Vancouver, B.C., Canada, August 2004.

Stephenson, W. J., Hartzell, S., Frankel, A. D., Asten, M., Carver, D. L., and Kim, W. Y. (2009). Site Characterization for Urban Seismic Hazards in Lower Manhattan, New York City, from Microtremor Array Analysis. *Geophys. Res. Lett.* 36 (3), 444. doi:10.1029/2008gl036444

Stucchi, M., Meletti, C., Montaldo, V., Crowley, H., Calvi, G. M., and Boschi, E. (2011). Seismic Hazard Assessment (2003-2009) for the Italian Building Code. *Bull. Seismol. Soc. Am.* 101, 1885-1911. doi:10.1785/0120100130

Sun, C.-G., Kim, H.-S., Chung, C.-K., and Chi, H.-C. (2014). Spatial Zonations for Regional Assessment of Seismic Site Effects in the Seoul Metropolitan Area. *Soil Dyn. Earthq. Eng.* 56, 44-56. doi:10.1016/j.soildyn.2013.10.003

Tarabusi, G., and Caputo, R. (2017). The Use of HVSR Measurements for Investigating Buried Tectonic Structures: the Mirandola Anticline, Northern Italy, as a Case Study. *Int. J. Earth Sci. Geol. Rundsch* 106 (1), 341–353. doi:10.1007/s00531-016-1322-3

Tragni, N., Calamita, G., Lastilla, L., Belloni, V., Ravanelli, R., Lupo, M., Salvia, V., and Gallipoli, M. R. (2021). Sharing Soil and Building Geophysical Data for Seismic Characterization of Cities Using CLARA WebGIS: A Case Study of Matera (Southern Italy). *Appl. Sci.* 11 (9), 4254. doi:10.3390/app11094254

ViDEPI project (2009). *Visibilità dei dati afferenti all'attività di esplorazione petrolifera in Italia Ministero dello sviluppo economico DGRME - Società Geologica Italiana, Assomineraria*. Available at: https://www.videpi.com

Vuan, A., Klin, P., Laurenzano, G., and Priolo, E. (2011). Far-Source Long-Period Displacement Response Spectra in the Po and Venetian Plains (Italy) from 3D Wavefield Simulations. *Bull. Seismol. Soc. Am.* 101 (3), 1055-1072. doi:10.1785/0120090371

Wang, J. H. (2008). Urban Seismology in the Taipei Metropolitan Area: Review and Prospective. *TAO Terr. Atmos. Ocean. Sci.* 19 (3), 2. doi:10.3319/tao.2008.19.3.213(t)

Wathelet, M., Chatelain, J.-L., Cornou, C., Giulio, G. D., Guillier, B., Ohrnberger, M., et al. (2020). Geopsy: A User-Friendly Open-Source Tool Set for Ambient Vibration Processing. *Seism. Soc. Lett.* 91 (3), 1878–1889. doi:10.1785/0220190360

Wathelet, M., Jongmans, D., Ohrnberger, M., and Bonnefoy-Laudet, S. (2008). Array Performances for Ambient Vibrations on a Shallow Structure and Consequences over V S Inversion. *J. Seismol.* 12, 1–19. doi:10.1007/s10950-007-9067-x

Wessel, P., and Smith, W. H. F. (1995). New Version of the Generic Mapping Tools. *EOS Trans. AGU* 76, 329. doi:10.1029/95eo00198

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

**Publisher's Note:** All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Copyright © 2022 Massa, Lovati, Puglia, Brunelli, Lorenzetti, Mascandola, Felicetta, Pacor and Luzi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
<figure><img src="image_81.png" /></figure>

Check for
updates

# Building Stock Classification Using Machine Learning: A Case Study for Oslo, Norway

Federica Ghione<fnref n="1" />,²*, Steffen Mæland<fnref n="3" />, Abdelghani Meslem<fnref n="1" />,⁴ and Volker Oye<fnref n="1" />,²

This paper describes a new concept to automatically characterize building types in urban areas based on publicly available image databases, making parts of seismic risk assessment more time and cost-effective, and improving the reliability of seismic risk assessment, especially in regions where building stock information is currently not documented. One of the main steps in evaluating potential human and economic losses in a seismic risk assessment, is the development of inventory databases for existing building stocks in terms of load-resisting structural systems and material characteristics (building typologies classification). The common approach for building stock model classification is to perform extensive fieldwork and walk-down surveys in representative areas of a city, and in some cases using random sample surveys of geounits. This procedure is time and cost consuming, and subject to personal interpretation: to mitigate these costs, we have introduced a machine learning methodology to automate this classification based on publicly available image databases. We here use a Convolutional Neural Network (CNN) to automatically identify the different building typologies in the city of Oslo, Norway, based on facade images taken from *in-situ* fieldwork and from Google Street View. We use transfer learning of state-of-the-art pretrained CNNs to predict the Model Building Typology. The present article attempts to categorize Oslo's building stock in five main building typologies: timber, unreinforced masonry, reinforced concrete, composite (steel-reinforced concrete) and steel. This method results in 89% accuracy score for timber buildings, though only 35% success score for steel-reinforced concrete buildings. We here classify and define for the first time a relevant set of five typologies for the Norwegian building typologies as observed in Oslo and applicable at national level. In addition, this study shows that CNNs can significantly contribute in terms of developing a cost-effective building stock model.

**Keywords:** building stock model, convolutional neural network, machine learning, seismic risk assessment, Oslo (Norway)

## 1 INTRODUCTION

Within the last century, earthquakes, flooding and droughts have been the dominating phenomena responsible for causing fatalities and economic losses worldwide (see Data Availability Statement section). In the last few decades though, earthquakes are the dominating phenomena standing for most fatalities (Wallemacq and House, 2018), providing the motivation and requirement for more

OPEN ACCESS

**Edited by:**

Roberto Paolucci,
Politecnico di Milano, Italy

**Reviewed by:**

Roberto Gentile,
University College London,
United Kingdom

Dennis Wagenaar,
Nanyang Technological University,
Singapore

*Correspondence:

Federica Ghione
federica.ghione@norsar.no

**Specialty section:**

This article was submitted to
Geohazards and Georisks,
a section of the journal
Frontiers in Earth Science

**Received:** 28 February 2022

**Accepted:** 26 May 2022

**Published:** 16 June 2022

**Citation:**

Ghione F, Mæland S, Meslem A and
Oye V (2022) Building Stock
Classification Using Machine Learning:
A Case Study for Oslo, Norway.
Front. Earth Sci. 10:886145.
doi: 10.3389/feart.2022.886145

<footnote n="1">¹Department of Applied Seismology, NORSAR, Kjeller, Norway, ²Department of Geosciences, University of Oslo, Oslo, Norway,</footnote>

<footnote n="3">³Department of Holistic Systems, SimulaMet, Oslo, Norway, ⁴Norwegian University of Life Sciences (NMBU), Ås, Norway</footnote>
detailed and effective seismic risk assessment studies. A seismic risk assessment estimates the probability of losses if an earthquake occurs, can assist through updated building regulations and to initiate other mitigation actions to avoid likelihood for casualties and economic losses. To do that, information regarding seismic hazard, vulnerability and exposure models of the area are needed (Silva et al., 2014). In particular, the exposure model contains information about all the buildings in a specific area, infrastructures, and population data. A seismic risk assessment for a building stock in a city needs a classification of the buildings in accordance with their structural typologies, an important parameter to define the building performance under seismic load.

The development of a building stock model is always challenging and time consuming, especially when the area of investigation is large. Until today it is a common practice to use both Google Earth and in-situ surveys to get information on building structural systems and material characteristics. However, machine learning methods to analyze visual imagery have recently been applied to classify building stocks, using online available façade images, obtained through Google Street View (GSV), as key information. To our knowledge, this automatic image analysis using Convolutional Neural Network (CNN) was previously only tested in two places to automatically detect building materials and types of lateral-load resistance: Gonzalez et al. (2020) for Medellin city, Colombia and Aravena Pelizari et al. (2021) for Santiago de Chile. Those two examples are similar to our study, in terms of applying CNNs to automatically detect and identify the different building typologies for the two cities using GSV images. In this study, however, we define for the first time a Model Building Typology (MBT) for Oslo, and we conduct a survey amongst experts within earthquake engineering. Based on the survey's results, the MBT identified in Oslo can further be applied at the national scale (Norway). Another example by GFDRR (2018) shows that using a combination of deep learning and GSV together with drone images may lead to remotely identify buildings that require further inspection and possible retrofitting/strengthening for Guatemala city, Guatemala.

This paper presents the building stock model for Oslo, capital of Norway. The general information of the buildings is obtained from the public cadastre from www.kartverket.no. This information contains the total number of buildings in Oslo with the corresponding coordinates, number of stories, number of housing units, usable area, and total area.

To perform seismic risk studies, the availability of detailed building stock typologies is a necessity. This is usually a very time-consuming process, and our innovative approach using transfer learning and CNNs for the classification may provide this required input in much shorter time and at lower cost. This could be a solution for most of the urbanized regions in the world to develop seismic risk assessments and incorporate earthquake preparedness actions.

The paper is organized as follow: in **Section 2** we describe the methodology followed and the CNN approach; **Section 3** describes the study area, the identified building typologies, and the data used for this work; **Section 4** introduces the results and **Section 5** presents the main conclusion.

# 2 METHODOLOGY

The methodology is based on the use of a CNN for automating the typology classification step. While the classification itself requires no manual involvement, the main workload goes into compiling a labelled image dataset (discussed in **Section 3**) and training the parameters of the CNN (described in this section). Once a CNN structure is selected and the CNN is suitably trained, classifying a set of unseen images is done quickly and at low computational cost. A diagram depicting the workflow is shown in **Figure 1**. First, a representative set of façade images must be collected and labelled, which entails the manual part of the process. Using these images to train a CNN, subsequent images of unknown typologies can be downloaded online and classified in bulk and used for seismic risk assessment. There are two challenges when using CNNs for this task: the need for a large set of labelled images, and the computational load of training them. We mitigate both by the use of *transfer learning*, which will be described in this section.

## 2.1 Image Recognition

Current state-of-the-art methods for recognizing objects in images are based on variants of CNNs, which are artificial neural networks that use sliding filters (or kernels) to process their inputs (LeCun et al., 2015). Unlike traditional image processing, these filters are learnt from examples during a training phase, and they are typically stacked so that the output from one convolutional layer is used as input to the next one. A benefit of the convolutional filters is that they are invariant under translations, meaning that the positioning of objects in an image does not matter (Goodfellow et al., 2016). Stacking convolutional layers allows for learning pattern hierarchies, where the first layers recognize simple shapes such as vertical or horizontal lines, while the last layers recognize compositions of these patterns, for instance the shape of a building. In most approaches, the convolutional layers are interspersed with pooling layers, which reduce the dimensions by downsampling. These reduce the required number of learnable parameters, and at the same time introduce invariance to rotation and scaling (Goodfellow et al., 2016). In addition, it is common to add special layers and mechanisms that facilitate the learning process, such as dropout (Srivastava et al., 2014) and residual (skip) connections (He et al., 2016a). The composition of layers is referred to as the *architecture* of the network, while we use *model* to indicate a network with a particular set of optimized parameters. For classification purposes, the output of the final convolution layer is input to one or several fully connected layers, which ultimately output a mutually exclusive prediction for which class an input image belongs to. The learnable parameters of a model include both the convolution filters and the weights of the classification layer.

## 2.2 Transfer Learning and Network Architectures

Modern CNNs typically have millions of free learnable parameters and optimizing them requires large corpora of