<figure><img src="image_82.png" /><figcaption>**FIGURE 1** | Workflow diagram. Using manually annotated façade images from fieldwork, from GSV, or a combination of the two, a CNN is trained to classify building typologies. A large set of images is then downloaded from GSV, and automatically classified to obtain the building stock for the area of interest.</figcaption></figure>

data. In image recognition, the standard dataset for training and performance evaluation is the ImageNet (ILSVRC) database (Russakovsky et al., 2015), which contains nearly 1.3 million labelled images of various common objects. When training a CNN on such a diverse dataset, the first layers will be sensitive to simple shapes like straight or curved lines, while the middle layers are sensitive to different compositions of these basic patterns. Since only the last layers are highly specific to the particular dataset the CNN is trained on, applications such as ours can benefit from transfer learning, where one first trains the CNN on a large, generic dataset such as ImageNet, and subsequently re-trains (or fine-tunes) the last layers on the application-specific data. This allows for re-using the knowledge contained in the first and middle layers and is particularly useful in cases like ours where the application-specific dataset is comparatively small.

To select the optimal CNN architecture for our case, we compute four performance metrics (described in next paragraph) for eleven candidates, all pre-trained on ImageNet data. These are available in the Keras framework (Chollet, 2015). For each candidate, we set all parameters to remain constant, but replace the final classification layer by a new layer with six output nodes, corresponding to the number of building typologies under consideration. As a measure against overfitting, we apply dropout (Srivastava et al., 2014) before the classification layer. This means temporarily removing nodes from the network during training, which improves robustness by reducing the nodes’ reliance on each other. We randomly drop nodes at a 20% probability. We then train on our data, updating only the parameters of the classification layer. We perform 4-fold cross-validation, meaning we divide the training data into four equally sized parts, and use three parts for training and one for computing metrics. The parts are then rotated so that all parts are used for computing metrics. Since the weights of the final layer are randomly initialized for each training round, potentially affecting the final performance, we repeat the cross-validation process twice, and report the

average results for the in total eight trained models per architecture. Training is done for up to 100 epochs, where an epoch is a single pass over all the training data, followed by a metric evaluation. If learning fails to improve over five consecutive epochs, the learning rate is reduced by half, thereby taking shorter steps toward the optimal solution. If learning fails to improve over ten consecutive epochs, the training is stopped.

The four computed metrics are: accuracy, precision, recall, and area under the receiver operating characteristic curve. Accuracy is defined as the fraction of samples where the predicted class exactly matches the true class; in our case of 6 distinct classes, a random classifier would have an accuracy of $1/6 = 0.17$, assuming that the number of samples of each class are the same. Precision is defined as the number of true positives divided by the sum of true positives and false positives, indicating the quality of the positive predictions. Recall, on the other hand, is defined as the number of true positives divided by the sum of true positives and false negatives, indicating the completeness of the positive predictions. In addition, we compute the area under the receiver operating curve (ROC), which is obtained by plotting the true positive rate against the false positive rate. Since true and false predictions relate to a binary classification problem, we compute precision, recall and ROC individually for each class, in a one-vs-all fashion, and then report the average across all classes. The results are listed in **Table 1**. For all metrics, the best possible result is 1, while the worst possible result is 0.

## 2.3 Fine-Tuning

From **Table 1** we observe that most CNN architectures yield similar results, but the three members of the DenseNet (Huang et al., 2017) family stand out as the highest performers across all evaluation metrics. They are structurally similar but differ in the number of convolution layers, which is given by the number at the end of their names (121, 169 and 201, respectively). The DenseNet201 architecture excels on all metrics and is therefore
<table><caption>**TABLE 1** | Metrics for different CNN architectures. Highlighted in bold are the best observed values. The last column indicates the total number of parameters for each architecture.</caption><thead><tr><th>Architecture</th><th>Accuracy</th><th>Precision</th><th>Recall</th><th>ROC AUC</th><th>Parameters</th></tr></thead><tbody><tr><td>Xception Chollet. (2017)</td><td>0.759</td><td>0.762</td><td>0.759</td><td>0.920</td><td>20,873,774</td></tr><tr><td>VGG16 Simonyan and Zisserman. (2015)</td><td>0.740</td><td>0.739</td><td>0.740</td><td>0.906</td><td>14,717,766</td></tr><tr><td>VGG19 Simonyan and Zisserman. (2015)</td><td>0.732</td><td>0.733</td><td>0.732</td><td>0.904</td><td>20,027,462</td></tr><tr><td>ResNet50V2 He et al. (2016b)</td><td>0.759</td><td>0.757</td><td>0.759</td><td>0.917</td><td>23,577,094</td></tr><tr><td>ResNet101V2 He et al. (2016b)</td><td>0.762</td><td>0.760</td><td>0.762</td><td>0.917</td><td>42,638,854</td></tr><tr><td>ResNet152V2 He et al. (2016b)</td><td>0.757</td><td>0.757</td><td>0.757</td><td>0.912</td><td>58,343,942</td></tr><tr><td>InceptionV3 Szegedy et al. (2016)</td><td>0.762</td><td>0.762</td><td>0.762</td><td>0.922</td><td>21,815,078</td></tr><tr><td>InceptionResNetV2 Szegedy et al. (2017)</td><td>0.750</td><td>0.752</td><td>0.750</td><td>0.912</td><td>54,345,958</td></tr><tr><td>DenseNet121 Huang et al. (2017)</td><td>0.764</td><td>0.764</td><td>0.764</td><td>0.922</td><td>7,043,654</td></tr><tr><td>DenseNet169 Huang et al. (2017)</td><td>0.770</td><td>0.772</td><td>0.770</td><td>0.926</td><td>12,652,870</td></tr><tr><td>DenseNet201 Huang et al. (2017)</td><td><strong>0.778</strong></td><td><strong>0.780</strong></td><td><strong>0.778</strong></td><td><strong>0.932</strong></td><td>18,333,510</td></tr></tbody></table>

selected for further optimization. While being largest in the DenseNet family, we note that this family of architectures is the smallest in terms of number of learnable parameters, compared to the other ones tested. This indicates that raw model size, as listed in the last column of **Table 1**, cannot be directly considered as a proxy for performance, and comparisons are required to find the best architecture for a given application. We find, however, that there is an approximately linear relationship between the model size and the training time. Compared to DenseNet201, the smallest architecture (DenseNet121) takes 12% shorter time to train a single epoch, while the largest architecture (ResNet151V2) takes 31% longer, when training on a Nvidia V100 (NVIDIA, 2020) graphics card.

Proceeding with the DenseNet201 architecture, we fine-tune it to our data in a two-step procedure. Using the full training dataset, apart from 20% of images that are set aside as validation data to monitor the training progress, we first train only the final classification layer of the model. Like before, we halve the learning rate if training has stagnated for five epochs, and end training if it has still not improved after ten epochs. In the second step, we now additionally unfreeze the parameters of the last 64 convolution layers. This allows for adapting a larger part of the model specifically to our data, while still retaining the low-level pattern recognition of the first convolution layers. The DenseNet201 architecture is organized into five blocks of multiple convolution layers, where each block has a complicated internal structure, but has a single connection to the next block. Therefore, we only consider it useful to free the parameters of entire blocks at the time; freeing the last 64 convolution layers equals the entire last block. Doing so means we now train 38% of all parameters in the model (approximately 7 million out of 18 million). Experiments with freeing a larger number of blocks resulted in overfitting, where the accuracy on the validation data diverges greatly from that of the training data. As for the first step, training is run until improvement has failed to improve for ten epochs.

## 3 DATA

In order to conduct seismic risk assessment a building stock model needs to be available, and the key input data to develop

<figure><img src="image_83.png" /><figcaption>**FIGURE 2** | Building distribution (shown with orange dots) in Oslo municipality (study area).</figcaption></figure>

such building stock models are the building typologies as extracted in the relevant region. Building typologies are related to the ability of a building to resist lateral loads that mainly affects the structural system and depends on material and height. The lateral load-resisting system and its material can be identified only from the blueprints (two-dimensional set of technical/engineering drawings that specify a building's dimensions, construction materials, and the exact placement of all its components) or by direct expert observations. Unfortunately, structural blueprints are not always available. Therefore, expert
<figure><img src="image_84.png" /><figcaption>**FIGURE 3** | In the figure, the fieldwork area is represented with dashed line.</figcaption></figure>

opinions seem to be the best option to assemble building inventories and fundamental information. However, in most of the cases, it is timewise and economically not possible to survey each asset. For this reason, many assumptions need to be made to establish a building stock model and interpolate data from surveyed areas to not surveyed neighborhoods. Recently, using GSV (Google Street View) and Google Earth has become an alternative way to carry out fieldwork for visualizing façade and material information remotely (GFDRR, 2018; Kang et al., 2018; Gonzalez et al., 2020; Aravena Pelizari et al., 2021). However, uncertainties are remaining in the building typology classification. In this paper, we use the EMS-98 building taxonomy (Grünthal, 1998) with some extension related to specific typology that are not present in the original classification [e.g. composite (steel-reinforced concrete) typology].

## 3.1 Study Area

Our study area is the city of Oslo, the capital of Norway (Figure 2). The oldest settlements in the area are from around 11000 BC, but Oslo was only founded in 1000 AC. The oldest part of the city is the eastern part, called Bjørvika in Gamle Oslo (see Figure 3). In 1624, a disastrous fire destroyed most of the city, and built on the ashes, the new town was called Christiania. After that devastating fire, a ban on wooden houses was introduced allowing only solid brick, and half-timbered brick houses (Eriksson et al., 2016). In 1769 Christiania had about 7500

inhabitants. During the 1600s the demand for wood rose, and Christiania became an important harbour for trading wood. Outside the city centre, many small wooden houses were built, and some of those still exist today. During the industrial period of the 1840s, many factories were built along the river Akerselva and the population increased significantly to about 113.000 inhabitants. In the late 1880s, many multi-stories brick tenements were constructed to fulfil the request for more living space. In 1925 the city changed its name back to Oslo and the Ring Road was introduced; in 1948 another important step for the city was the development of the new subway system (Eriksson et al., 2016). In 2020 Oslo had a population count of 697,549 and by 2040 the expected number is around 926,000 (Eriksson et al., 2016). This prediction has a direct consequence for the city, i.e. it will need to accommodate the new residents with new buildings. Oslo municipality covers an area of 480 km² and it is subdivided in 17 boroughs or bydel in Norwegian.

One key element to consider during the planning of a city expansion with respect to new infrastructures is an accurate evaluation of the seismic risk, which is strongly depending on the local seismic hazard conditions, in turn influenced by soil amplification. Norway shows low to medium seismicity that leads to a lower level of seismic hazard compared to other Southern European countries (Danciu et al., 2021). Compared on a national scale, the city of Oslo falls within a zone of intermediate seismic hazard. Oslo was hit by a significant earthquake (5.4 Mw) on the 23rd of October 1904 (Bungum et al., 2009). Although the epicenter was located 115 km south of Oslo, the event (known as 1904 Oslofjord earthquake) generated ground motions that propagated on both sides of the Oslo fjord from the south of Fredrikstad/Tønsberg to the north of Oslo. The earthquake was felt over an area of 800.000 km² from Namsos in the north to Poland, and across southern Norway to Helsinki in the east (Bungum et al., 2009). The maximum intensity on the Mercalli scale in Oslo was reported to VI and major damages were reported, mainly for wooden and unreinforced masonry buildings in the Oslo area. Although the severe damages to the buildings, no relevant and significant mitigation measures were put in place after this event. The first seismic standard introduced in Norway was by the end of the 1980s, and it was mostly destined for seismic design of offshore structures. Within the same period there was also another document destined for buildings, but it was introduced and used only as recommendation. In 2004 a new standard on the design of structures considering loads from seismic influence was adopted (NS 3491-12). In 2008, Norway adopted the Eurocode and since then it has been the only standard for seismic design for all types of structures and infrastructures.

## 3.2 Model Building Typology

Building stocks usually vary from locality to locality and certainly even more between different countries due to different construction practices, material availability and construction period.

The following steps are followed to recognize the different building typologies for Oslo and to develop a Model Building Typology (MBT):
Timber (T)

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

Unreinforced Masonry (MUR)

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

Reinforced Concrete (CR)

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

Composite (SRC)

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

Steel (S)

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

Other

<figure><img src="image_90.png" /><figcaption>**FIGURE 4** | Building typologies for the building stock in Oslo: timber (T), unreinforced masonry (MUR), reinforced concrete (CR), composite (SRC), steel (S) and other.</figcaption></figure>

1) a first overview of the building typologies is obtained from Google Earth for the different neighborhoods in the city. The building stock is observed at different scales (from small to large):

- the inbuilt 3D-building option allows us to identify the lateral load bearing systems, predominant work material for walls, façade decoration, flooring and roofing types and number of stories.

- The Street View mode gives us a ground-level view of the different buildings and it helps us to identify similarity and heterogeneity in the same geosit.

2) With the plan view, we understand the general distribution of the building stock in the city looking from the top. It allows us to identify different roof materials and characteristics that can be potentially linked to different building typologies. Through the first evaluation using Google Earth, some area of the city center (Grünerløkka, Oslo central station, Bjørvika) and Alna (see **Figure 3**) are chosen for detailed *in-situ* fieldwork. Those areas contain a good representation of all the typologies and are good candidates to test the machine learning methodology to automatically identify different building typologies.

3) During 5 days of fieldwork during the winter of 2021, about 350 pictures of the facade buildings are taken manually, and information related to structural system and material related data are collected (e.g. lateral loadbearing systems observed, material type, flooring/roofing system, number of stories, usage/type of activity). These form key information to define a MBT.

4) After steps 1 and 2, we define an initial building typology classification, and we divide the observations into five groups with the corresponding typology.

5) To confirm and validate the initial evaluation of the typologies, a survey questionnaire regarding seismic

vulnerability assessment of the existing building stock in Oslo is sent out to experts (mainly engineers working in Norway) in order to validate the building typologies that they have been identified. The survey was divided into eleven different sections, with questions related to the structural system and material characteristics of the building stock in Oslo, date of practice of a given typology and general practices. The expert's opinion confirmed the initial evaluation and they agreed that the building typologies recognized in Oslo are applicable at national scale. The questionnaire shows results compatible to the preliminary assessment.

6) Combining results from both fieldwork and survey questionnaire, a final Model Building Typology (MBT) is defined, and it represents the existing building typologies in Oslo city and Norway.

7) Images from GSV are downloaded for all the districts of the city. A total of 5074 pictures are manually labelled using the MBT previously defined. The pictures are used to train the Convolutional Neural Network (CNN). As Gonzalez et al. (2020), Aravena Pelizari et al. (2021), also here the CNN methodology is applied to automatically detect the different building typologies in Oslo.

The results of the survey questionnaire combined with the field survey have shown that the existing building stock in Oslo can be divided into three main construction periods:

1) Buildings built before 1950, estimated to represent roughly 35% of the total existing building stock, and mostly made of timber and unreinforced masonry.

2) Buildings built between 1950 and 1998, estimated to represent 25% of the total number of buildings, and they are mainly made of reinforced concrete.
<table><caption>**TABLE 2** | Confusion matrix computed on the test data set. Rows show the true class, and the columns the CNN-identified class. Correct predictions follow the diagonal, highlighted in bold. As an example, 73% of CR buildings are correctly identified, but 3% are classified as MUR, 3% are classified as S, 2% SRC, 9% as T and 10% as other. In parentheses are the actual number of images.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="6">Predicted</th></tr><tr><th>CR</th><th>MUR</th><th>S</th><th>SRC</th><th>T</th><th>Other</th></tr></thead><tbody><tr><td rowspan="6">True</td><td>CR</td><td><strong>0.73 (86)</strong></td><td>0.03 (4)</td><td>0.03 (3)</td><td>0.02 (2)</td><td>0.09 (11)</td><td>0.10 (12)</td></tr><tr><td>MUR</td><td>0.05 (5)</td><td><strong>0.71 (76)</strong></td><td>0.00 (0)</td><td>0.00 (0)</td><td>0.20 (21)</td><td>0.05 (5)</td></tr><tr><td>S</td><td>0.18 (5)</td><td>0.00 (0)</td><td><strong>0.57 (16)</strong></td><td>0.00 (0)</td><td>0.11 (3)</td><td>0.14 (4)</td></tr><tr><td>SRC</td><td>0.60 (12)</td><td>0.00 (0)</td><td>0.05 (1)</td><td><strong>0.35 (7)</strong></td><td>0.00 (0)</td><td>0.00 (0)</td></tr><tr><td>T</td><td>0.01 (3)</td><td>0.01 (6)</td><td>0.00 (0)</td><td>0.00 (0)</td><td><strong>0.89 (424)</strong></td><td>0.09 (42)</td></tr><tr><td>Other</td><td>0.02 (4)</td><td>0.00 (1)</td><td>0.01 (2)</td><td>0.00 (0)</td><td>0.12 (32)</td><td><strong>0.87 (232)</strong></td></tr></tbody></table>

3) Buildings built between 1998 and up to today, estimated to cover 40% and they are found as reinforced concrete, composite (steel-reinforced concrete) and steel.

The following typologies are identified in Oslo and can be applied at national scale (see **Figure 4**):

* Timber (T): the construction with timber has been in practice for many years (before 1950), and the timber buildings mostly are frame structures consisting of wooden frames with solid or plate timber members. This typology is mainly used for housing and typically can be one to two stories, but there are also some 3 to 4 story-high buildings that can also be found as residential apartments and public buildings for various activities, e.g., schools. The timber construction is still in practice and starting from 2000 there have been new modern timber buildings that can also be found as high-rise buildings with more than 8 stories. As per the situation today, timber buildings represent around 20% of the existing building stock in Oslo.

* Unreinforced masonry (MUR): unreinforced masonry constructions are load-bearing wall system structures made of material that vary depending on the construction's age, and which can be burnt clay bricks, stone masonry blocks, concrete blocks or mixtures of other materials. Most of the existing buildings of this type were constructed before 1950. The majority can be found as low/ mid-rise structures with 2-5 stories, and very few with more than 6 stories. They are used for residential, commercial and other general activities. In terms of location, this type of construction mostly can be found in the center of the city, and as a rough estimate, they represent 30% of the current building stock in Oslo.

* Reinforced concrete (RC): the construction of reinforced concrete buildings has been in practice since 1950, and the construction procedure of this typology consists of reinforced concrete frames (columns and beams), cast-in-place. But starting from 2000, the pre-cast practice has become more frequent in the construction of this type of structures. The existing RC buildings can cover different ranges of height classes (low, mid and high-rise) and can be

found in most of the districts and zones of the city. It is estimated that this type of buildings represents 35% of the total existing buildings in the city.

* Composite (steel-reinforced concrete) (SRC): the composite construction has been in practice since 2000s, and it is estimated that it represents almost 10% of the total existing building stock of Oslo. Typically, the number of stories for this type of building can range between 2 to more than 15 stories for a few of the most recent buildings. The construction procedure of this typology consists of steel frames and cast-in-place concrete frames (columns and beams) and/or concrete shear walls. Buildings of this typology are mostly used as residential apartments, offices, commercial activities, and they are found in most of the districts and zones of the city.

* Steel (S): the existing steel structures are load-bearing steel moment frame constructions, and mostly are found outside the city center. This type of structures is, in general, used for industrial activities, also as big grocery stores, supermarkets, malls, parking lots, and hangars. Due to the nature of the utilization, story heights can be up to 6 m, and the number of stories can be up to 2 stories. The construction in steel has recently increased with the urban development, and currently it is estimated that this type of construction covers about 5% of the existing buildings in the city of Oslo.

In addition to the above typologies, one extra category "other" is added to account for the case that the algorithm could not recognize the typology or that no building could be identified in the picture at all.

## 3.3 Image Dataset

The dataset of all buildings in Oslo is obtained from the public cadastre (downloaded in December 2020) and it contains the number and the coordinates of all free-standing buildings of more than 50 m². In addition to façade photographs collected manually, images for each building position are downloaded automatically using the Google Street View API (GSV). The Street View service provides near-continuous street level imagery of most of the world's cities, and the service's API allows for direct download of images for a given set of coordinates. In order to accept an image to be related to a building, we set a threshold that the image needs to be taken within 40 m from the requested building coordinate. For the
A

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

true: T, predicted: T

B

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

true: SRC, predicted: CR

<caption>**FIGURE 5** | Examples of correctly (A) and falsely (B) classified images.</caption>

<table><caption>**TABLE 3** | Distribution in terms of numbers and percentages of the predicted building typologies.</caption><thead><tr><td>Typologies</td><td>Number of buildings</td><td>%</td></tr></thead><tbody><tr><td>T</td><td>56,305</td><td>56.7</td></tr><tr><td>MUR</td><td>7,979</td><td>8.0</td></tr><tr><td>CR</td><td>7,536</td><td>7.6</td></tr><tr><td>SRC</td><td>399</td><td>0.4</td></tr><tr><td>S</td><td>515</td><td>0.5</td></tr><tr><td>other</td><td>26,605</td><td>26.8</td></tr></tbody></table>

134,432 buildings in Oslo, 74% of them fulfil this requirement and it has façade imagery available. Coverage is better in the city center where images for 94% of the buildings are available but can be as low as 50% in distant suburbs.

We do not manually check the image quality, such as how many of the requested buildings are included in the image, or whether it is occluded by trees, passing vehicles, scaffolding or similar. Such cases, where the building typology cannot be identified from the image, are labeled as “other”, and the occurrence in the fieldwork area is approximately 27%. An example of this category is shown in **Figure 4**. Labelling of these cases is often difficult and can be subjective, and is expected to be a considerable source of uncertainty for training the automatic labelling procedure.

A total of 5,074 images from fieldwork and from GSV are manually labelled, and they constitute the dataset used for training and validating the CNN model. We set aside 20% (1,019 images) as a test dataset for the final performance evaluation, leaving 4055 images to use for training. Of the latter, we will also set aside a fraction of them for monitoring the training progress that is described in **Section 2**.

For CNN training and prediction, all images are downsampled to 224 × 224 pixels, which is the input resolution of the CNN model. Images from GSV are downloaded in a square format, while rectangular images from fieldwork are center cropped before downsampling. In order to artificially increase the size of the training data, we augment the images by applying the following

<figure><img src="image_93.png" /><figcaption>**FIGURE 6** | Distribution of the predicted building typologies in Oslo, also including the categories “other” and “no images available”.</figcaption></figure>

transformations: randomly zooming in by up to 20%, randomly rotating by up to 25°, and randomly mirroring along the vertical axis. This is done only during training, and it is a standard procedure to improve the model’s ability to generalize.
<figure><img src="image_94.png" /><figcaption>**FIGURE 7** | (A): Google Earth 3D building's view of the example area in Oslo; (B), the predicted MBT for the same example area. The observed and true MBT for the selected area coincides with the predicted one.</figcaption></figure>

# 4 RESULTS

Having completed the fine-tuning, we compute the final performance metrics on the test dataset. The model achieves an accuracy of 0.825, a precision of 0.825, and a recall of 0.825. Without fine-tuning, we obtain an accuracy of 0.763, precision of 0.775 and recall of 0.769, showing that fine-tuning is greatly beneficial to performance. To investigate the classification performance per typology, we present a confusion matrix in **Table 2**. For each true typology (given by the rows), the confusion matrix shows the rate of test images that are assigned to each predicted label (given by the columns). The best performance is seen for timber (T), where 89% of images are classified correctly. This is expected, because timber is the most common typology in the training data, and therefore the category which the CNN has seen the largest variation of. It is also desirable, since we consider the training (and testing) data to be representative of the entire study area, and a high accuracy for the most common category will necessarily lead to a high overall accuracy.

The poorest performance is seen for steel-reinforced concrete (SRC), where 35% of the images are classified correctly, while 60% of them are classified as reinforced concrete (CR). These two typologies can be difficult to distinguish even for experts, as the façade characteristics may not unambiguously determine the typology. The classification of the “other” groups is in general quite successful with 87%. Within the 13% of misclassifications, almost all images in the “other” category are identified as timber buildings, which can be explained by the frequent presence of fences, rooftops and parked cars in these images, which are typical elements in residential areas common for timber buildings. Some examples of images and their predicted labels are shown in **Figure 5**. In **Figure 5** we can observe an example of a timber building correctly classified and a steel-reinforced concrete building misclassified as reinforced concrete.

Concerning the confusion matrix, we note that we did not apply class-specific weights to the images during training, i.e., all

images were given equal priority. It is an option, however, to weight the images by the inverse of the class prevalence, so that rare classes are given higher priority than common ones. This would lead to a more equal classification accuracy across the typologies, instead of having a very high accuracy for the most prominent class (timber) and low accuracy for the least prominent one (steel-reinforced concrete). At the same time, it would also reduce the overall accuracy, and thereby the quality of the final result, which is why we decided against it.

The predicted building typologies are shown in **Table 3**, in terms of number of buildings and percentages, and in **Figure 6** in terms of spatial distribution. It is important to mention that the total number of buildings of Oslo (134,432) also includes buildings where no GSV images are available to perform the classification (“no images available” category includes 35,093 buildings, shown in **Figure 6**). This category is not shown in **Table 3** because we want to present the percentages of the classified buildings that used GSV images for the classification using machine learning. The “no image available” category represents 26.1% of the total number of buildings, and together with the “other” typology, the two categories sum to 45.9% of the total number of buildings. This means that we were able to automatically attribute a model building typology to 54.1% of the buildings in Oslo.

The predominant typology identified in Oslo is timber, that represents 56.7% of the classified buildings. In **Figure 6** we can observe the distribution of the predicted building typologies: most of the classified buildings are localized in the city center and in the urbanized area. GSV images are not available in remote areas (as forest) in Nordmarka and Østmarka: based on our local knowledge and through supervised learning in the area with Google Earth, the main typology identified is timber, represented by private cabins.

**Figure 7** shows the predicted building typologies for an example area in Oslo. This area is chosen because we have a
good representation of different typologies in a small area. On the left side of the figure, the Google Earth 3D building's view is allowing to identify the MBT of the buildings under investigation. The right side of the figure shows the predicted MBT, that coincide with the observed MBT.

# 5 DISCUSSION AND CONCLUSION

This work shows the potential for combining machine learning and publicly available street-level imagery to automate the process of classifying Model Building Typologies for large-scale seismic risk assessment, using Oslo city (Norway) as a case study. Using a state-of-the-art Convolutional Neural Network (CNN) pre-trained on the ImageNet database, we developed a model that classifies typology in unseen images with 83% accuracy, using only data sources available online: the public cadastre and Google Street View. Supplemental high-quality images taken during field work were also used but are not required. Our workflow shows how seismic risk assessment can be highly automated and performed quickly without or limiting time-consuming and costly on-site surveys.

We observe that the classification accuracy varies with building typology: from 89% for timber and down to 35% for steel-reinforced concrete. The reason for this is twofold:

1. first, the distribution of the building mass in Oslo is heavily skewed towards timber buildings, with few steel or steel-reinforced concrete (SRC) buildings to use in the CNN training.

2. Secondly, SRC buildings are typically misclassified as reinforced concrete (CR), which is unfortunately also often the case in manual classification analysis. These misclassifications were verified by breaking down the separate typology predictions for an unseen test dataset by the true labels.

Misclassifications also occur between typologies that share the surrounding environment; in particular, we observe timber buildings surrounded by trees and fences being classified as “other”, and vice versa. This indicates that the CNN to some extent incorporates the environment into the prediction and does not rely on the building properties alone. There exist several techniques (Simonyan et al., 2013; Samek et al., 2017; Selvaraju et al., 2017) to investigate this on a per-image basis, which should be considered in future studies.

The performance of the classifier can be improved by adding more labelled images, which is a matter of additional analyst time. Still, using a pre-trained model reduces the need for training data drastically (Yosinski et al., 2014), as well as lowering the computational cost compared to training a model from scratch. Our method is limited by the availability and good quality of GSV images, which do not offer complete geographical coverage and may be blocked by trees or vehicles. Hence, supplemental field work may be necessary for certain areas. Other online image providers can also complement GSV.

Given the success of this transfer learning approach, also demonstrated by Gonzalez et al. (2020), future work should investigate how well the CNN methodology generalizes and

can hence be applied to other Nordic cities, both with and without additional re-training. We expect that this should be possible at least within Norway, likely also further to Sweden and Finland with generally similar building stocks.

Until recently, the identification of building typologies and the development of a building stock model for seismic risk assessment were limited regarding their spatial coverage as well as financial resources and the lack of representative *in-situ* information. With the data and methodology presented in this paper, these limitations can now be overcome for many areas of the world.

Future work within this topic could include a semi-automatic pre-classification of certain neighborhoods. For instance, we could feed the CNN with pre-conditioning data about the type of the current neighborhood, e.g. residential, commercial, industrial area, which will likely increase the success rate. In addition, one could think about merging information from google-maps roof aspect-ratio with the façade information as done in our current work.

# DATA AVAILABILITY STATEMENT

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/NorwegianSeismicArray/ML-for-building-typologies.

# AUTHOR CONTRIBUTIONS

FG is the main author of the article. She focused in the abstract, introduction, data, results, discussion and conclusion sections. She prepared all the pictures in the article. Her main contribution to the research work is represented by fieldwork in Oslo, taking manually and labelling all the pictures, to develop the different building typologies and to prepare and to analyse the results of the survey sent to the experts. SM, as second author, focused on the methodology, results, discussion and conclusion. He implemented and trained the CNN model, compared the different architectures, and prepared the classification results. AM helped to recognize and identify the different building typologies, to write the description of the typologies and to review the paper. VO helped in the writing and review processes. He discussed with the other authors regarding the best approach to develop the research and the article itself.

# FUNDING

The results and the research showed in this article are connected to the GEObyIT project, funded by the Research Council of Norway (grant number 311596).

# ACKNOWLEDGMENTS

The authors are grateful for the economic support provided by the Research Council of Norway.
REFERENCES

Aravena Pelizari, P., Geiß, C., Aguirre, P., Santa María, H., Merino Peña, Y., and Taubenböck, H. (2021). Automated Building Characterization for Seismic Risk Assessment Using Street-Level Imagery and Deep Learning. *ISPRS J. Photogrammetry Remote Sens.* 180, 370–386. doi:10.1016/j.isprsjprs.2021.07.004

Bungum, H., Pettenati, F., Schweitzer, J., Sirovich, L., and Faleide, J. I. (2009). The 23 October 1904 MS 5.4 Oslofjord Earthquake: Reanalysis Based on Macroseismic and Instrumental Data. *Bull. Seismol. Soc. Am.* 99, 2836–2854. doi:10.1785/0120080357

Chollet, F. (2015). Keras. GitHub. Available at: https://github.com/fchollet/keras.

Chollet, F. (2017). Xception: Deep Learning with Depthwise Separable Convolutions. in Proceedings - 30th IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2017. doi:10.1109/CVPR.2017.195

Danciu, L., Nandan, S., Reyes, C., Basili, R., Weatherill, G., Beauval, C., et al. (2021). The 2020 Update of the European Seismic Hazard Model - ESHM20: Model Overview. *EFEHR Tech. Rep.* 001 v1 0.0, 1–121. doi:10.12686/a15

Eriksson, L., Borchgrevink, J., Sæther, M. M., Daviknes, H. K., Adamou, S., and Andresen, L. (2016). TU1206 COST Sub-urban WG1 Report. Available at: https://static1.squarespace.com/static/542bc753e4b0a87901dd6258/t/5707869ae707eb820b4118a5/1460111042910/TU1206-WG1-012+Oslo+City+Case+Study.pdf.

GFDRR (2018). Machine Learning for Disaster Risk Management. Available at: http://doi.wiley.com/.51, doi:10.1111/j.1468-0394.1988.tb00341.x

Gonzalez, D., Rueda-Plata, D., Acevedo, A. B., Duque, J. C., Ramos-Pollán, R., Betancourt, A., et al. (2020). Automatic Detection of Building Typology Using Deep Learning Methods on Street Level Images. *Build. Environ.* 177, 106805–106812. doi:10.1016/j.buildenv.2020.106805

Grünthal, G. (1998). *European Macroseismic Scale 1998*.

He, K., Zhang, X., Ren, S., and Sun, J. (2016a). Deep Residual Learning for Image Recognition. in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, Las Vegas, NV, USA, 27-30 June 2016 770–778. doi:10.1109/CVPR.2016.90

He, K., Zhang, X., Ren, S., and Sun, J. (2016b). "Identity Mappings in Deep Residual Networks," in *Lecture Notes in Computer Science (Including Subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)*, 630–645. doi:10.1007/978-3-319-46493-0_38

Heaton, J., Bengio, Y., and Courville, A. (2018). *Deep Learning*. MIT Press. Cambridge, Massachusetts, United States, doi:10.1007/s10710-017-9314-z

Huang, G., Liu, Z., Van Der Maaten, L., and Weinberger, K. Q. (2017). Densely Connected Convolutional Networks. in Proceedings - 30th IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2017. doi:10.1109/CVPR.2017.243

Kang, J., Körner, M., Wang, Y., Taubenböck, H., and Zhu, X. X. (2018). Building Instance Classification Using Street View Images. *ISPRS J. Photogrammetry Remote Sens.* 145, 44–59. doi:10.1016/j.isprsjprs.2018.02.006

LeCun, Y., Bengio, Y., and Hinton, G. (2015). Deep Learning. Y. LeCun, Y. Bengio and G. Hinton. *Nature* 521, 436–444. doi:10.1038/nature14539

NVIDIA (2020). NVIDIA V100 TENSOR CORE GPU. Available at: https://images.nvidia.com/content/technologies/volta/pdf/volta-v100-datasheet-update-us-1165301-r5.pdf.

Russakovsky, O., Deng, J., Su, H., Krause, J., Satheesh, S., Ma, S., et al. (2015). ImageNet Large Scale Visual Recognition Challenge. *Int. J. Comput. Vis.* 115, 211–252. doi:10.1007/s11263-015-0816-y

Samek, W., Wiegand, T., and Müller, K.-R. (2017). Explainable Artificial Intelligence: Understanding, Visualizing and Interpreting Deep Learning Models. Available at: http://arxiv.org/abs/1708.08296.

Selvaraju, R. R., Cogswell, M., Das, A., Vedantam, R., Parikh, D., and Batra, D. (2017). Grad-CAM: Visual Explanations from Deep Networks via Gradient-Based Localization. in Proceedings of the IEEE International Conference on Computer Vision, Venice, Italy, 22-29 Oct. 2017 doi:10.1109/ICCV.2017.74

Silva, V., Crowley, H., Pagani, M., Monelli, D., and Pinho, R. (2014). Development of the OpenQuake Engine, the Global Earthquake Model’s Open-Source Software for Seismic Risk Assessment. *Nat. Hazards* 72, 1409–1427. doi:10.1007/s11069-013-0618-x

Simonyan, K., Vedaldi, A., and Zisserman, A. (2013). Deep inside Convolutional Networks: Visualising Image Classification Models and Saliency Maps. Preprint at http://arxiv.org/abs/1312.6034 arXiv.org.

Simonyan, K., and Zisserman, A. (2015). Very Deep Convolutional Networks for Large-Scale Image Recognition. in 3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings.

Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R. (2014). Dropout: A Simple Way to Prevent Neural Networks from Overfitting. *J. Mach. Learn. Res.* 15, 1929–1958.

Szegedy, C., Ioffe, S., Vanhoucke, V., and Alemi, A. A. (2017). Inception-v4, Inception-ResNet and the Impact of Residual Connections on Learning. in 31st AAAI Conference on Artificial Intelligence, AAAI 2017.

Szegedy, C., Vanhoucke, V., Ioffe, S., Shlens, J., and Wojna, Z. (2016). Rethinking the Inception Architecture for Computer Vision. in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition. doi:10.1109/CVPR.2016.308

Wallemacq, P., and House, R. (2018). Economic Losses, Poverty & Disasters: 1998-2017. *Nat. Resour.* Available at: https://www.researchgate.net/profile/Pascaline-Wallemacq-3/publication/331642958_Economic_Losses_Poverty_and_Disasters_1998-2017/links/5c859a6d92851c69506b1f0f/Economic-Losses-Poverty-and-Disasters-1998-2017.pdf?origin=publication_detail

Yosinski, J., Clune, J., Bengio, Y., and Lipson, H. (2014). "How Transferable Are Features in Deep Neural Networks?" in *Advances in Neural Information Processing Systems*.

**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 Ghione, Møland, Meslem and Oye. 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_95.png" /></figure>

Check for
updates

# Evaluating the Minimum Number of Earthquakes in Empirical Site Response Assessment: Input for New Requirements for Microzonation in the Swiss Building Codes

Vincent Perron*, Paolo Bergamo and Donat Fäh

Swiss Seismological Service, Swiss Federal Institute of Technology of Zürich (ETHZ), Zürich, Switzerland

OPEN ACCESS

**Edited by:**
Simone Barani,
University of Genoa, Italy

**Reviewed by:**
Giovanni Forte,
University of Naples Federico II, Italy
Enrico Paolucci,
University of Siena, Italy

**\*Correspondence:**
Vincent Perron
vincent.perron.mail@gmail.com

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

**Received:** 11 April 2022
**Accepted:** 01 June 2022
**Published:** 13 July 2022

**Citation:**
Perron V, Bergamo P and Fäh D (2022) Evaluating the Minimum Number of Earthquakes in Empirical Site Response Assessment: Input for New Requirements for Microzonation in the Swiss Building Codes. Front. Earth Sci. 10:917855. doi: 10.3389/feart.2022.917855

Site-specific hazard analyses and microzonation are important products for densely populated areas and facilities of special risk. The empirical amplification function is classically estimated using the standard spectral ratio (SSR) approach. The SSR simply consists in comparing earthquake recordings on soil sites with the recording of the same earthquake on a close-by rock reference. Recording a statistically significant number of earthquakes to apply the SSR can however be difficult, especially in low seismicity areas and noisy urban environments. On the contrary, computing the SSR from too few earthquakes can lead to an uncertain evaluation of the mean amplification function. Defining the minimum number of earthquake recordings in empirical site response assessment is thus important. We compute empirical amplification functions at 60 KiKnet sites in Japan from several hundred earthquakes and three Swiss sites from several tens of earthquakes. We performed statistical analysis on the amplification functions to estimate the geometric mean and standard deviation and more importantly to determine the distribution law of the amplification factors as a function of the number of recordings. Independent to the site and to the frequency, we find that the log-normal distribution is a very good approximation for the site response. Based on that, we develop a strategy to estimate the minimum number of earthquakes from the confidence interval definition. We find that 10 samples are the best compromise between minimizing the number of recordings and having a good statistical significance of the results. As a general rule, a minimum of 10 uncorrelated earthquakes should be considered, but the higher the number of earthquakes, the lower the uncertainty on the geometric mean of the site amplification function. Moreover, the linear site response is observed to be independent to the intensity of the ground motion level for the analyzed dataset.

**Keywords:** seismic hazard, site effects, microzonation, statistic, signal processing
# INTRODUCTION

Site effects can significantly increase the seismic hazard and risk locally. Unconsolidated deposits such as thick and soft sediments in sedimentary basins are prone to strongly amplify the ground motion. Site effects are caused, among others, by the seismic impedance between rock and sediments, the 1D, 2D and 3D resonances, and the edge-generated surface waves. In turn, the site response can vary significantly from one site to another (site-to-site variability, e.g., Bindi et al., 2009; Hollender et al., 2015; Bindi et al., 2017; Imtiaz et al., 2018; Perron et al., 2018) and from one earthquake to another (within-site variability, e.g., Thompson et al., 2012; Ktenidou et al., 2016; Ktenidou et al., 2017; Maufroy et al., 2017; Perron, 2017; Zhu et al., 2018; Zhu et al., 2022). At large ground motion levels, non-linear effects in specific soils will increase the site response uncertainty as well (Régnier et al., 2013; Régnier et al., 2016). Understanding and reducing the ground motion estimation uncertainty is important for Probabilistic Seismic Hazard Assessment, especially at a long return period (Bommer and Abrahamson, 2006). The site-to-site and within-site variabilities have practical implications for site-specific study and microzonation, for instance on the spatial resolution and required duration of the instrumentation.

The within-site variability is very small when estimated from 1D SH site response analysis because it is a strong simplification of the real phenomena. On the contrary, approaches based on direct observations from real earthquake recordings are appropriate for analyzing the variability of the site response. One of the most commonly used approaches to measure the empirical amplification function is the standard spectral ratio (SSR) introduced by Borcherdt (1970). It consists in performing the ratio in the Fourier domain between the signal recorded at one station on sediments and the signal obtained at another station located nearby on a stiffer site condition (i.e., a rock site) for the same earthquake. However, in noisy urban areas in regions of low-to-moderate seismicity, recording earthquakes with a good signal-to-noise ratio (SNR) can require several months, if not years. It is thus important to estimate the number of earthquakes that should be recorded at the sites to evaluate the empirical site amplification function based on the desired accuracy.

The main goal of this work is to define, for the specification in the Swiss building SIA 261/1 (SIA, 2020), the minimum number of earthquakes in empirical site effect assessment. We first evaluate the stability and validity of the mean amplification as a function of the number of earthquake recordings used to compute it. The variations of the mean amplification are expected to be directly related to the within-site variability at each site. To verify that, we estimate the SSR and surface-to-borehole spectral ratio (SBSR) amplification function for stations of the Swiss strong motion network and of the Japanese KiKnet network having recorded hundreds of earthquakes. We use this large amount of data to determine the statistical distribution of the site amplification. Based on the statistical distribution, we propose an analytic equation predicting the variation of the mean amplification according to the standard deviation and to the number of recorded events. We also determined the dependence

on the mean amplification functions of the ground motion intensity, measured as the peak ground acceleration (PGA).

# METHOD AND RESOURCES

In Switzerland, we developed a waveform database covering the time period from January 1998 to September 2019. Waveforms at each Swiss site were selected according to a magnitude-distance filter. In Japan, the database is covering the time period from October 1997 to March 2016. The SSR is computed for each component individually or the mean of the two horizontal components and can be noted as follows:

$$SSR_i(f) = \frac{FAS_{Si}(f)}{FAS_{Ri}(f)}, \quad (1)$$

where SSR_i(f) is the SSR for the ith component as a function of frequency f and FAS_Si and FAS_Ri are respectively the Fourier amplitude spectra (FAS) at the site and at the reference computed over the ith component. The SSR approach is based on the assumption that the earthquake source and wave propagation along the path are the same between the site and the reference and thus canceled out when performing the spectral ratio between the two. This assumption is valid if the site-to-reference distance (R_STA) is much smaller than the hypocentral distance (Rh). In practice, adopting Rh > 10R_STA is considered to be enough, even though a certain part of the SSR variability can probably be explained by a remaining influence of the source and of the path (Borcherdt, 1970; Perron, 2017). The ground motion amplification at the reference station is assumed to be negligible, that is to say, equal to one at every frequency. In practice, it is never the case (Hollender et al., 2017; Hollender et al., 2018; Hobiger et al., 2021), so the SSR-based amplification factors are not absolute but are always relative to the considered reference. One of the main limitations of the SSR is of having a rock outcropping susceptible to be used for the reference site located not too far from the considered site of interest. An alternative to the classical SSR is to deploy one station at the earth's surface on sediments and the second at the same location but in a borehole deep enough to reach the geophysical bedrock. This so-called SBSR approach has the advantage of solving the between-station distance limitation but introduces some new difficulties because of the seismic wave reflection at the earth's surface. The upgoing and downgoing waves are indeed fully constructive at the earth's surface, although they can be destructive at certain frequencies at depth (Cadet et al., 2012). However, in the context of analyzing only the variability of the site response, the downgoing wave interaction can reasonably be neglected (Cadet et al., 2012; Hollender et al., 2018). We followed the same procedure for every computation of the site response in Switzerland and Japan. This procedure is as follows:

1) Automatic quality checks of earthquake recordings and automatic picking of the P and S wave arrival (T_P, T_S) through a time-frequency analysis;
A

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

B

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

C

D

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

E

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

G

**FIGURE 1** | Example of the surface-to-borehole spectral ratio (SBSR) computation at KiKnet station IBRH12 in Japan. (A) The map shows the location of the site (green triangle) and the epicenters of the selected earthquakes (yellow-to-red dot according to the earthquake magnitude). Panels (B) and (C) present the power spectral density (PSD) for the noise (black lines) and for the earthquake recording on the horizontal mean component at the site (blue lines) and at the reference (green lines). Panel (D) indicates the SNR at the site (blue lines) and at the reference (green lines), as well as the number of earthquakes spectrum with SNR > 5 (red line) as a function of frequency. The distribution of the SBSR as a function of frequency for the horizontal mean component (E), for the horizontal as a function of the azimuth (F), and for the vertical (G) component. The color scale indicates the density of lines, each line corresponding to the SBSR of one single earthquake.

2) Selection of earthquakes with hypocentral distance at least five times the interstation distance ($R_{STA}$);

3) Selection of the signal window between $T_p$ and the coda defined by $3.3T_s-2.3T_p$ (Perron et al., 2017) and of the noise window before $T_p$ and of the same duration as the signal window. Site and reference use the same time windows;

4) Computation of the FAS for the noise and the signal window;

5) Computation of the horizontal mean FAS using the quadratic mean: $\sqrt{\frac{N^2+E^2}{2}}$;

6) Smoothing and resampling of the horizontal mean FAS on a logarithmic scale using the Konno and Ohmachi (1998) approach with a b-value of 50;

7) Estimation of the SNR;

8) Selection of earthquakes with SNR > 5 over at least a two-octave frequency band window both at the site and at the reference;

9) Spectral ratio computation between the horizontal mean FAS at the site and at the reference for each earthquake;

10) Estimation of the within-site events geometric mean and standard deviation at each frequency;

11) Detection of outliers as a group of samples of probability <0.1% over a frequency band larger than one octave;

12) Outliers are discarded, and the geometric mean and standard deviation are recomputed

Figure 1 shows an example of the SBSR computation in Japan.

## STANDARD SPECTRAL RATIO AND SURFACE-TO-BOREHOLE SPECTRAL RATIO RESULTS

In total, SSR is estimated from three pairs of stations where approximately 100 good-quality earthquakes have been recorded in Switzerland, and SBSR is computed from 60 pairs of surface-to-borehole stations with up to 2000 good-quality earthquakes in Japan. Figure 2 and Figure 3 respectively show the distribution of the SBSR for 60 pairs of surface-to-borehole stations in Japan and the SSR for the three pairs of surface stations in Switzerland. Figure 4 provides a summary of the number of good-quality earthquake recordings, geometric mean, and geometric standard deviation as a function of frequency in Japan (gray curves) and Switzerland (red curves).

Figure 2, Figure 3, and Figure 4 clearly show that the amplification functions are different from one site to another, both in terms of mean and standard deviation. It reflects the differences in the geological conditions of the sites, which determine, among others, the fundamental resonance
<figure><img src="image_100.png" /><figcaption>**FIGURE 2** | Amplification function computed from the SBSR between 60 pairs of stations in Japan. The color from dark blue to light green indicates an increasing density of curves, each curve corresponding to one single earthquake.</figcaption></figure>
<figure><img src="image_101.png" /></figure>

**FIGURE 3** | Amplification function computed from the standard spectral ratio between 3 pairs of stations in Switzerland. The color from dark blue to light blue indicates an increasing density of curves, each curve corresponding to one single earthquake.

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

**FIGURE 4** | Number of good-quality earthquakes (left panel), within-site geometric mean (central panel), and within-site geometric standard deviation (right panel) as a function of frequency for 60 surface-to-borehole spectral ratios in Japan (gray curves) and three standard spectral ratio in Switzerland (red curves).

frequency of the site ($f_0$), here corresponding to the first peak on the amplification function. The within-site standard deviation can also vary drastically from one site to another and depending on the frequency. In Japan, we can separate the amplification functions into two groups: the first group with $f_0 > 0.5$ Hz, with an amplification function equal to one and low standard deviation (close to 1.05) for frequency below $f_0$; a second group with $f_0$ below the minimum frequency of the analysis here (0.1 Hz), and having significant amplification (above one) and high variability at low frequency. It is also clear that the variability of the site response is on average higher in Switzerland than that in Japan. For the Swiss sites, this is probably because of the SSR method imposing relatively high site-to-reference distances and non-negligible site effects at the surface reference station. In Japan, we can observe some anomalies (eye shapes departing from the log-normal distribution) in the amplification function at high frequency (e.g., for stations: KiK-IBRH13; KiK-IBRH17; KiK-TCGH16). It is not possible to clearly determine its origin, but from our experience, this is very probably an artificial artifact because of coupling issues of the borehole instrumentation or

because of a modification on the instrumentation at some point due to maintenance of the station for instance.

## DISTRIBUTION OF THE WITHIN-SITE VARIABILITY

As we have seen in the previous section, both the mean and standard deviation of the amplification function as a function of frequency are dependent to the geological characteristics of the site itself. However, the nature of the site response distribution is the same independently to the site or to the frequency and has been shown to be well modeled by a log-normal distribution (Ktenidou et al., 2011). In other words, the distribution of the logarithm of the relative amplification of the ground motion between two sites is Gaussian. To qualitatively verify the log-normal distribution of the site response at every frequency, the quantile-quantile (Q-Q) plot and the histogram are represented at frequencies 0.5, 1.0, 2.5, 5.1, 9.9, and 20.6 Hz in **Figure 5**. The shape of the histograms of the logarithm of the amplification
<figure><img src="image_103.png" /><figcaption>**FIGURE 5** | Quantile-quantile plot of the logarithm of the amplification factors for 60 surface-to-borehole spectral ratios in Japan (gray curves) and three standard spectral ratios in Switzerland (red curves) at six different frequencies (one panel per frequency). On each panel, the histogram (gray area) of the standard normal distribution computed from the logarithm of the amplification factors at all sites at the corresponding frequency is compared with the best normal distribution fit (green curve).</figcaption></figure>

factors represents a Gaussian and Q-Q curves of every site at every frequency are well aligned along the 1/1 line, in particular in the interval ± 2σ to the mean. These indicate that the site response is very well approximated by log-normal distribution at least up to ± 2σ. Beyond 2σ, the few non-natural outliers and the limited number of samples increase the scatter of the Q-Q curves, meaning that the log-normal distribution is still valid but interpretations made out of it are less reliable.

Proving the log-normal distribution of the amplification function is important because then peculiar statistical properties apply. For example, if a variable x is normally distributed then the distribution of sample means (x̄_n) computed from subsets of n samples also are normally distributed. One major output of that is the confidence interval (Ic). Given that a sample mean (x̄_n) and unbiased standard deviation (s_n) have been estimated from a finite number of samples (n), the confidence interval is the interval inside which the population mean (μ) for an infinite number of samples has a certain confidence level to be included in. It is defined as follows:

$$Ic_{1-\alpha\%} = \left[ \bar{x}_n - Z_{\alpha/2} \frac{s_n}{\sqrt{n}}; \bar{x}_n + Z_{\alpha/2} \frac{s_n}{\sqrt{n}} \right], \quad (2)$$

where Z_α/2 is the critical value that defines the confidence level (1 - α). For a normal distribution and a confidence level of 95%, Z_0.025 is equal to 1.96. However, because the number of samples can be sometimes very limited (i.e., only a few earthquakes have been recorded), it is preferable to use the Student distribution,

also called t-distribution. This distribution correctly accounts for a small number of samples and tends to be a normal distribution as the number of samples increases. For a Student distribution, the formulation of Ic_1-α% is the same (Eq. 2), but the estimation of Z_α/2 is different, as it now also depends on n. The evolution of Z_α/2,n as a function of n and for the confidence levels 68, 95, 99, and 99.9% is given in Figure 6, left panel. In the following, we will keep using the notation x̄_n and s_n for the measured sample geometric mean and standard deviation, whereas μ and σ represent the population geometric mean and standard deviation of the distribution. For an infinite number of samples, the two notations become equivalent: x̄_∞ = μ and s_n = σ. Moreover, we will only focus on the confidence level of 95%, because the 95% confidence interval corresponds approximately to the interval comprised between [-1.96σ, 1.96σ], which in turn corresponds to the portion where the Q-Q plot best fit the 1/1 line (Figure 5). As the distribution is not normal but log-normal, we accordingly modified the confidence interval formulation. The 95% confidence interval for a log-Student distribution is finally:

$$Ic_{95\%} = \left[ \bar{x}_n * \frac{1}{\exp\left(Z_{0.025, n} \frac{\ln(s_n)}{\sqrt{n}}\right)}; \bar{x}_n * \exp\left(Z_{0.025, n} \frac{\ln(s_n)}{\sqrt{n}}\right) \right], \quad (3)$$

with x̄_n and s_n respectively the sample geometric mean and standard deviation computed as
<figure><img src="image_104.png" /><figcaption>FIGURE 6 | Critical value Z (left panel) and confidence interval (right panel) as a function of the number of samples for the confidence levels 68%, 95%, 99%, and 99.9% for the standard normal distribution (dashed lines), and the standard Student distribution (solid lines).</figcaption></figure>

$$ \bar{x}_n = \exp\left(\frac{1}{n}\sum_{i=1}^{n} \ln(x_i)\right), \quad (4) $$

$$ s_n = \exp\left(\sqrt{\frac{1}{n-1}\sum_{i=1}^{n}(\ln(x_i) - \ln(\bar{x}_n))^2}\right), \quad (5) $$

Figure 6 (right panel) shows the evolution of $Ic_{68\%}$, $Ic_{95\%}$, $Ic_{99\%}$, and $Ic_{99.9\%}$ for a standard normal and standard Student distribution ($\mu = 0$; $\sigma = 1$). It illustrates the very rapid reduction of the confidence interval as the number of samples increases, from more than 10 $\sigma$ when $n < 10$ to less than 1 $\sigma$ when $n > 10$.

## VALIDITY OF THE CONFIDENCE INTERVAL PREDICTIONS

After demonstrating the validity of the log-normal assumption, we verified the validity of the prediction of $Ic_{95\%}$ for a Student distribution as a function of the number of earthquakes $n$ by comparing $Ic_{95\%}$ with the observations in Switzerland and Japan. First, we defined two different confidence intervals:

$$ Ic_{95N}(n) = \left[ \bar{x}_N * \frac{1}{\exp\left(Z_{0.025, N} \frac{\ln(s_N)}{\sqrt{n}}\right)} ; \bar{x}_N * \exp\left(Z_{0.025, N} \frac{\ln(s_N)}{\sqrt{n}}\right) \right] \quad (6) $$

$$ Ic_{95n}(n) = \left[ \bar{x}_n * \frac{1}{\exp\left(Z_{0.025, n} \frac{\ln(s_n)}{\sqrt{n}}\right)} ; \bar{x}_n * \exp\left(Z_{0.025, n} \frac{\ln(s_n)}{\sqrt{n}}\right) \right] \quad (7) $$

where $\bar{x}_N$ and $s_N$ are respectively the total geometric mean and standard deviation computed over the entire dataset of N events. $\bar{x}_n$ and $s_n$ are respectively the local geometric mean and standard

deviation computed over a subset of $n$ randomly selected events. $Ic_{95N}(n)$ is the total 95% confidence interval used to predict the variation of any local mean $\bar{x}_n$ computed from $n$ events. $Ic_{95n}(n)$ is the local 95% confidence interval used to predict the interval of variation of the total mean $\bar{x}_N$. This assumes that $\bar{x}_N = \mu$ and $s_N = \sigma$, which is reasonably correct here since $N$ is most of the time much higher than 100 earthquakes.

To estimate the reliability of the confidence interval more quantitatively, we bootstrapped the amplification factors at each frequency over 1000 random selections of $n$ events, with $n = [2 3 4 6 8 10 14 18 24 32]$. We evaluated the proportion of local means included inside the total confidence interval ($P_1 = \bar{x}_n \subset Ic_{95N}(n)$), and the proportion of total means included inside the local confidence interval ($P_2 = \bar{x}_N \subset Ic_{95n}(n)$). Following Eq. 3, $P_1$ and $P_2$ can be written:

$$ P_1(f, n) = \frac{1}{1000} \sum_{k=1}^{1000} \left( \bar{x}_N(f) * \frac{1}{\exp\left(Z_{0.025, N} \frac{\ln(s_N(f))}{\sqrt{n}}\right)} \le \bar{x}_{nk}(f) \right. \\ \left. \le \bar{x}_N(f) * \exp\left(Z_{0.025, N} \frac{\ln(s_N(f))}{\sqrt{n}}\right) \right) \quad (8) $$

$$ P_2(f, n) = \frac{1}{1000} \sum_{k=1}^{1000} \left( \bar{x}_{nk}(f) * \frac{1}{\exp\left(Z_{0.025, n} \frac{\ln(s_{nk}(f))}{\sqrt{n}}\right)} \right. \\ \left. \le \bar{x}_N(f) \le \bar{x}_{nk}(f) * \exp\left(Z_{0.025, n} \frac{\ln(s_{nk}(f))}{\sqrt{n}}\right) \right) \quad (9) $$

with the inequation equal to 1 when it is true and 0 otherwise. $\bar{x}_{nk}$ and $s_{nk}$ are respectively the $k^{th}$ local geometric mean and standard
<figure><img src="image_105.png" /></figure>

**FIGURE 7** | Evaluation of P1 (dark brown line) and P2 (light brown line) on the standard spectral ratio computed at Swiss station SIOO/SIOV from 1000 randomly selected subsets of *n* = {2 3 4 6 8 10 14 18 24 32} earthquakes (top-left to bottom-right panel). On each panel, the left axis provides the amplification scale and the right axis indicates P1 and P2 proportion in percentages. The 1000 local means $\bar{x}_n$ are represented according to their density of points from dark blue to light green. The total 95% confidence interval prediction for *n* events ($IC_{95\%}$) is represented with orange lines. P1 can be easily visualized by looking at the proportion of $\bar{x}_n$ points exceeding the $IC_{95\%}$ (blue points outside the orange lines). There is no way to simply represent P2 here. The number of events, and the mean P1 and P2 over the frequency range are written on each panel.

deviation computed over a subset of *n* randomly selected events. If the distribution is perfectly normal, then both $P_1$ and $P_2$ are equal to 95%. However, we do not expect the site response distribution to be perfectly normal at every site and for all frequencies, so a certain convergence to 95% should be observed as the number of events *n* increases.

**Figure 7** shows the bootstrap estimation of $P_1$ and $P_2$ from the amplification function of SIOO/SIOV in Switzerland. First, it is clear that the variability between the 1000 $\bar{x}_{nk}$ decreases (blue points) as *n* increases (from top-left to bottom-right panel). This decay seems well predicted by $IC_{95\%}$ (orange lines). This observation is also supported by $P_1$ which is relatively close to the value of 95% at all frequency and for any $n > 2$. For $n = 2$ we can observe that $P_1$ is slightly higher than 95% between 10 and 20 Hz. In contrast, $P_2$ shows some significant low values for any

$n < 10$. However, $P_2$ shows a better agreement with the 95% value as *n* increases. This observation confirms the good approximation of using the log-normal distribution to model the site amplification variability. $IC_{95\%}$ makes a relatively good prediction of the observed variability of $\bar{x}_n$, even when the number of samples is low.

Now, we follow the same procedure for every three SSR in Switzerland and SBSR in Japan. The corresponding results are given in **Figure 8**. We can make a similar observation as in **Figure 7**, $P_1$ is the average equal to the 95% value at all frequency and for every number of events. For the Swiss SSR, we can, however, observe a stronger scatter when the numbers of events are minimum ($n < 4$). Again, we observe a stronger deviation of 95% in $P_2$ both in Switzerland and in Japan. In Switzerland, the discrepancy of $P_2$ is higher, especially close to 1 Hz and for $n < 6$.
P₁: proportion of $\bar{x}_n \subset Ic_{95N}$ (%)

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

P₂: proportion of $\bar{x}_N \subset Ic_{95n}$ (%)

<figure><img src="image_107.png" /><figcaption>**FIGURE 8** | Evaluation of P1 and P2 as a function of the number of events *n* at 6 frequencies (0.5, 1.0, 2.5, 5.1, 9.9, and 20.6 Hz) for three Swiss sites and 60 Japanese sites.</figcaption></figure>

$P_2$ is an average lower than 95% but tends to it as *n* increases. A good agreement is found for *n* > 6 and a complete stabilization is observed above 14 events. In Japan, we observed a different behavior, with $P_2$ being too low when *n* = 2, and then too high when $2 < n < 8$ mainly at low frequency (f < 2 Hz). For *n* > 8, we observed a good stabilization of $P_2$ with mean values slightly below 95%.

The confidence interval computed from a large site response dataset is a good estimator of what is going to be the behavior of the mean computed from much smaller subsets of even only three earthquakes and for any frequency. However, it is clear that using 10 recordings of earthquakes or above greatly improves the quality of the prediction and the significance of the results. In conclusion, at least 10 events should be considered to have a good
statistical significance and to make good use of the confidence interval predicting power.

## VARIABILITY OF THE MEAN AMPLIFICATION FUNCTION AS A FUNCTION OF THE NUMBER OF EVENTS

Some questions which arise when evaluating the amplification function at a specific site are as follows: Is the number of earthquake recordings sufficient to accurately estimate the amplification function? Which minimal number of earthquakes ($n_{min}$) should be used to evaluate the site response? Based on the confidence interval definition (Eq. 3), it is clear that the variability of $\bar{x}_n$ depends both on $s_n$ and $n$. Because $s_n$ is site- and frequency-dependent (Figure 4), $n_{min}$ is by consequence also site- and frequency-dependent. In other words, there is no unique value of $n_{min}$ which can be considered for every site response analysis in the world. On the other hand, the property of the site response to be log-normally distributed can be supposed as universal. It is then possible to determine $n_{min}$ for any site response analysis, based on the log-normal distribution assumption and the use of the confidence interval definition.

Provided that the geometric mean $\bar{x}_n$ and standard deviation $s_n$ of the site response has been measured at a particular site over a certain number of earthquakes $n$, it is possible to determine in which confidence interval the population mean for an infinite number of events $\mu$ has a certain confidence level (here 95%) to be included in. It is also possible to predict what will be the reduction of this interval if the number of earthquake observations increases. In the same way, it is possible to determine the number of earthquakes required to limit to a certain level the width of the interval where $\mu$ has a 95% probability to be found within. The width of the interval is independent to the $\bar{x}_n$ and can be defined from Eq. 3 by

$$C_{95\%} = \exp\left(Z_{0.025, n} \frac{\ln(s_n)}{\sqrt{n}}\right), \quad (10)$$

$C_{95\%}$ is the coefficient of variation between $\mu$ and $\bar{x}_n$ such as $\frac{\bar{x}_n}{C_{95\%}} \le \mu \le C_{95\%}\bar{x}_n$ with a 95% probability. It is now possible to estimate the minimum number of earthquakes required to limit the variation between $\mu$ and $\bar{x}_n$ below a certain coefficient $C_{95\%}$ as

$$n_{min} = \left(Z_{0.025, n} \frac{\ln(s_n)}{\ln(C_{95\%})}\right)^2 \quad (11)$$

For example, if the amplification at 1 Hz has been measured from $n = 10$ earthquakes with a geometric standard deviation of $s_{10} = 1.5$, we can estimate the minimum number of earthquake $n_{min}$ to have $C_{95\%} = 1.2$ (20% of variation) with a probability of 95% as

$$n_{min} = \left(Z_{0.025,10} \frac{\ln(s_{10})}{\ln(C_{95\%})}\right)^2 = \left(2.26 \frac{\ln(1.50)}{\ln(1.20)}\right)^2 \\ = 25.31 \rightarrow 26 \text{ earthquakes}$$

It is important to note that for a Student distribution, $Z_{0.025, n}$ is the function of $n$. $Z_{0.025, n}$ will decrease very

rapidly as the number of measured earthquakes increases (Figure 6). Using Eq. 11 and measured $s_n$ (Figure 4), $n_{min}$ is computed for every site in Switzerland and Japan, and at every frequency. The results are reported in Figure 9. As already discussed, $n_{min}$ is dependent on $s_n$, so it is variable for the different sites and frequency. Swiss SSRs have the highest uncertainty and logically required the highest number of earthquakes for a given coefficient of variation $C_{95\%}$. Table 1 summarizes the minimum number of earthquakes which is valid for 99, 95, and 84% of our sites and frequencies as a function of $C_{95\%}$. For 10 earthquakes recorded, the estimation of the mean is only 40% accurate approximately ($C_{95\%} = 1.4$). It is possible to reduce this uncertainty to 25% by recording 20 events ($C_{95\%} = 1.25$). Depending on the desired limit for the coefficient of variation of the mean, one can make own estimations of the minimum number of earthquakes using Eq. 7.

It has to be highlighted that $s_n$ is the key parameter for the estimation of $n_{min}$. If $s_n$ is wrongly determined, so will be $n_{min}$. One difficulty to have a representative determination of $s_n$ is how to deal with the outliers. Including erratic outliers will artificially increase $s_n$, while removing natural outliers from rare events will truncate the true distribution and reduce $s_n$. Another difficulty is that looking only at the value of $n_{min}$ might not be enough for all sites. One could claim that because the site response has been measured from 30 earthquakes, the statistical significance of the result is good and the coefficient of variation of the mean is low. However, if all the events present the same characteristic and location because they belong to the same cluster of events, then the significance of the results is not good and the true variability of the site response might be strongly underestimated. For instance, Perron (2017) showed that approximately 50% of the within-site variability in 2D and 3D basins comes from the lighting effect, which strongly depends on the source location. This implies that both the number of events and their spatial distribution around the site should be considered in site response analysis.

## DEPENDENCE OF THE SITE RESPONSE VARIABILITY ON THE INTENSITY OF THE GROUND MOTION

The dependence of the site response on the intensity of the ground motion is a complex research topic that interests the community for several decades (e.g., Sánchez-sesma, 1987; Aki, 1993). The non-linear behavior of unconsolidated soil to strong ground motion solicitations is of major interest in engineering seismology. Non-linearity tends to reduce the fundamental resonance frequency of the site, leading to an increase of the hazard at low frequency and a decrease at high frequency (Régnier et al., 2016). In extreme cases, it can also lead to liquefaction phenomena.

One question often arises when speaking about empirical site effects assessment which is: is the measured amplification function from weak ground motion representative of site
<figure><img src="image_108.png" /><figcaption>FIGURE 9 | Minimum number of earthquakes as a function of frequency for the coefficient of variation C<sub>95%</sub> equal to 1.05, 1.10, 1.15, 1.20, 1.25, 1.30, 1.4, and 1.5 (panels). KiKnet stations with f<sub>0</sub> > 0.5 Hz are represented in black, KiKnet stations with f<sub>0</sub> < 0.1 Hz are represented in gray, and Swiss stations are represented in red.</figcaption></figure>

<table><caption>TABLE 1 | Minimum number of earthquakes *n*<sub>min</sub> as a function of the coefficient of variation C<sub>95%</sub>.</caption><thead><tr><th>C<sub>95%</sub></th><th>1.05<br>(5%)</th><th>1.10<br>(10%)</th><th>1.15<br>(15%)</th><th>1.20<br>(20%)</th><th>1.25<br>(25%)</th><th>1.30<br>(30%)</th><th>1.40<br>(40%)</th><th>1.50<br>(50%)</th></tr></thead><tbody><tr><th>n<sub>min 99%</sub></th><td>403</td><td>106</td><td>50</td><td>29</td><td>20</td><td>14</td><td>9</td><td>6</td></tr><tr><th>n<sub>min 95%</sub></th><td>214</td><td>56</td><td>26</td><td>16</td><td>11</td><td>8</td><td>5</td><td>4</td></tr><tr><th>n<sub>min 84%</sub></th><td>109</td><td>29</td><td>14</td><td>8</td><td>6</td><td>4</td><td>3</td><td>2</td></tr></tbody></table>

response to strong ground motion? To address this question, we compute the equivalent of the standard normal distribution (μ = 0, σ = 1) for every individual amplification function at all sites in Switzerland and Japan as

$$Z_i = \frac{\ln(xi) - \ln(\bar{x}_n)}{\ln(s_n)}, \qquad (12)$$

This common standard normal distribution formulation allows using the site response of every site together. *Zi*(*f*) represent the ith normalized amplification function normally distributed with x̄<sub>n</sub> = 0 and s<sub>n</sub> = 1. Together, it represents about 28,000 normalized amplification functions obtained from thousands of earthquakes recorded at 63 pairs of stations (three Swiss sites and 60 Japanese sites). For each normalized

amplification function, we computed on the corresponding waveforms the horizontal mean PGA.

**Figure 10** shows the number of events per frequency, the distribution of the PGA and the normalized amplification function for four PGA bins [(0.001 0.01), (0.01 0.1), (0.1 1), and (1 10) m/s²]. First, it should be mentioned that the number of events varies strongly from one PGA bin to another. This explains the apparent differences when looking at the normalized amplification function (black curves) of the different bin. We observe that the normalized amplification function for every PGA bin can be explained by the standard normal distribution, which indicates that no non-linear behavior is observed here. The mean is fairly equal to 0 and the standard deviation is equal to 1 for every frequency of every bin. That demonstrates, first, that the linear behavior characterizes the vast majority of the sites, and second, that the linear site response is independent to the ground motion intensity. Therefore, if we consider a specific site having a linear behavior, the amplification function observed from the weak motion of a small magnitude earthquake will be the same as the one for the strong motion of a large magnitude earthquake, all other things being the same. This highlights the importance and the validity of using the recording of low-to-moderate earthquakes to assess the anelastic amplification functions for larger earthquakes as long as there is no significant non-linear site response at the site of interest.
<figure><img src="image_109.png" /><figcaption>**FIGURE 10** | Top-left panel: Total number of normalized amplification functions obtain from 3 Swiss SSR distribution and 60 SBSR Japanese distribution and as a function of frequency. Top right: Histogram of the peak ground acceleration (PGA) distribution. From middle left to bottom-right panel: normalized amplification function for four PGA bins and as a function of frequency. The mean and mean plus/minus standard deviation are represented with solid red lines and dotted red lines respectively.</figcaption></figure>

CONCLUSION

Site effect is a major contributor to the seismic hazard, and its
evaluation at specific sites of interest generally requires the
recording of several earthquakes. We address here the
question of the site response variability and of the minimum
necessary number of earthquakes to be recorded.

To address this question, we carefully compute empirical
amplification functions at 60 KiKnet sites from several
hundred earthquakes and three Swiss sites from several tens of
earthquakes. We performed statistical analysis on the
amplification function to estimate the geometric mean and
standard deviation, and more importantly to determine the
distribution law of the amplification factor at each frequency.
Independent to the site and to the frequency, we found that the log-normal distribution is a very good approximation for the site response. Based on that we developed a strategy to estimate the minimum number of earthquakes from the confidence interval definition. We first demonstrate the validity of the use of the confidence interval to model the uncertainty of the geometric mean estimation. We found that between 8 and 14 earthquakes are necessary to have a good prediction by the confidence interval, that is to say, a good statistical significance. For most of the sites, 10 samples seem to be the best compromise between minimizing the number of recordings and having a good statistical significance of the results. Based on the confidence interval, we provide the analytic formula to estimate the minimum number of earthquakes to be recorded, as a function of the within-site standard deviation (Eq. 11). We used it on the Swiss and Japanese amplification function and determine, among others, that with a 95% probability: the mean varies by less than 40% for 10 earthquakes, and less than 25% for 20 events.

It is very important to point out that satisfying the minimal number of earthquakes by itself is not sufficient. The selected earthquakes should be uncorrelated and as much evenly distributed around the site as possible to cover the entire variability of the site response. Therefore, one should not use only earthquakes belonging to a single cluster of events. In our dataset, the linear site response is observed to be independent to the intensity of the ground motion. In other words, assessing the site response from the recording of low PGA and low magnitude earthquakes, provides the same amplification functions as from recording of high PGA and large magnitude earthquakes, as far as the soil behaves linearly.

As a general rule, a minimum of 10 uncorrelated earthquakes should be considered, but the higher the number of earthquakes, the lower the uncertainty on the geometric mean site response assessment. Based on our results, the specification in the Swiss building SIA 261/1 recommends taking a minimum of 10 uncorrelated earthquakes to perform site-specific studies.

## REFERENCES

* Aki, K. (1993). Local Site Effects on Weak and Strong Ground Motion. *Tectonophysics* 218, 93-111. doi:10.1016/0040-1951(93)90262-i

* Bindi, D., Luzi, L., and Pacor, F. (2009). Interevent and Interstation Variability Computed for the Italian Accelerometric Archive (ITACA). *Bull. Seismol. Soc. Am.* 99, 2471-2488. doi:10.1785/0120080209

* Bindi, D., Spallarossa, D., and Pacor, F. (2017). Between-Event and Between-Station Variability Observed in the Fourier and Response Spectra Domains: Comparison with Seismological Models. *Geophys. J. Int.* 210, 1092-1104. doi:10.1093/gji/ggx217

* Bommer, J. J., and Abrahamson, N. A. (2006). Why Do Modern Probabilistic Seismic-Hazard Analyses Often Lead to Increased Hazard Estimates? *Bull. Seismol. Soc. Am.* 96, 1967-1977. doi:10.1785/0120060043

* Borcherdt, R. D. (1970). Effects of Local Geology on Ground Motion Near San Francisco Bay. *Bull. Seismol. Soc. Am.* 60, 29-61.

* Cadet, H., Pierre-Yves, B., and Adrian, R.-M. (2012). Site Effect Assessment Using KiK-Net Data: Part 1. A Simple Correction Procedure for Surface/Downhole Spectral Ratios. *Bull. Earthq. Eng.* 10, 421-448. doi:10.1007/s10518-011-9283-1

## DATA AVAILABILITY STATEMENT

Publicly available datasets were analyzed in this study. This data can be found here: The data of the Swiss stations (CH, https://doi.org/10.12686/sed/networks/ch) can be accessed following the instruction on the webpage http://www.seismo.ethz.ch/en/research-and-teaching/products-software/waveform-data/ (last accessed April 2022). The data of the station site characterization can be accessed at http://stations.seismo.ethz.ch (last accessed April 2022). The figures were produced using MATLAB which is available at www.mathworks.com/products/matlab (last accessed August 2021). Japaness KiKnet data are available at https://www.kyoshin.bosai.go.jp/ (last accessed June 2022).

## AUTHOR CONTRIBUTIONS

VP did most of the work, as well as the writing of the article. PB helped with the statistic and reviewing work. DF is the main supervisor and fund provider.

## FUNDING

This work received financial contributions from the Swiss Federal Office for the Environment (FOEN), the Swiss Federal Office for Civil Protection (FOCP), and the Swiss Federal Institute of Technology Zurich (ETHZ). Open access funding is provided by ETH Zurich.

## ACKNOWLEDGMENTS

This work was made in the framework of the Earthquake Risk Model for Switzerland project financed by contributions from the Swiss Federal Office for the Environment (FOEN), Swiss Federal Office for Civil Protection (FOCP), and Swiss Federal Institute of Technology Zurich (ETHZ).

* Hobiger, M., Bergamo, P., Imperatori, W., Panzera, F., Marrios Lontsi, A., Perron, V., et al. (2021). Site Characterization of Swiss Strong-Motion Stations: The Benefit of Advanced Processing Algorithms. *Bull. Seismol. Soc. Am.* 111, 1713-1739. doi:10.1785/0120200316

* Hollender, F., Cornou, C., Dechamp, A., Oghalaei, K., Renalier, F., Maufroy, E., et al. (2017). Characterization of Site Conditions (Soil Class, VS30, Velocity Profiles) for 33 Stations from the French Permanent Accelerometric Network (RAP) Using Surface-Wave Methods. *Bull. Earthq. Eng.* 16, 2337-2365. doi:10.1007/s10518-017-0135-5

* Hollender, F., Perron, V., Imtiaz, A., Svay, A., Mariscal, A., Bard, P.-Y., et al. (2015). "Close to the Lair of Odysseus Cyclops: The SINAPS@ Post-seismic Campaign and Accelerometric Network Installation on Kefalonia Island," in 9ème Colloque National AFPS, Marne-la-Vallée, France, 30 Nov-2 Dec 2015.

* Hollender, F., Roumelioti, Z., Regnier, J., Perron, V., and Bard, P.-Y. (2018). "Respective Advantages of Surface and Downhole Reference Stations for Site Effect Studies: Lessons Learnt from the Argonet (Cephalonia Island, Greece) and Cadarache (Provence, France) Vertical Arrays," in 16th European Conference on Earthquake Engineering (16ECEE), Thessaloniki, Greece, 18-21June 2018, 13.
Imtiaz, A., Perron, V., Hollender, F., Bard, P. Y., Cornou, C., Svay, A., et al. (2018). Wavefield Characteristics and Spatial Incoherency: A Comparative Study from Argostoli Rock- and Soil-Site Dense Seismic Arrays. *Bull. Seismol. Soc. Am.* 108, 2839–2853. doi:10.1785/0120180025

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

Ktenidou, O.-J., Chávez-García, F.-J., Raptakis, D., and Pitilakis, K. D. (2016). Directional Dependence of Site Effects Observed Near a Basin Edge at Aegion, Greece. *Bull. Earthq. Eng.* 14, 623–645. doi:10.1007/s10518-015-9843-x

Ktenidou, O.-J., Chavez-Garcia, F. J., and Pitilakis, K. D. (2011). Variance Reduction and Signal-To-Noise Ratio: Reducing Uncertainty in Spectral Ratios. *Bull. Seismol. Soc. Am.* 101, 619–634. doi:10.1785/0120100036

Ktenidou, O.-J., Roumelioti, Z., Abrahamson, N., Cotton, F., Pitilakis, K., and Hollender, F. (2017). Understanding Single-Station Ground Motion Variability and Uncertainty (Sigma): Lessons Learnt from EUROSEISTEST. *Bull. Earthq. Eng.* 16, 2311–2336. doi:10.1007/s10518-017-0098-6

Maufroy, E., Chaljub, E., Theodoulidis, N. P., Roumelioti, Z., Hollender, F., Bard, P. Y., et al. (2017). Source-Related Variability of Site Response in the Mygdonian Basin (Greece) from Accelerometric Recordings and 3D Numerical Simulations. *Bull. Seismol. Soc. Am.* 107, 787–808. doi:10.1785/0120160107

Perron, V. (2017). Apport des enregistrements de séismes et de bruit de fond pour l'évaluation site-spécifique de l'aléa sismique en zone de sismicité faible à modérée. PhD thesis. Grenoble, France: Université Grenoble Alpes.

Perron, V., Gélis, C., Froment, B., Hollender, F., Bard, P.-Y., Cultrera, G., et al. (2018). Can Broad-Band Earthquake Site Responses Be Predicted by the Ambient Noise Spectral Ratio? Insight from Observations at Two Sedimentary Basins. *Geophys. J. Int.* 215, 1442–1454. doi:10.1093/gji/ggy355

Perron, V., Laurendeau, A., Hollender, F., Bard, P.-Y., Gélis, C., Traversa, P., et al. (2017). Selecting Time Windows of Seismic Phases and Noise for Engineering Seismology Applications: A Versatile Methodology and Algorithm. *Bull. Earthq. Eng.* 16, 2211–2225. doi:10.1007/s10518-017-0131-9

Régnier, J., Cadet, H., and Bard, P. Y. (2016). Empirical Quantification of the Impact of Nonlinear Soil Behavior on Site Response. *Bull. Seismol. Soc. Am.* 106, 1710–1719. doi:10.1785/0120150199

Régnier, J., Cadet, H., Bonilla, L. F., Bertrand, E., and Semblat, J.-F. (2013). Assessing Nonlinear Behavior of Soils in Seismic Site Response: Statistical Analysis on KiK-Net Strong-Motion Data. *Bull. Seismol. Soc. Am.* 103, 1750–1770. doi:10.1785/0120120240

Sánchez-sesma, F. J. (1987). Site Effects on Strong Ground Motion. *Soil Dyn. Earthq. Eng.* 6, 124–132. doi:10.1016/0267-7261(87)90022-4

SIA (2020). SIA 261 Einwirkungen auf Tragwerke. Zurich, Switzerland (in German): Schweizerischer Ingenieur- und Architekentverein, Swiss Society of Engineers and Architects. Available at: http://shop.sia.ch/normenwerk/ingenieur/sia%20261/d/2020/D/Product/(Accessed April 1, 2022).

Thompson, E. M., Baise, L. G., Tanaka, Y., and Kayen, R. E. (2012). A Taxonomy of Site Response Complexity. *Soil Dyn. Earthq. Eng.* 41, 32–43. doi:10.1016/j.soildyn.2012.04.005

Zhu, C., Chávez-García, F. J., Thambiratnam, D., and Gallage, C. (2018). Quantifying the Edge-Induced Seismic Aggravation in Shallow Basins Relative to the 1D SH Modelling. *Soil Dyn. Earthq. Eng.* 115, 402–412. doi:10.1016/j.soildyn.2018.08.025

Zhu, C., Cotton, F., Kwak, D.-Y., Ji, K., Kawase, H., and Pilz, M. (2022). Within-Site Variability in Earthquake Site Response. *Geophys. J. Int.* 229, 1268–1281. doi:10.1093/gji/ggab481

**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 Perron, Bergamo 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_110.png" /></figure>

Check for
updates

# Seismic Risk Assessment of Multiple Cities: Biases in the Vulnerability Derivation Methods for Urban Areas With Different Hazard Levels

Maria Camila Hoyos¹* and Andres Felipe Hernandez²

¹Independent Seismic Risk Consultant, Bogota, Colombia, ²Institut des Sciences de la Terre, Université Grenoble Alpes, Grenoble, France

Previously, it has been shown that probabilistic seismic risk assessments (PSRAs) at urban scale present important discrepancies when compared with analyses conducted using methodologies from regional or national PSRA. However, conducting site-specific urban-scale analyses for a considerable number of cities may not be feasible due to limitations in time, resources, and in some cases availability of information, and thus more general models or methodologies are used. This brings into the picture the importance of identifying and quantifying the possible biases, discrepancies, and uncertainties when using different methods, both in the hazard and vulnerability components. Regarding the latter, several sources of uncertainty and biases have been identified in 1) the selection of ground motion records, either by using a general pool of records such as the ones from FEMA P695 or by performing a site- or hazard-specific analysis that requires a significant effort, especially in areas with a poor history of seismic instrumentation and even more in regions with no evidence of previous PSRA at all; and 2) the fragility or vulnerability derivation nonlinear dynamic methods: incremental dynamic analysis [IDA], cloud analysis [CA], and multi-stripe Analysis [MSA], among others. Focusing on these sources of uncertainty and bias, and with the challenge to bring solutions for places with scarce information, in this study, we aim to explore the use of different vulnerability derivation assumptions for the three principal cities of Colombia: Bogota, Medellín, and Cali, where most of the economic growth is concentrated. This considers the different seismic hazard levels and tectonic environment contributions in each city. Afterward, a comparison between the results of the analysis without the hazard-specific record selection and the site-specific one for each city is performed to establish the cases in which the former is applicable without adding more biases or uncertainties in the process.

**Keywords:** seismic risk assessment, fragility function derivation, record selection, multi-site analysis, urban risk

## 1 INTRODUCTION

Probabilistic seismic risk assessment (PSRA) at the urban scale has become a critical tool for local governments, disaster risk reduction offices, and insurers in the definition of long-term actions for risk management in cities such as the response planning after an event, the definition of city-specific insurance plans, and the overall urban planning. This considering that cities concentrate 80% of the GDP generation

**OPEN ACCESS**

**Edited by:**

Dimitris Pitilakis,
Aristotle University of Thessaloniki,
Greece

**Reviewed by:**

Anastasia Nekrasova,
Institute of Earthquake Prediction
Theory and Mathematical Geophysics
(RAS), Russia

M. Amin Haniri-Ardebili,
National Institute of Standards and
Technology (NIST), United States

**\*Correspondence:**

Maria Camila Hoyos
mc.hoyos.ramirez@gmail.com

**Specialty section:**

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

**Received:** 31 March 2022

**Accepted:** 17 June 2022

**Published:** 04 August 2022

**Citation:**

Hoyos MC and Hernandez AF (2022)
Seismic Risk Assessment of Multiple
Cities: Biases in the Vulnerability
Derivation Methods for Urban Areas
With Different Hazard Levels.
Front. Earth Sci. 10:910118
doi: 10.3389/feart.2022.910118
worldwide and 56% of the global population, with an estimation to increase to 68% by mid-century (United Nations, 2019).

In recent decades, earthquake disaster risks in cities have increased mainly due to the high rate of urbanization, lack of urban planning, and inadequate or uncontrolled construction practices, among others. Because of this, special attention has been given to the analysis of urban environments with the inclusion of specific targets, indicators, or goals for urban resilience in the sustainable development goals (Goal 11), the Sendai framework for disaster risk reduction (Indicator E2), and The New Urban Agenda (United Nations, 2015; United Nations Office for Disaster Risk Reduction, 2015; United Nations, 2016).

Considering this new demand for urban-specific models, many studies have downsized the methodologies followed for global, regional, or national PSRAs, without much consideration of the implications and biases generated in the process. From the three main inputs of a seismic risk model, hazard, exposure, and vulnerability, this is particularly true for the vulnerability component. As stated by Kohrangi et al. (2017b), when deriving vulnerability models for portfolios in multiple cities, even if they have different hazard characteristics, the common engineering practice is to use a pool of general records regardless of their consistency with the hazard at each site, which has been shown to lead to potentially biased risk estimates, even under a “sufficient” intensity measure (IM). For this reason, some studies have proposed methods to reduce this bias by improving the record selection (Jayaram et al., 2011; Lin et al., 2013b) by modifying the general pool of records to include in some way the expected spectral shape of the site (Haselton et al., 2011) or by using more “sufficient” and “efficient” IMs as the average spectral acceleration—AvgSa (De Biasio et al., 2014; Eads et al., 2015; Kohrangi et al., 2017a).

Among those studies, Kohrangi et al. (2017b) established that identical buildings should be characterized by different vulnerability functions at different sites because the magnitude of the earthquakes around the site and the distance to the nearby faults, among others, can modify the type of ground motions that could occur at each site. In this regard, there has been an ever-growing group of studies comparing the results in terms of fragility or vulnerability curves or even risk estimates between analyses considering general pools of records such as the FEMA P-695 far-field set (Federal Emergency Management Agency, 2009) with the ones using more detailed site-specific record selection methodologies (Lin et al., 2013a; Kohrangi et al., 2017a; Kohrangi et al., 2020). Overall, it has been shown that there is no general rule to follow and that the use of a general pool of records or a specific one should be based on more careful consideration of the characteristics of the hazard in the site, the expected ground motions, the structures to analyze, and the availability of resources.

Another important source of uncertainty comes from using different fragility derivation methodologies as each one considers different assumptions, follows different methodologies, and needs specific inputs that could add to the uncertainties of the fragility derivation process. Among the nonlinear dynamic models that use time-history analyses within the analytical fragility derivation methods, incremental dynamic analysis, multi-stripe analysis, and cloud analysis have been given much attention. As in the previous case, some studies have conducted comparisons between the

different derivation methods, illustrating the variabilities that can be encountered in the process and giving recommendations on when each one should be considered based on the limitations and benefits of each method (Jalayer, 2003; Mackie and Stojadinovi, 2005; Baker, 2015; Jalayer et al., 2017). However, as with its ground motion selection counterpart, many studies are based on specific case studies and thus their considerations apply only in similar scenarios to those that have been studied. For this reason, there is a need to study in more regions and under different assumptions if the conditions to use one or the other are fulfilled or if there needs to be special care under some circumstances when using one procedure over the other.

Considering this and the small number of studies that have been presented analyzing nonengineered structures located in urban environments exposed to different hazard levels and with different tectonic regime contributions, we aim to provide a reference to establish if the methods currently used in the engineering practice are applicable under these conditions. For the case study, we will derive fragility functions for the three phenotype buildings previously presented in Hoyos and Hernández (2021) for the three principal cities of Colombia, Bogota, Medellin, and Cali, where each one is exposed to a different hazard level and different tectonic regime contributions. In the following sections, we will present the methodologies that will be used and compared within the scope of this research, followed by the specifics of the case study included in a more detailed manner, the characteristics of the sites, and the structures to be analyzed. Later, the record selection process is described, and the fragility curves derived following the different record selection procedures and nonlinear dynamic fragility derivation methodologies are shown and compared, establishing the biases incurred in each case and giving some recommendations on when the analyses of some of the structures should be performed.

# 2 METHODOLOGY

## 2.1 Record Selection

A set of records must be used to perform the fragility analysis and for this, two different approaches were determined, using a hazard-specific selection through the conditioned spectrum (CS) as well as using the most common set of records from the FEMA P695 (Federal Emergency Management Agency, 2009).

### 2.1.1 Conditioned Spectrum

The CS method (Jayaram et al., 2011) was used to select the ground motion records by computing the mean scenario (i.e., mean magnitude, M, mean distance, R, and mean epsilon, ε) that best represents the site of analysis for the selected intensity measure level types (Harmsen, 2001) and tectonic regimes. There are several variations of this approach (Lin et al., 2013), but what we followed is the so-called “approximate” method where the CS is estimated using the mean values of magnitude and distance as well as the logic-tree weights of the ground motion prediction models. To perform this method, correlation models for the tectonic regimes of active shallow crust and subduction were considered (Jayaram and Baker, 2009; Jaimes and Candia, 2019).
<figure><img src="image_111.png" /><figcaption>**FIGURE 1** | Location of Bogota, Medellin, and Cali as well as some characteristics of the seismicity in the country. Dash-dotted lines: active faults traces. Continuous lines: boundary tectonic plates.</figcaption></figure>

<table><caption>**TABLE 1** | Main characteristics of the different typologies considered in the analysis.</caption><thead><tr><th>Typology</th><th>Inter-storey height</th><th>Structural period [s]</th><th>Γ</th><th>Yield drift [%]</th><th>Ultimate drift [%]</th><th>References</th></tr></thead><tbody><tr><td>MUR-H2</td><td>2.4</td><td>0.25</td><td>1.2</td><td>0.1</td><td>0.5</td><td>Acevedo et al. (2017)</td></tr><tr><td>CR-H4</td><td>2.95</td><td>0.5</td><td>1.3</td><td>0.15</td><td>0.8</td><td>Sinisterra (2017)</td></tr><tr><td>CR-H8</td><td>2.95</td><td>1.0</td><td>1.3</td><td>0.15</td><td>0.8</td><td>Sinisterra (2017)</td></tr></tbody></table>

Adapted from Table 4 from (Hoyos and Hernández, 2021).

The advantage of this method is the definition not only of the target intensities but also the spectral shape of all the records selected, which in comparison to other scaling methods provides a much more accurate and realistic spectral shape.

### 2.1.2 FEMA

The site-independent far-field set of 44 records from the FEMA P695 was used as a comparison method because it is a common practice to perform nonlinear dynamic analysis with them, but

without the proper assessment of which cases, these scenarios are suitable to be used.

## 2.2 Nonlinear Dynamic Methods for Fragility Curve Derivation

Among the literature on the topic, there are three main nonlinear dynamic methods for the derivation of analytical fragility models: the incremental dynamic analysis—IDA (Vamvatsikos and
<figure><img src="image_112.png" /><figcaption>**FIGURE 2** | Capacity curves for the three typologies considered. Left: unreinforced masonry two-storeys (MUR-H2). Center: reinforced concrete four-storeys (CR-H4). Right: reinforced concrete eight-storeys (CR-H8).</figcaption></figure>

<figure><img src="image_113.png" /><figcaption>**FIGURE 3** | Tectonic regime contribution in the three cities: Bogotá (top), Medellín (center), and Cali (bottom) for the three-building types: two-story unreinforced masonry (left), four-story reinforced concrete (center), and eight-story reinforced concrete (right), for different return periods.</figcaption></figure>
<figure><img src="image_114.png" /><figcaption>**FIGURE 4** | Uniform hazard spectra for 10% probability of exceedance (poe) for Sa (0.25 s) left, Sa (0.5 s) center, and Sa (1.0 s) right for the three cities: Bogotá (top), Medellín (center), and Cali (bottom).</figcaption></figure>

Cornell, 2002), the multi-stripe analysis—MSA (Jalayer, 2003; Jalayer and Cornell, 2009), and the cloud analysis—CA (Cornell et al., 2002) which are the ones considered in this study, clarifying that the adaptative IDA version—AIDA (Lin and Baker, 2013) was considered instead of the IDA when the CS selected set of records was used, and the Bayesian CA (Jalayer et al., 2015) following the procedure reported in Martins and Silva (2020) was used instead of the original CA. For the derivation of fragility curves, the IDA considers a fixed suite of records scaled successively to higher intensity levels, where each ground motion in the suite is scaled until it causes structural collapse. This method was the one followed when deriving the fragility curves using the FEMA P695 far-field record set. The AIDA follows a similar principle but considers different suites of records that are scaled to multiple IMs while they fulfill the hazard characteristics at said IM level, and thus not one but many curves are generated across the multiple bins selected by connecting the records that are repeated from one level to the other.

The MSA, on the other hand, establishes some discrete intensity levels and includes in each a different set of ground motions selected based on site-specific characteristics considering

methods such as the CS. This is done since the representative ground motions (and thus their properties) change at each IM level and with the objective to avoid one of the criticisms of the IDA methodology, the over-scaling of the records. This consideration of the local hazard characteristics makes it the most accurate or reliable method if a significant number of records is available for each IM level. Finally, the CA considers performing a regression fitting between the IM and the engineering demand parameter—EDP—in the logarithmic space. Its updated Bayesian version goes one step ahead and includes a differentiation in the regression between the cases where the structure exceeds the collapse limit state and indeed collapses and where it does not.

## 3 CASE STUDY

Colombia is located in the northwestern part of South America where the interaction between the tectonic plates of Nazca, the Caribbean, and South American, as well as the Panama and North Andes blocks produces deformations that trigger events in the country of different characteristics, some of them with destructive
<table><caption>**TABLE 2** | Percentages of the number of events selected for each city, structural period, and tectonic regime type.</caption><thead><tr><th rowspan="2">Tectonic regime type</th><th colspan="3">Bogota</th><th colspan="3">Medellín</th><th colspan="3">Cali</th></tr><tr><th>Sa (0.25)</th><th>Sa (0.5)</th><th>Sa (1.0)</th><th>Sa (0.25)</th><th>Sa (0.5)</th><th>Sa (1.0)</th><th>Sa (0.25)</th><th>Sa (0.5)</th><th>Sa (1.0)</th></tr></thead><tbody><tr><td>Active shallow crust</td><td>100%</td><td>100%</td><td>100%</td><td>37.3%</td><td>34.6%</td><td>36.0%</td><td>7.5%</td><td>12.9%</td><td>14.5%</td></tr><tr><td>Subduction interface</td><td>0%</td><td>0%</td><td>0%</td><td>46.1%</td><td>58.8%</td><td>57.3%</td><td>19.5%</td><td>65.4%</td><td>72.9%</td></tr><tr><td>Subduction intraslab</td><td>0%</td><td>0%</td><td>0%</td><td>16.7%</td><td>6.6%</td><td>6.7%</td><td>73.0%</td><td>21.8%</td><td>12.6%</td></tr></tbody></table>

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

consequences in the past. This complex tectonic environment generates several geologic faults along the three branches of the Andes Mountains from Ecuador to Venezuela.

Then, to study the effect of different hazard intensity levels and the contributions of different tectonic regimes and thus different hazard scenarios, the three main cities in Colombia, Bogotá, Medellín, and Cali, are included in the case study (Figure 1). These are the cities with the largest contribution to the GDP in the country and account for more than 20% of the total population of the country (DANE, 2018). All of them are located in the central part of the country; Bogota is in the eastern Cordillera, while

Medellín and Cali are between the central and western branches. The proximity of these cities with some traces of active shallow faults as well as the subduction events generated on the Pacific coast, and the history of destructive events in the past such as the M<sub>w</sub> 6.1 Armenia earthquake on 25 January 1999, is sufficient enough to develop and execute plans for the prevention and mitigation of seismic risk in the country.

The seismic hazard model (SHM) used is provided by the Colombian Geological Survey which was made in collaboration with the Global Earthquake Model Foundation and takes into account the seismological and geologic studies in the country as
<figure><img src="image_116.png" /></figure>

**FIGURE 6** | Conditioned spectra (CS), response spectra from the records selected, and uniform hazard spectra (UHS) for subduction intraslab at the intensity measure level of 10% probability of exceedance in 50 years for Sa (0.25 s) left, Sa (0.5 s) center, and Sa (1.0 s) right for the cities: Medellín (top) and Cali (bottom). Continuous blue thick line: median CS. Dash-lines: 2.5 and 97.5 CS percentiles. Continuous gray thin lines: records selected.

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

**FIGURE 7** | Conditioned spectra (CS), response spectra from the records selected, and uniform hazard spectra (UHS) for subduction interface at the intensity measure level of 10% probability of exceedance in 50 years for Sa (0.25 s) left, Sa (0.5 s) center, and Sa (1.0 s) right for the cities: Medellín (top) and Cali (bottom). Continuous blue thick line: median CS. Dash-lines: 2.5 and 97.5 CS percentiles. Continuous gray thin lines: records selected.

well as a homogeneous earthquake catalog from the national seismological network updated until 2020. A detailed study of the intensity levels registered by the national accelerometric network allowed the authors of the SHM to establish a selection of ground

motion prediction models that, using a logic tree, take into account the uncertainty in the intensity levels estimation. The SHM has four different tectonic regimes that contribute to the seismic hazard levels in the country: active shallow crust,
<figure><img src="image_118.png" /><figcaption>**FIGURE 8** | Scaled records from FEMA to the intensity level of 0.45 g for the three different intensity measures of Sa (0.25 s) left, Sa (0.5 s) center, and Sa (1.0 s) right. Gray lines: scaled records. Continuous blue line: mean. Dashed blue lines: 2.5 and 9.7 percentiles.</figcaption></figure>

<table><caption>**TABLE 3** | Damage states used in the derivation of fragility curves.</caption><thead><tr><th>DS1: Slight</th><th>DS2: Moderate</th><th>DS3: Extensive</th><th>DS4: Collapse</th></tr></thead><tbody><tr><td>0.75Sdy</td><td>0.5Sdy + 0.33Sdu</td><td>0.25Sdy + 0.67Sdu</td><td>Sdu</td></tr></tbody></table>

subduction interface, subduction intraslab, and deep seismicity. However, for the cities of Bogota, Medellin, and Cali, the seismic sources of active shallow crust and subduction are the most relevant. For more details on the seismic hazard model, the reader could refer to (SGC, 2018; Arcila et al., 2020).

Three different structural types, representative of some of the most common vulnerable building classes in Colombia and ranging from the short to the medium structural period, are considered in the case study: the two-storey unreinforced masonry structures (MUR-H2) and the four-storey (CR-H4) and eight-storey (CR-H8) pre-code reinforced concrete buildings. The structural characteristics of the phenotype buildings were taken considering the local 1 type of structures reported in Hoyos and Hernández (2021), which are presented in **Table 1**.

Taking the values of the characteristics reported in **Table 1** and using Eqs 1, 2 taken from Villar-Vega et al., (2017), the capacity curves of the equivalent single degree of freedom—SDOF—model in terms of Sa and Sd were derived. These curves are presented in **Figure 2**.

$$S_d = \frac{N_{storeys} \times h_{storey} \times \theta_{global}}{\Gamma} \quad (1)$$

(Villar-Vega et al., 2017)

$$Sa_y = Sd_y \left(\frac{2\pi}{T_y}\right) \quad \text{and} \quad Sa_u = C*Sa_y, \quad (2)$$

(adapted Villar-Vega et al., 2017) where C = 1.0 for the two-story masonry structure and C = 0.8 for the concrete structures.

## 4 RECORD SELECTION

The seismic hazard levels or target levels used to scale the ground motion records were computed for the three cities using the seismic

hazard model by Arcila et al. (2020), which is developed in OpenQuake software (Pagani et al., 2014; Silva et al., 2014). The contributions per tectonic environment were estimated for the sites of analysis at different return periods, and as a result, the dominant sources that contribute to the hazard of those cities come from the active shallow crust and subduction (interface and intraslab) tectonic environments. **Figure 3** presents the contribution to the hazard for each city and structural period studied at ten different intensity measure levels or return periods (from 73 to 100,000 years). As it can be seen, the contribution of deep seismicity sources is minimum for all cities; so from now on, it will not be taken into account in further analysis. Likewise, in Bogota, only sources of the active shallow crust will be considered because they represent almost all the contributions to the hazard. On the contrary, in Cali, the higher contribution comes from the subduction sources, but a small proportion of the active shallow crust is still present. Medellin, on the other hand, presents a balanced proportion between sources of active shallow crust and subduction.

Ten different intensity measure levels (IMLs) were established to obtain the target intensity values at the intensities measured close to the fundamental periods, T1, of the systems presented in **Table 1**, to increase the efficiency and ensure lower uncertainties in the response predictions (Luco and Cornell, 2007). **Figure 4** shows the hazard curves of the three cities for the intensity measures (IMs) considered, Sa(0.25s), Sa(0.5s), and Sa(1.0s). The figure depicts in each case the total hazard curve as well as the contribution of active shallow crust (ASC), subduction interface (SUB Interface), and subduction intraslab (SUB Inslab). As it can be seen, Cali is the city with the highest hazard levels for all the IMs due to the proximity to subduction sources on the Pacific coast, while Bogota as stated before only has an important contribution of active shallow sources.

It is worth mentioning that the target intensity levels were estimated per each case for the different tectonic environments and not for the total hazard, as well as considering a response on rock (soil type B). In the view of the authors, using target values much higher than that could occur according to the hazard model, leading to an overestimation of the intensity levels in the records selected. Therefore, the record selection was performed independently for active shallow crust and subduction (interface and intraslab).
<figure><img src="image_119.png" /><figcaption>**FIGURE 9** | Fragility curves derived using different methods: cloud (dotted), AIDA (dashed line), MSA (continuous line), and FEMA (dotted-dashed line); in the three cities: Bogotá (top), Medellín (center), and Cali (bottom); for the three-building types: two-story unreinforced masonry (left), four-story reinforced concrete (center), and eight-story reinforced concrete (right).</figcaption></figure>

The accelerograms were collected from several ground motion databases worldwide that include events for the tectonic environments presented in the analysis, such as Pacific Earthquake Engineering Research (PEER) NGA-West2 (Ancheta et al., 2014), NGA-Sub (Bozorgnia, 2020), Colombian Geological Survey (SGC, 2022), K-NET and KiK-net networks (NIED, 2022), National Seismological Service of Mexico (SSN-UNAM, 2022), and the SIBER-RISK strong motion database of Chilean earthquakes (Castro et al., 2020). A set of about 200 records was found per each case (i.e., city and structural period), ensuring that when no records were available for a

particular target intensity value, a maximum scaling factor of 5.0 was set as a threshold. Then, to unify all the records, we used the contributions of the hazard excluding the deep seismicity sources to create a set of about 1800 records in total for all the cases. **Table 2** presents the proportion of the number of events selected per case.

## 4.1 Conditioned Spectra

**Figure 5** shows the response spectrum from the records selected and scaled for the intensity measure level of 10% probability of exceedance (i.e., 475 years return period) in Bogota, Medellín,
<figure><img src="image_120.png" /></figure>

and Cali for the active shallow crust tectonic regime type, on each of the three intensity measures of interest, Sa (0.25 s), Sa (0.5 s), and Sa (1.0 s). As it can be seen from the figure, all the response spectra follow the shape of the conditioned spectrum (CS) and more or less are between the 2.5 to 97.5 percentiles of the CS. **Figures 6 and 7** present the same results as **Figure 5** but for the cases of subduction intraslab and subduction interface. In those cases, only records from Medellin and Cali were selected.

## 4.2 FEMA

The FEMA records were scaled by the mean of all the 44 ground motion records to different target intensity levels. **Figure 8** shows the response spectra of all the records scaled and the mean spectrum scaled to an intensity level of 0.45 g for the different intensity measures of Sa (0.25 s), Sa (0.5 s), and Sa (1.0 s). This procedure was repeated ten times from 0.15 to 1.5 g with increments of 0.15 g.

## 5 FRAGILITY DERIVATION

The SDOF models for each of the structural types whose capacity curves are presented in **Figure 2** were then subjected to two nonlinear time history analyses—NLTHA; the first used input ground motions those of the FEMA P695 far-field record set described in **Section 2.1.2** and the second one using the CS record set for each structural type and each city as described in **Section 2.1.1**. The maximum IM and EDP measures for each record were then recorded and later used to derive fragility functions using the methods explained in **Section 2.2**, considering the damage states reported in **Figure 6** by Martins and Silva (2020), which are presented in **Table 3**.

For each method, a different procedure was followed. For the cloud analysis, a linear regression is conducted between an IM and EDP values in the logarithmic space, differentiating the noncollapse and collapse cases by censoring the last ones. The median and standard deviation of each damage state can be calculated using the procedure described in Martins and Silva (2020), which is available

in the Vulnerability Modelers Toolkit open code (Martins et al., 2021). On the other hand, the IDA procedure follows the derivation of IDA curves that also provide an IM-EDP relationship. For calculating the mean and standard deviation, following an IM-based procedure, Eqs 3, 4 can be used for each damage state

$$ \theta = \frac{\sum_{i=1}^{n} \ln Sa^{DCR=1}}{n}, \qquad (3) $$

$$ \beta_{Sa|DS=dsi} = \sqrt{\frac{\sum_{i=1}^{n} (\ln Sa^{DCR=1} - \theta)^2}{n - 1}}, \qquad (4) $$

where DCR is the ratio between the recorded EDP over the DS. Finally, for the MSA procedure, considering the damage matrix, for each IM used as the target, it is possible to establish how many records can cause collapse (or for other damage states how many exceed the value) and using the maximum likelihood estimator, it is possible to fit lognormal parameters for each damage state. The comparison of the fragility curves obtained from each of the methods, for each structural type in each city, is presented in **Figure 9**. Additionally, the comparison of the MSA curves between the three cities is presented in **Figure 10**. The comparison of the fragility curves obtained from each of the methods, for each structural type in each city, is presented in **Figure 9**. Additionally, the comparison of the MSA curves between the three cities is presented in **Figure 10**.

## 6 ANALYSIS OF RESULTS

There is an important variability in the results obtained for the three studied building types shown in **Figures 9 and 10**. In the case of the eight-story reinforced concrete building (CR-H8), the derived fragilities seem to be unaffected by the record selection approach as can be observed in the right plot of **Figure 10**. In this case, the differences between the curves for the three cities compared with the one derived with the FEMA P695 records are insignificant, and also the differences between the fragility
<figure><img src="image_121.png" /><figcaption>**FIGURE 11** | Comparison of the CS and mean FEMA spectra for 10% probability of exceedance in 50 years for Sa (0.25 s) left, Sa (0.5 s) center, and Sa (1.0 s) right for the three cities: Bogotá (top), Medellín (center), and Cali (bottom).</figcaption></figure>

derivation methods are minimal, as shown in the three plots on the right of **Figure 9**. Thus, for this particular building type, it could be said that the site-specific record selection, which can be more time-consuming, could be neglected as the general pool of records will give accurate results.

On the opposite side, however, the results for the two-story unreinforced masonry structure (MUR-H2) are shown to be extremely sensitive to both the record selection procedure and the fragility curve derivation methodology. Considering the former, for the three cities under study, the fragilities derived using the FEMA P695 record set largely underestimate the capacity of the structures and thus can cause extremely conservative risk results, as they are the most fragile curves in the left plots of **Figure 9**. On the other hand, regarding the fragility curve derivation method, there is a large dispersion in the results in particular for the severe and collapse damage states. For all cases, the MSA method gives the less conservative fragilities, while the AIDA method tends to give the most conservative ones, having the CA landing usually between both. It is interesting to see that even when the AIDA and not the original IDA is used, the results still tend to be more conservative when using this procedure than the other two methods. However, at this point, it is important to mention that given the large IMs at which the records

needed to be scaled to reach collapse, there can be a considerable reduction in the number of available records that can be used for these high IMs without reaching the limit of 5 in the scale factors, and thus the results were expected to produce larger variabilities for the final damage states.

Regarding the four-story reinforced concrete structure CR-H4, there is considerable variability in the fragility curves derived using different methods for all damage states but the slight one, however, not as pronounced as in the MUR-H2 case. In this case, there is not an overall behavior present among all cities when comparing the site-specific curves with those derived from FEMA (as in the case of MUR-H2 and CR-H8). As it can be appreciated in the center plots of **Figure 9**, there is a similarity with the FEMA fragilities for the Medellín case; however, for Cali, the IDA using FEMA P695 far-field records tends to overestimate the capacity of the structures, while in Bogotá it underestimates it.

Considering this last point, an interesting result could be seen in the comparison between cities. As previously stated, both the MUR-H2 and the CR_H4 buildings present very different fragilities among cities (see **Figure 10**). In the CR-H4 case, this behavior seems to be counterintuitive when looking at the total hazard curves presented in **Figure 4**, as the total hazard in
Cali is shown to be higher than that of Bogotá and Medellín. At this point, it is important to mention that within the selection of records using the CS procedure, given that tectonic regime-specific records were searched for, the conditioning intensity from the spectrum was that of the specific tectonic regime spectrum and not the one from the total spectrum. This considering that within the hazard model that was used, the sources of each tectonic regime cannot experience or produce values above those specifically calculated from them, and thus using the total hazard spectrum instead of the tectonic regime specific one could cause an overestimation of the demand of the site when compared with its known hazard. In this way, based on the percentages of records by the contribution of the tectonic regime in each city, that were obtained considering the disaggregation in each site and that are presented in Table 2 and looking at Figure 4, it can be seen that in some cases, the contribution of the modeled active shallow crust in Bogotá can be higher than any of the contributions of a specific regime type in Medellín and Cali. Considering this issue, further studies should be conducted to see the sensitivity and accuracy of fragility curves derived in regions with contributions of multiple sources to see if there is a considerable overestimation when assuming the total hazard spectra instead of that of each specific source.

Finally, to try to understand the overall behavior of the fragility curves, an analysis of the mean conditional spectra for the site-specific record sets and the FEMA recordset was conducted, looking into the spectral shapes of the mean response for each. **Figure 11** shows the mean conditional spectra representing 10% of exceedance in 50 years for all cases. It is very interesting to see that the spectral shape of the CS in all cities for the CR-H8 structure is very similar to that of the FEMA P695 far-field recordset. This agreement in the spectra can be the cause of the similarities of the fragility curves in this case and brings up a possible analysis that can be conducted to see if the FEMA P695 recordset can indeed be used in the derivation of fragility curves at a site.

Continuing with the analyses of **Figure 11**, for the MUR-H2 case, it can be seen that in all cases, the mean spectra from the FEMA dataset differ in spectral shape from that of any of the site-specific CS in the cities and tend to be higher after the structural period of analysis (0.25s). This is why for the larger damage states, when the structures start to present cracks and elongate their period, the FEMA fragilities are more fragile than the other ones. Finally considering the CR-H4 case, a very similar spectral shape can be seen for the Medellín case (the reason why the fragilities were also so similar in that case), while for Bogotá, the mean for the FEMA records after the structural period (0.5 s) tends to fall while the ones in Cali rise.

# 7 CONCLUSION

There is no general rule to follow to choose either a record selection approach or a fragility derivation method as it should be evaluated in every specific case based on a more careful consideration of the characteristics of the hazard at the site, the expected ground motions, the structures to analyze, and the availability of resources. However, it could be seen that the

spectral shape does seem to be an important parameter in the record selection method and the set of records that should be included in the analysis. Likewise, the ground motion prediction model (GMPM) plays an important role in the spectral shape of the CS and should be given more attention; in particular in sites where there is no locally derived GMPM. In this way, the use of a general set of records such as the FEMA P-695 far-field set seems appropriate for the fragility analysis of structures with medium-to-high fundamental periods. On the contrary, for a short fundamental period of vulnerable building types, a more detailed record selection should be performed as the FEMA P-695 most likely underestimates the capacity of the structures.

The use of target values much higher than what could occur according to the hazard model leads to an overestimation of the intensity levels in the records selected. Based on the hazard model used, the sources of each tectonic regime cannot experience or produce values above those specifically calculated from them, and thus using the total hazard spectrum instead of the tectonic regime specific one could cause an overestimation of the demand of the site when compared with its known hazard.

With all of these in mind, further studies should be conducted to see the sensitivity and accuracy of 1) fragility curves derived in regions with contributions from multiple sources to see if there is a considerable overestimation when assuming the total hazard spectra instead of that of each specific source; 2) the use of different GMPM to compute the conditioned spectra and the impact on the spectral shape in the record selection; 3) the consideration of using different methods in the computation of the conditioned spectra, with one or multiple GMPMs, different logic tree weights, different disaggregation weights, and by including all the possible scenarios (magnitudes and distances).

# DATA AVAILABILITY STATEMENT

The original contributions presented in the study are included in the article; further inquiries can be directed to the corresponding author.

# AUTHOR CONTRIBUTIONS

AH CS, record selection 50% MH NLTHA, Derivation fragility curves 50%.

# FUNDING

This work was funded by the EU Horizon 2020 program under Grant Agreement Number 813137, ITN-MSCA URBASIS project.

# ACKNOWLEDGMENTS

The authors thank the Colombian Geological Survey for sharing the National Seismic Hazard Model of Colombia.
REFERENCES

Acevedo, A. B., Jaramillo, J. D., Yepes, C., Silva, V., Osorio, F. A., and Villar, M. (2017). Evaluation of the Seismic Risk of the Unreinforced Masonry Building Stock in Antioquia, Colombia. Nat. Hazards 86 (S1), 31–54. doi:10.1007/s11069-016-2647-8

Ancheta, T. D., Darragh, R. B., Stewart, J. P., Seyhan, E., Silva, W. J., Chiou, B. S.-J., et al. (2014). NGA-West2 Database. Earthq. Spectra 30 (3), 989–1005. doi:10.1193/070913EQS197M

Arcila, M., Garcia, J., Montejo, J., Eraso, J., Valcacer, J., Mora, M., et al. (2020). *Modelo Nacional De Amenaza Sísmica Para Colombia*. Servicio Geológico Colombiano. doi:10.32685/9789585279469

Baker, J. W. (2015). Efficient Analytical Fragility Function Fitting Using Dynamic Structural Analysis. Earthq. Spectra 31 (1), 579–599. doi:10.1193/021113EQS025M

Bozorgnia, Y. (2020). *Data Resources for NGA-Subduction Project*. NGA-Subduction Principal Investigator. Pacific Earthquake Engineering Research Center.

Castro, S., Benavente, R., Crempien, J. G., Candia, G., and de la Llera, J. C. (2020). A Consistently Processed Strong-Motion Database for Chilean Earthquakes. Seismol. Res. Lett.

Cornell, C. A., Jalayer, F., Hamburger, R. O., and FoutchFoutch, D. A. (2002). Probabilistic Basis for 2000 SAC Federal Emergency Management Agency Steel Moment Frame Guidelines. J. Struct. Eng. 128 (4), 526–533. doi:10.1061/(asce)0733-9445(2002)128:4(526)

DANE (2018). *Censo Nacional de Población y Vivienda - CNPV*. Base de Datos. Departamento Administrativo Nacional de Estadística. Available at: https://www.datos.gov.co/Estad-sticas-Nacionales/Censo-Nacional-de-Poblacion-y-Vivienda-CNPV-2018/qzc6-q9qw (Accessed February 15, 2022).

De Biasio, M., Grange, S., Dufour, F., Allain, F., and Petre-Lazar, I. (2014). A Simple and Efficient Intensity Measure to Account for Nonlinear Structural Behavior. Earthq. Spectra 30 (4), 1403–1426. doi:10.1193/010614EQS006M

Eads, L., Miranda, E., and Lignos, D. G. (2015). Average Spectral Acceleration as an Intensity Measure for Collapse Risk Assessment. Earthq. Eng. Struct. Dyn. 44 (12), 2057–2073. doi:10.1002/eqe.2575

Federal Emergency Management Agency (2009). Quantification of Building Seismic Performance Factors. Available at: www.ATCouncil.org (Accessed January 31, 2022).

Harmsen, S. C. (2001). Mean and Modal in the Deaggregation of Probabilistic Ground Motion. Bull. Seismol. Soc. Am. 91 (6), 1537–1552. doi:10.1785/0120000289

Haselton, C. B., Baker, J. W., Liel, A. B., and Deierlein, G. G. (2011). Accounting for Ground-Motion Spectral Shape Characteristics in Structural Collapse Assessment through an Adjustment for Epsilon. J. Struct. Eng. 137 (3), 332–344. doi:10.1061/(asce)st.1943-541x.0000103

Hoyos, M. C., and Hernández, A. F. (2021). Impact of Vulnerability Assumptions and Input Parameters in Urban Seismic Risk Assessment. Bull. Earthq. Eng. 19 (11), 4407–4434. doi:10.1007/s10518-021-01140-x

Jaimes, M. A., and Candia, G. (2019). Interperiod Correlation Model for Mexican Interface Earthquakes. Earthq. Spectra 35 (3), 1351–1365. doi:10.1193/080918EQS200M

Jalayer, F., and Cornell, C. A. (2009). Alternative Non-linear Demand Estimation Methods for Probability-Based Seismic Assessments. Earthq. Eng. Struct. Dyn. 38 (8), 951–972. doi:10.1002/eqe.876

Jalayer, F., De Risi, R., and Manfredi, G. (2015). Bayesian Cloud Analysis: Efficient Structural Fragility Assessment Using Linear Regression. Bull. Earthq. Eng. 13 (4), 1183–1203. doi:10.1007/s10518-014-9692-z

Jalayer, F., Ebrahimian, H., Miano, A., Manfredi, G., and Sezen, H. (2017). Analytical Fragility Assessment Using Unscaled Ground Motion Records. Earthq. Eng. Struct. Dyn. 46 (15), 2639–2663. doi:10.1002/eqe.2922

Jalayer, F. (2003). *Direct Probability Seismic Analysis: Implementing Non-linear Dynamic Assessment*. Stanford University.

Jayaram, N., and Baker, J. W. (2009). Correlation Model for Spatially Distributed Ground-Motion Intensities. Earthq. Eng. Struct. Dyn. 38 (15), 1687–1708. doi:10.1002/eqe.922

Jayaram, N., Lin, T., and Baker, J. W. (2011). A Computationally Efficient Ground-Motion Selection Algorithm for Matching a Target Response Spectrum Mean and Variance. Earthq. Spectra 27 (3), 797–815. doi:10.1193/1.3608002

Kohrangi, M., Bazzurro, P., Vamvatsikos, D., and Spillatura, A. (2017a). Conditional Spectrum-Based Ground Motion Record Selection Using Average Spectral Acceleration. Earthq. Eng. Struct. Dyn. 46 (10), 1667–1685. doi:10.1002/eqe.2876

Kohrangi, M., Vamvatsikos, D., and Bazzurro, P. (2017b). Site Dependence and Record Selection Schemes for Building Fragility and Regional Loss Assessment. Earthq. Eng. Struct. Dyn. 46 (10), 1625–1643. doi:10.1002/eqe.2873

Kohrangi, M., Vamvatsikos, D., and Bazzurro, P. (2020). Multi-Level Conditional Spectrum-Based Record Selection for IDA. Earthq. Spectra 36 (4), 1976–1994. doi:10.1177/8755293020919425

Lin, T., and Baker, J. W. (2013). "Introducing Adaptive Incremental Dynamic Analysis: A New Tool for Linking Ground Motion Selection and Structural Response Assessment," in 11th International Conference on Structural Safety & Reliability, New York, NY, 805–811.

Lin, T., Harmsen, S. C., Baker, J. W., and Luco, N. (2013). Conditional Spectrum Computation Incorporating Multiple Causal Earthquakes and Ground-Motion Prediction Models. Bull. Seismol. Soc. Am. 103 (2A), 1103–1116. doi:10.1785/0120110293

Lin, T., Haselton, C. B., and Baker, J. W. (2013a). Conditional Spectrum-Based Ground Motion Selection. Part II: Intensity-Based Assessments and Evaluation of Alternative Target Spectra. Earthq. Eng. Struct. Dyn. 42, 1867–1884. doi:10.1002/eqe.2303

Lin, T., Haselton, C. B., and Baker, J. W. (2013b). Conditional Spectrum-Based Ground Motion Selection. Part I: Hazard Consistency for Risk-Based Assessments. Earthq. Eng. Struct. Dyn. 42 (12), 1847–1865. doi:10.1002/eqe.2301

Luco, N., and Cornell, C. A. (2007). Structure-Specific Scalar Intensity Measures for Near-Source and Ordinary Earthquake Ground Motions. Earthq. Spectra 23 (2), 357–392. doi:10.1193/1.2723158

Mackie, K. R., and Stojadinovi, B. (2005). "Comparison of Incremental Dynamic, Cloud, and Stripe Methods for Computing Probabilistic Seismic Demand Models," in Structures Congress 2005 (ASCE Library). doi:10.1061/40753(171)184

Martins, L., and Silva, V. (2020). Development of a Fragility and Vulnerability Model for Global Seismic Risk Analyses. Bull. Earthq. Eng. 19, 6719–6745. doi:10.1007/s10518-020-00885-1

Martins, L., Silva, V., Crowley, H., and Cavalieri, F. (2021). Vulnerability Modellers Toolkit, an Open-Source Platform for Vulnerability Analysis. Bull. Earthq. Eng. 19 (13), 5691–5709. doi:10.1007/s10518-021-01187-w

NIED (2022). NIED K-NET, KiK-Net. Available at: https://www.doi.org/10.17598/NIED.0004 (Accessed November 15, 2021).

Pagani, M., Monelli, D., Weatherill, G., Danciu, L., Crowley, H., Silva, V., et al. (2014). OpenQuake Engine: An Open Hazard (And Risk) Software for the Global Earthquake Model. Seismol. Res. Lett. 85 (3), 692–702. doi:10.1785/0220130087

SGC (2018). Modelo de Amenaza Sismica de Colombia. Available at: https://amenazasismica.sgc.gov.co/ (Accessed November 15, 2021).

SGC (2022). *Catálogo de Aceleraciones.* Dirección de Geoamenazas. Servicio Geológico Colombiano.

Silva, V., Crowley, H., Pagani, M., Monelli, D., and Pinho, R. (2014). Development of the OpenQuake Engine, the Global Earthquake Model's Open-Source Software for Seismic Risk Assessment. Nat. Hazards 72 (3), 1409–1427. doi:10.1007/s11069-013-0618-x

Sinisterra, H. (2017). *Determinación Del Desempeño Sísmico Para Edificios En Pórticos de Concreto Reforzado Diseñados Bajo La Norma de Diseño Sismo Resistente C.C.C.S.R-84*. Cali: Universidad Javeriana.

SSN-UNAM (2022). *Catalogo de Sismos*. IDG, Servicio Sismologico Nacional.

United Nations Office for Disaster Risk Reduction (2015). Sendai Framework for Disaster Risk Reduction 2015-2030. Geneva, Switzerland: UNDRR.

United Nations (2015). *Transforming Our World: The 2030 Agenda for Sustainable Development*. New York, United States: General Assembly September 25, 2015.

United Nations (2016). New Urban Agenda. Available at: http://uploads.habitat3.org/hb3/NUA-English.pdf (Accessed February 20, 2022).
United Nations (2019). *World Urbanization Prospects the 2018 Revision*. New York, United States: Department of Economic and Social Affairs.

Vamvatsikos, D., and Cornell, C. A. (2002). *Incremental Dynamic Analysis. Earthq. Eng. Struct. Dyn.* 31 (3), 491–514. doi:10.1002/eqe.141

Villar-Vega, M., Silva, V., Crowley, H., Yepes, C., Tarque, N., Acevedo, A. B., et al. (2017). Development of a Fragility Model for the Residential Building Stock in South America. *Earthq. Spectra* 33 (2), 581–604. doi:10.1193/010716EQS005M

**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 Hoyos and Hernandez. 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_122.png" /></figure>

Check for updates

OPEN ACCESS

EDITED BY
Simone Barani,
University of Genoa, Italy

REVIEWED BY
Paolo Boncio,
University of Studies G. d'Annunzio
Chieti and Pescara, Italy
Snježana Markušić,
University of Zagreb, Croatia

*CORRESPONDENCE
Daniela Famiani,
daniela.famiani@ingv.it

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

RECEIVED 06 May 2022
ACCEPTED 15 July 2022
PUBLISHED 19 August 2022

CITATION
Famiani D, Cara F, Di Giulio G, Vassallo M and Milana G (2022), Detection of hidden faults within the Fucino basin from single-station ambient noise measurements: The case study of the Trasacco fault system.
*Front. Earth Sci.* 10:937848.
doi: 10.3389/feart.2022.937848

COPYRIGHT
© 2022 Famiani, Cara, Di Giulio, Vassallo and Milana. 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.

# Detection of hidden faults within the Fucino basin from single-station ambient noise measurements: The case study of the Trasacco fault system

Daniela Famiani<fnref n="1" />*, Fabrizio Cara<fnref n="1" />, Giuseppe Di Giulio<fnref n="2" />, Maurizio Vassallo<fnref n="2" /> and Giuliano Milana<fnref n="1" />

The presence of normal fault systems in central Italy, outcropping or hidden below Quaternary covers in intra-mountain basins, is the expression of the Neogene-Quaternary evolution of the area, characterized by an extensional tectonic regime following the fold and thrust structuring of the Apennine orogen. Italian urban settlements of central Italy are developed on hills or mountains but also in lowland areas, which are often set up in sedimentary basins. In this framework, urban centers found close to fault lines are common, with strong implications on the seismic risk of the area. In this work, we performed a dense seismological passive survey (88 single-station ambient noise measurements) and used the horizontal-to-vertical spectral ratio (HVNSR) technique to investigate hidden faults in the Trasacco municipality located in the southern part of the Fucino Basin (central Italy), where microzonation studies pointed out hypothetical fault lines crossing the urban area with the Apennine orientation. These hidden structures were only suggested by previous studies based on commercial seismic lines and aerial photogrammetry; their presence in the basin area is confirmed by our measurements. This case study shows the potentiality of using the HVNSR technique in fault areas to have a preliminary indication of anomalous behaviors, to be investigated later with specific geophysical techniques. Our approach can support microzonation studies whenever fault zones are involved, especially in urban areas or in places designated for future developments.

KEYWORDS

microzonation, HVNSR, hidden faults, lateral heterogeneities, subsoil reconstruction

<footnote n="1">¹Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Sismologia e Tettonofisica, Roma, Italy,</footnote>

<footnote n="2">²Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Sismologia e Tettonofisica, L'Aquila, Italy</footnote>
## Introduction

Seismic microzonation studies are often carried out at a basic level, sometimes leaving open scientific points that require further *ad hoc* investigations to be solved. In this study, we face one of these cases regarding the Fucino Basin (hereinafter FB), a large intra-mountain continental basin in central Italy. The FB was the greatest lake of central Italy until its complete drainage at the end of the 19th century and has a tectonic origin related to the presence of important normal fault systems. Here, we focus on the southern part of the FB, specifically close to Trasacco village (red square in Figure 1). In 2013, the Trasacco municipality (about 6,000 inhabitants) has been the object of a seismic microzonation study, and one of the main objectives was to clarify the position of the Trasacco fault, an active and capable NW-SE oriented normal structure following the elongation of the Vallelonga valley (see zoom in Figure 1 for the location). In fact, the issue of characterizing and zoning active and capable faults is crucial for urban planning purposes and also for ensuring safe conditions to the pre-existing buildings. Despite the large amounts of data (geognostic, geophysical, and geotechnical) collected by professional geologists during the microzonation activities, no evidence of the fault was found in the urbanized area. The spatial continuation of this fault in the northern direction of the basin area is also uncertain.

Hence, we decided to perform a geophysical study devoted to identifying the extension of the Trasacco fault in the basin area and other possible tectonic structures in it. In particular, we performed a wide campaign of single-station ambient noise recordings, using the well-known horizontal-to-vertical spectral ratio (HVNSR) technique as the method of analysis (Nakamura, 1989). Some of these data, which were recorded simultaneously at different locations, have been also used to roughly estimate the shear-wave velocities of the sedimentary layers. The HVNSR results are discussed under two different perspectives: first, we adopt a classical interpretation scheme to get a first-order reconstruction of the subsoil of the study area through the identification of the fundamental resonance frequency (f0) of the investigated sites; and second, we perform a directional analysis on HVNSR, for e.g., by rotating the horizontal components of the ground motion for sites close to tectonic elements. The purpose of this latter study was to identify a sort of signature related to specific complexities of the investigated sites and possibly linked to the presence of the Trasacco fault and all the other tectonic structures. The correlation between these observations and the geology is supported by the availability of other relevant previous studies in the Fucino Basin.

## Geological setting

The FB is an intra-mountain lacustrine basin with a tectonic origin. Some of the faults that originated in the basin are able to

generate large seismic events, such as the 13 January 1915 M 7.0 Avezzano earthquake, which completely destroyed the town of Avezzano and caused about 30,000 casualties. The FB is infilled by Upper Pliocene-Holocene lacustrine and alluvial deposits (Bard and Bouchon, 1980a; Bard and Bouchon, 1980b; Giraudi, 1988; Bosi et al., 1995), which unconformably overlie Meso-Cenozoic carbonate and Neogene terrigenous successions (Cavinato et al., 2002). Quaternary deposits also include alluvial fans made by the dismantling of the surrounding reliefs. FB is the largest intra-mountain endoreic depression of the Apennines, filled by a well-preserved thick Quaternary continental succession (Mannella et al., 2019), representing an important record to reconstruct the geological evolution of central Italy (volcanic eruptions and climate changes). According to Galadini et al. (1995) and Cavinato et al. (2002), from a neo-tectonic point of view, the evolution of the FB has been mainly affected by the activity of two systems of faults that border the northern and eastern sectors of the area. These two can be described as the “active” sides of the basin, and they can be easily distinguished, for the evidence of structural elements and outcropping geological units, from the southern and western parts of the basin which, conversely, are supposed to have a passive role in the genesis of the plain.

The total subsidence of the area may be defined as the sum of the hanging-wall subsidence generated by the normal fault and sediment compaction, and the regional subsidence or uplift (Doglioni et al., 1998).

From the analysis of surface geological data and the interpretation of seismic lines, Cavinato et al. (2002) divided the sedimentary infilling of FB in two different groups of stratigraphic units: Lower and Upper Units. The Lower Units (Upper Pliocene) crop out in the northern and north-eastern margins of the basin and consist of breccia, fluvial, and marginal to open-lacustrine deposits. The Upper Units (Lower Pleistocene-Holocene) are represented by marginal lacustrine/ fluvial deposits; thick coarse-grained fan-delta deposits are inter-fingered with fluvial-lacustrine deposits at the foot of the main relief. It is worth mentioning that some of the articles published for the area aimed at studying the characteristics of the recent units more in detail: general descriptions of continental deposits (Zarlenga, 1987), more specific analysis of morphotectonic (Bard and Bouchon, 1985; Blumetti et al., 1993) or geomorphological features, as alluvial fans (Frezzotti and Giraudi, 1992) or terraces (Accordi 1975; Messina 1996).

Focusing more on the structural elements present in FB, the main faults are (right side in Figure 1): the Tremonti-Celano-Aielli Fault (TCAF) (WSW-ENE) and S. Potito-Celano Faults (NW-SE) in the north; the Luco Fault (LF) in the west; and the Trasacco Fault (TF) and Villavallelonga Fault (VF), the Pescina-Celano Fault (PCF), and the Serrone Fault (NW-SE) in the southern and southeastern parts of the basin. Recently, Lanari et al. (2021) proposed a structural analysis of the FB making inferences on
<figure><img src="image_123.png" /></figure>

how sediment loading/unloading influences the dynamics of fault systems, demonstrating positive feedback between sedimentation and faulting.

An important contribution for reconstructing the FB and understanding the role played by every structural element in its general evolution is given by the interpretation of industrial seismic reflection profiles [location by Patruno and Scisciani, (2021) displayed in Figure 2] which can identify the subsurface geometry of outcropping faults and also image the subsoil layers and their thickness.

For some authors, most of the thickness of the fluvio-lacustrine succession (more than 800 m) (Cavinato et al., 2002; Patacca et al., 2008; Cara et al., 2011) is generated by the San Benedetto-Serrone-Gioia dei Marsi fault (hereinafter SSGF) system, located in the eastern margin of the FB. Some others (Patruno and Scisciani 2021), conversely, acknowledged the maximum thickness of the Quaternary deposits (~1,750 m) in two depocentral areas, for the combined activity of the TCAF and San Benedetto-Serrone-Gioia dei Marsi Fault (SSGF).

In addition to the contribution of seismic lines, some authors (Cella et al., 2021; Mancinelli et al., 2021) recently performed studies on the magnetic anomalies of the FB. Cella et al. (2021) not only carried out a 3D gravity model for evaluating the depth of the Meso-Cenozoic carbonate bedrock but also provided

additional constraints on the position of the known hidden faults. Both the aforementioned research works discuss the paleogeographic and paleotectonic evolutions of the FB confirming the presence of a major depocenter as the effect of the activity of the SSGF in the eastern part of the FB, and a second depocenter whose evolution is ascribable to the activity of the Trasacco fault, on its western side.

Mancinelli et al. (2021) reconstructed 2D gravity models along a grid of seven sections crossing the basin, in which is shown the influence of the Miocene flysch deposits in the observed residual general anomaly present in the FB. They also highlighted the thickness of Miocene deposits, as modeled in the sections. A higher thickness of flysch deposits (around 1,100 m) is found in the western part of the basin, consistent with a syn-orogenic emplacement of the deposits, and decreases toward the south and east.

Concerning the shallower part of the subsoil, paleo-seismological studies in FB started during the '80s (Giraudi 1986; Serva et al., 1986; Giraudi 1988; Galadini and Galli 1996; Michetti et al., 1996; Galadini et al., 1997; Galadini et al., 1997), and a collection of all the results is summarized in Galadini and Galli (1999). Amoroso et al. (2016) also presented a synoptic review of all the paleoseismic data collected for the FB, including information related to the slip
<figure><img src="image_124.png" /><figcaption>**FIGURE 2**
Map of previous investigations: the fault lines from the ITHACA database are in red, and boreholes from ISPRA and Ente Fucino are in green. The study area is delimited by a black box.</figcaption></figure>

rates of all the tectonic structures, deriving a possible overall extension rate across the basin of about 3–3.5 mm/yr.

The study area of this work is located in the Trasacco municipality, going from the urban area toward the central part of the basin. The target area is interesting due to the presence of alluvial fans in the urbanized area and of buried faults in the most rural area. Although there is no evidence of fault scarps in this area, aerial photos show a straight NW-SE-oriented lineament (Oddone 1915; Galadini and Galli 1999) which represents the contact between soils with different lithologies. This feature addressed the excavation of a huge number of hand boreholes to precisely locate the fault in the field. Paleo-seismological trenches were performed in the area by Galli et al. (2012); Galadini and Galli, (1996), together with many boreholes and trenches, studied the Trasacco fault and reconstructed the vertical displacement along its length and associated them to specific historical seismic events. Results of

their study highlighted the multiple activation of this fault during the Holocene: the most recent event is related to the 1915 Fucino earthquake and the previous one is dated to fifth–sixth A.D. The slip rate of the fault obtained from available data decreases toward its NW edge, reflecting the trend of its offsets toward the center of the basin.

# Methods

## Data

The single-station noise measurements were performed by using two different seismic stations: the Reftek-130 and Lennartz MarsLite digitizers modified by SARA electronic instruments connected to Lennartz 3D/5s velocimetric sensors and Terrabot (SARA electronic instruments) units, and all-in-one 24-bit
<figure><img src="image_125.png" /></figure>

digitizers with internal 4.5 Hz sensors. For all the stations, time synchronization was ensured by GPS antennas.

Overall, the survey consisted of 88 noise measurements performed in different time periods (zoom in Figure 2), of which 12 have been visually inspected to evaluate the quality of data, by the analysis of the Fourier spectra. For some sites with bad quality data, we have repeated the measurement in different time periods (eg., TB** and FN** measures located along Section 1 in Figure 3, zoomed location in Figure 11). The minimum length of recording was fixed at 90 min, in order to allow a good statistic of the results in the frequency domain of interest. The sampling rate was set to 200 or 250 sps. Some sets of measurements were collected simultaneously.

The positions of the measurement points were chosen with the aim to obliquely cross the hypothetical surface projection of the Trasacco faults (Apennine trend of the tectonic structures), as evidenced from the seismic lines' interpretation.

Overall, the final geometry of the survey can be seen as a set of transects trying to cross orthogonally to the supposed fault line. Actually, the survey covered the entire study area but with an irregular density and with special attention to the basin area in

the north of Trasacco where buried fault structures were suggested by previous studies. In fact, whenever a set of measurements ended, the following measurements were performed trying to refine the previous results. For this reason, the distance between stations of each transect was largely variable, ranging from 30 to 500 m.

One of the transects (corresponding to Section 2 in Figure 3) is 6 km long and largely exceeds the Trasacco area. The aim of this measurement line was to investigate and understand the general trend of the main impedance contrasts in the FB at a wider scale.

## HVNSR analyses

Single-station noise data were analyzed using Geopsy software (Wathelet et al., 2020) in terms of HVNSR curves (Nakamura, 1989). The analysis was also computed by rotating the horizontal components from 0 to 180° in steps of 10. The complete steps of analysis, developed within the SESAME project (Site Effects Using Ambient Excitations, SESAME 2004), consist of:
1) Applying an anti-trigger algorithm that selects 40s-long windows into the whole three-component recordings, in order to avoid short or anomalous transients often of anthropic origins. The chosen length of the windows, 40s, guaranteed a good frequency sampling according to the targets of the study;

2) removing mean and linear trends from each time window;

3) tapering the edges of the time windows;

4) rotating the horizontal components of the desired angle;

5) computing the Fourier spectra of the three components;

6) smoothing the Fourier spectra using the logarithmic Konno-Ohmachi algorithm (Konno and Ohmachi, 1998), with the bandwidth coefficient b equal to 40;

7) by visual inspection, identifying regular trends on the Fourier spectra and discarding windows affected by possible anthropic disturbances;

8) dividing the rotated horizontal spectra by the vertical spectrum for the selected windows;

9) geometrically averaging the spectra or the HVNSR obtained for each 40s-long window.

The classical computation of HVNSR does not imply the directional analysis. In this case, the two horizontal components are NS (0°) and EW (90°), and after points e) and f), the two horizontal spectra are simply merged using an arithmetic mean, in order to have a unique horizontal spectrum for each 40s-long window.

Results of the directional analysis of HVNSR are often shown in an image map, with frequency as x-axis, angle of rotation as y-axis, and the color of the map that is representative of the amplitude of HVNSR for each frequency/angle pair. Nevertheless, the same results can also be shown by plotting the rotated HVNSR curves all together. This kind of representation is useful to prove any possible anisotropy of the wavefield. The classical post-processing of HVNSR (considering only the merged NS and EW components) consists of looking at the shape of each curve and retrieving the main peak, also called the fundamental or resonance frequency (f0) of the selected site. This step is not always trivial. Sometimes, the peak is unique, very narrow, and with an amplitude greater than 2. Therefore, the resonance peak can be recognized easily. On the contrary, there are several peaks very often, or the main peak is very broad and/or weak. In these cases, the user experience on the interpretation of the results is really crucial because it requires a general overview of both the geological setting and the HVNSR curves on the surrounding area of the selected measurement point. Hence, we were able to associate a resonance frequency value for all stations, except for some flat HVNSR curves (Supplementary Appendix S1). The frequency peaks are selected between 0.2 and 4 Hz, discarding higher frequency values because they were considered irrelevant for the scope of this work.

A way of representing the results obtained by the HVNSR technique is to plot the resonance frequencies in a map (Figure 3) and interpolate them (Figure 4) to produce a contour plot (ISPRA, 1881; ITHACA Working Group, 2019; Mascandola et al., 2019). An inverse distance weighting (IDW) method, for e.g., a deterministic spatial interpolation model, was applied to interpolate f0 data in the map. The general premise of this method is that the attribute values of any given pair of points are related to each other, but their similarity is inversely proportional to the distance between the two locations (Lu and Wong, 2008).

Another way of visualizing information from HVNSR results is to consider the entire curves, not the resonance frequency alone. For each transect, we produced an image map interpolating the HVNSR curves along the distance (Figures 5, 6). The linear regression method was applied to data to find the best fitting line to plot, and then spatial coordinates (in meters) were projected on it. This kind of plot is helpful to imagine subsoil heterogeneities (Joyner et al., 1981; Famiani et al., 2020): the amplitude of spectral ratio curves is represented in a colored scale, while the XY axis shows the distance of the single noise measurement with respect to the beginning of the transect and the frequency values of HVNSR, respectively.

## Velocity profile estimation

Many of the noise measurements were collected in synchronous acquisition. Hence, the simultaneous recordings acted as a passive array of seismic stations, located in 1D (lines along transects) or 2D (when positioned in grids or irregular geometries) configurations. These kinds of data can then be used to retrieve the dispersion characteristics of the noise wavefield crossing the stations and finally estimating a velocity profile for the shallower layers of the area. Unfortunately, among all the possible simultaneous sets of stations, only the ones including stations from TB00 to TB08 (western part of Section 1 in Figure 11) gave reliable results.

To compute the dispersion curve from the noise data, we first calculated the cross-correlation (CC) functions between the vertical components of station pairs. The synchronized recordings of vertical components were first processed using the one-bit normalization and the spectral whitening (Bensen et al., 2007). Then, the CC functions were computed for each station pair. To compute the dispersion curve of the seismic signals emerging from the CC functions, we applied a velocity analysis to them. The method is similar to the constant velocity stack (CVS) analysis (Yilmaz, 1987; Yamanaka et al., 1994), which is very popular in active seismic reflection processing and has already been used on CC results for different Italian areas (Vassallo et al., 2019; Di Giulio et al., 2020). The CC functions were filtered in different frequency bands starting from 0.5 to 20 Hz. For each frequency band, the CC functions were shifted
<figure><img src="image_126.png" /><figcaption>FIGURE 4
Map of f0 interpolation. Red dots represent the noise measurements with the f0 value reported. The fault lines in red come from the ITHACA database. Green diamonds represent the boreholes considered for 1D modeling (Figure 8).</figcaption></figure>

back in time according to the theoretical surface travel times computed for different constant velocities starting from 50 m/s until 900 m/s using a velocity step of 10 m/s. Then, the phase-weighted Stack (PWS, Schimmel and Paulssen, 1997) was computed, and the absolute maximum of the PWS stack was used to estimate the presence of a horizontally aligned phase in the corrected seismic section. For each filter, the maximum stack function provides the velocity of the surface waves at the considered frequency.

The dispersion curves obtained through CC functions and the velocity analysis were finally inverted in order to obtain a preliminary 1D velocity profile of the subsoil. To improve the inversion process, the HVNSR results have also been used as constraints. In particular, the inversion process has been imposed to fit as best as possible to the part of the HVNSR curve around the resonance peak. The computation of the velocity models has been performed using Geopsy software, in particular, the dinker tool (Wathelet 2008): this code adopts a Monte Carlo-like approach, the so-called neighborhood algorithm, to search for subsoil velocity models that fit the data.

The initial parameterization of the subsoil consisted in a four-layered model over a half-space for all the tested sites, for which the shallower layer is assumed to have a linear increasing Vs. However, in order to avoid strict constraints to the process, a wide variability of the parameterization in terms of thickness and Vs of the layers was allowed.

# Results

## HVNSR results

The average spectral ratio curves of all the sites are plotted in Supplementary Appendix S1.

It is worth noting that most of the HVNSR curves showing a single peak, multiple peaks, or flat shapes (CMP2, CMP4, and CMP5) occur in the frequency range of interest (0.2–4 Hz). This observation suggests the presence of a unique strong impedance contrast in the subsoil, likely ascribable to the geological interface between the lacustrine filling and the bedrock of the area (Meso-
<figure><img src="image_127.png" /></figure>

Cenozoic carbonate and Neogene terrigenous successions). For the stations installed in correspondence to the deeper part of the basin (TF18, FN19, and FN20 in Section 1 and Figure 11 and from TF19 to TF26 in Section 2 and Figure 12), a second broad amplification peak after the fundamental one occurs.

Figures 3, 4 show the resonance frequency and the f0 value interpolation maps, respectively. The distribution of f0 follows a double trend: in the northern part of the maps (Sections 1, 2, and 3, location in Figure 4), there is a decrease of f0 moving toward north, whereas in the southern part (Sections 4, 5, and 6, location in Figure 4) the f0 values decrease toward NW. Contour plots of HVNSR curves along the sections are reported in Figure 5.

Sections from 1 to 6 are all W-E oriented and display HVNSR amplitude contouring from north to south of the study area. Assuming that the depth of the impedance contrast is linearly proportional to 1/f0, we can find a

reasonable geologic interpretation of f0. Hence, we can state that sections 1, 2, and 3, representative of the basin part of the Trasacco municipal territory, show a regular decrease in f0, which means the deepening of the impedance contrast toward north and, at the same time, a trend that resembles a structural high feature in the central parts, especially in sections 2 and 3. Moving from the eastern part of Section 2, some HVNSR measurements were performed to follow the general E-W trend of the basin toward the eastern depocenter of the FB (Section 2b—Figure 7). A second broad resonance peak is displayed for the eastern part of this section, which could suggest the presence of a shallower impedance contrast unreachable by all the boreholes available in the area.

Section 4 represents a transition zone where there is a change in the f0 trend compared to sections 1, 2, and 3. In fact, it shows a
<figure><img src="image_128.png" /><figcaption>FIGURE 6
HVNSR contour plots of Section 10 and Section 11 located along seismic line 4 (top) and line 2 (bottom) interpreted in Patruno and Scisciani, (2021). The correspondence between the location of HVNSR contour-plot sections and the seismic line is reported with a red line and A and B extremes of the sections. A sketch of the hypothetical fault is also reported with a dotted red line.</figcaption></figure>

more flat and continuous shape of the impedance contrast without any evidence of the structural high.

Sections 5–6 cross the alluvial fan, and the interpretation of the contouring of HVNSR is less straightforward because there is no information about the total thickness and shear-wave velocity of the alluvial fan deposits. When Section 5 shows, similar to Section 4, a continuous seismic interface, Section 6 reveals a typical valley shape, except for the low-frequency amplification for stations CMP6 and CMP7 located at the western border of the valley which might be due to 2D valley effects. Furthermore, in the eastern part of the section, the transect of ambient noise measurements (Figure 8) highlights a quick variation of f0 (from 2.5 to 1.5 Hz) in a very short distance (200 m) suggesting the

presence of a quite abrupt interruption of the lateral continuity of the main impedance contrast. This behavior could be due to the presence of a normal fault which was reported by the ITHACA database (Figure 1A).

Taking into account that, from available geophysical investigations (multichannel analysis of surface waves) performed on the alluvial fan area (Palombelli, 2014), the average Vs for the shallow depths is around 500 m/s; we can say that the average depth of the Vallelonga valley is between 100 and 150 m. Another interesting observation is that the only three measurements which have a flat HVNSR shape (CMP2, CMP4, and CMP5) are in the central part of Section 6. They then have a similar response in terms of HVNSR as that of the rock
<figure><img src="image_129.png" /><figcaption>**FIGURE 7**
HVNSR section plot of Section 2b. The colored circles indicate the position of each ambient noise measurement along the section and their f0 value. A sketch of the hypothetical fault is reported with a dotted red line.</figcaption></figure>

<figure><img src="image_130.png" /><figcaption>**FIGURE 8**
HVNSR curves of the eastern part of Section 6. The colors reported in the legend correspond to stations located from east (blue) to west (red) of the section, moving from the outcropping bedrock (almost a flat HVNSR average curve) toward the center of the valley (f0 from 2.5 to 1.5 Hz).</figcaption></figure>

sites, suggesting that the presence of alluvial fans can result in deep layers with shear-wave velocities lower than the ones above them (e.g., velocity inversions). In these cases, the capability of the HVNSR technique to retrieve the real resonance frequency of the area is really poor, as also shown by Castellaro and Mulargia, (2009).

Sections 7–8 are both NW-SE oriented and located in the hanging-wall and foot-wall of the Trasacco fault, respectively. They show a similar trend of deepening toward the north but the presence of the structural high in Section 8 is reflected in higher f0 values in the north-western part of the section than Section 7. On the contrary, Section 9 is SW-NE oriented, crossing
perpendicularly the main fault structures as reported for the area in the ITHACA database. In this contour-plot is evident that there is a sudden deepening toward NE of the main seismic impedance contrast at around 2,800 m from the beginning of the section which corresponds to the noise measurement TF16 (see location in Figures 3, 4).

Sections 10 and 11 reported in Figure 6 (location in Figure 4) are drawn along two industrial seismic lines interpreted by Patruno and Scisciani, (2021); being particularly important, we postpone the comparison between the sections and the seismic lines afterward, in the Discussion section. However, both the sections are N-S oriented and focused on areas close to buried fault lines: in particular, Section 10 is N-S oriented, following toward the center of the basin in the eastern border of the Valdelonga valley (from FN16 to FN19 in Figure 3). A sudden decrease of the resonance frequency is reported between FN16 and FN11–FN12. Section 11 highlights a deeper seismic interface than Section 10, but ending, in its northern part, to similar frequency values.

## 1D subsoil models

In order to reconstruct some 1D subsoil velocity models and make inferences on the stratigraphic structure of the investigated area, we consulted public borehole log data, mostly drilled for hydrologic exploration, and few geognostic drillings made available by independent professionals. The basic level of seismic microzonation studies for the Trasacco municipality (www.webms.it) provided shallow geological and geotechnical data mainly located close to the urban area at the southern border of the FB. Other data consist of old stratigraphic logs coming from ISPRA (Italian Institute for Environmental Protection and Research) and “Ente Fucino”, a local managing institution, for old wells drilled during the ‘50s, mainly consisting of stratigraphy of water wells.

We used them to fill the lack of geological information especially where no natural outcrops were available (e.g., in the urbanized area). Unfortunately, none of the available borehole logs reached the geologic bedrock in the Trasacco municipality and the only information that we took into account was the lithology and the thickness of the deposits, just to estimate qualitatively a reasonable average shear-wave velocity of the shallower fillings of the Valdelonga valley. The presence of a large alluvial fan inside the valley represents a big issue for the interpretation of our results. The noise measurements (TF17 and TF18 belonging to Section 1, FT06 and TF13 belonging to Section 2, and TF05 and FT18 on Section 3) that are close to boreholes (green diamonds in Figure 4) with stratigraphic information (Figure 9 left) were selected to make 1D preliminary models of the subsoil with the dinver tool of Geopsy. For the area where the points are located, the availability of the borehole logs,

although the description of the stratigraphic log is not very detailed (some of them were old water wells), allowed us to conclude that the sedimentation conditions are relatively simple. In fact, the shallow layer consists in strata of fine-grained sediments sometimes inter-layered by a gravelly layer (see Figure 9 for details), finally overlying a gray-clay layer. None of these boreholes reaches the depth of the geologic bedrock.

The parameterization of the input subsoil starting model was designed according to the available stratigraphic logs.

The Vs models (Figure 9) obtained by joint-inverting the empirical HVNSR and the dispersion curve obtained through the cross-correlation technique (Figure 10) give an idea of the main impedance contrasts for the area. During the joint inversion, the HVNSR curve was assumed as the ellipticity of the fundamental mode of the Rayleigh wave. The inversion is able to reproduce the main features of the field curves even if in some cases the frequency trough present in the HVNSR curve is not found (e.g., for TF05 and FN18). Due to the simplified approximations of our inversion on the target curves, the 1D velocity results can be considered rough models with an increasing uncertainty at larger depths. However, we observe two main velocity contrasts in the profiles obtained from the inversion (Figure 9): the first one is found at a depth between 40 and 70 m and the second one at a depth between 200 and 520 m. From these preliminary models, we can summarize the results as follows:

1. TF05, FN18, and TF13 sites seem to be located on a structural high with a main impedance contrast located less than 250 m below the ground surface;

2. TF17 and TF18 are in an intermediate setting, with the depth of the main interface being between 280 and 300 m;

3. FT06 is outside the structural high, with a deep interface at 500 m.

## Directional analysis of HVNSR

The directional analysis of HVNSR can highlight the possible 2D or 3D effects in the site response of some noise measurements. We explored specific features of these analyses trying to correlate them with other geological and geophysical information, retrieved from other studies. In particular, we focused on the points closer to the uncertain tectonic elements of the investigated area: the surface fault lines identified from the satellite imagery (ITHACA database) and the buried fault segments hypothesized by the interpretation of the seismic lines (Cavinato et al., 2002; Patacca et al., 2008; Patruno and Scisciani, 2021). The latter was useful to have an idea of the main seismic interfaces present in the area including fault lines and their mutual relations.

Directional HVNSR was computed for sites located along sections 1, 2, and 3 (Figure 4), and the results are plotted in Figures 11–13.
<figure><img src="image_131.png" /></figure>
<figure><img src="image_132.png" /><figcaption>FIGURE 10
Cross-correlation of the ambient noise of synchronous HVNSR measurements in the northern part of the Trasacco area (Section 1).</figcaption></figure>

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

For Section 3 (results in Figure 11, Figure 12, Figure 13), we observe a marked azimuthal dependence of the resonance frequency; moving from 80 to 170°, we notice a change in the f0 value from 0.85 to 1.2 Hz, respectively, for site TF08. A similar behavior is observed for stations FT06 and FT07 for Section 2 (results in Figure 12); moreover, in this case, the stations between TF12 and TF15 show a second peak close to the fundamental one not justifiable with abrupt changes of the stratigraphy.

For Section 1 (results in Figure 11, Figure 12, Figure 13), the directional analysis pointed out a clear azimuthal variation of HVNSR, for FN05 and TB04 stations, close to the Trasacco fault branch. In terms of HVNSR amplitudes, we observe at all the mentioned sites, minimum values around f0 for the azimuth parallel to the strike of the nearby fault. This is a feature also observed by Matsushima et al. (2014). Concerning the frequency of the HVNSR peaks, the results of Section 3 seem to be the clearer to interpret. In particular, the TF05 site, which is located on the structural high, shows a directional HVNSR revealed by the shift of the f0 peak from 1 Hz and the azimuth perpendicular to the strike of the fault (NW-SE direction), to the f0 peak of 1.3 Hz and the azimuth parallel to the strike of the fault. The TF08 results, conversely, highlight the opposite condition: the f0 peak value moving from values of 0.8 to 1 Hz for the azimuth parallel and perpendicular to the fault strike, respectively.

## Discussion and conclusion

The detection of subsoil geological structures is a common goal for many geoscience disciplines. It is in general faced by integrating results coming from different geophysical techniques. The HVNSR technique is particularly effective for the estimations of the resonance frequencies of a given site and, by using independent information, for the estimations of the thickness or the average shear-wave velocity of the surface layers below the measurement point. It works well in

simple geological settings, such as the sedimentary contexts, with horizontal and parallel low-velocity layers over a stiff seismic bedrock. The FB, of Quaternary age, only partially fulfills this geological simplicity, being interested by several faults, not all exactly known, that interrupt the homogeneity of the geological strata. Also, the presence of alluvial fans in large parts of the basin is a big issue for the right interpretation of HVNSR curves. In this work, we collected all the previous geological and geophysical information about the Trasacco area. This information not only drove the choice of the dense single-station noise survey but also helped in the interpretation of the HVNSR results. The starting motivation of this survey was to have an idea of the local structure of the continental basin but also to estimate the depth of the calcareous or siliciclastic bedrock. However, we were aware that the FB has a quite complex geologic setting: it has plenty of reverse and normal faults which, during their past activation, caused a relevant displacement of the geological layers and the interruption of their lateral continuity.

Nevertheless, we believe that the results of our noise survey added useful information for the subsoil reconstruction of the FB, shedding some lights on this complex topic, especially for the presence of buried structures in the basin area. The interpolation map of f0 allowed us to observe a double trend of the main impedance contrast of the subsoil. When the northern part of the study area shows a lowering of f0 values toward north (center of the basin), the f0 isolines (Figure 4) in the southern part show a slight change in the trend of the main impedance contrast which now dips toward NW; this variation can be explained by the coexistence of the Trasacco alluvial fan and the Vallelonga valley which could both locally influence the seismic response. Regarding the identity of the impedance contrast responsible for the fundamental resonance peak, we believe that it does not correspond to the top of calcareous rocks but to siliciclastic deposits, supported also by the high thickness values of flysch deposits reconstructed for the study area by Mancinelli et al. (2021).
<figure><img src="image_134.png" /></figure>

Organizing the measurements in regular linear geometries has allowed the comparisons with the interpreted seismic lines already available in the FB. Despite a partial disagreement between the authors (Cavinato et al., 2002; Patacca et al., 2008; Patruno and Scisciani 2021) on the interpretation of some parts of the commercial seismic lines available in FB, the HVNSR curves interpolated along the sections (Figure 6) revealed similar results to what were provided by seismic reflection lines for depths compatible with the resonance frequency of the area.

The extension of the surveys to the scale of the municipality gave the chance to promptly reveal the presence of abrupt subsoil lateral heterogeneities as the presence of fault lines. This is easily proved following the f0 isolines from the map of f0 interpolation (Figure 4), and checking the trend of the HVNSR curves along the transects. Moreover, the strong azimuthal variation of HVNSR for the sites close to the tectonic elements supports their real presence where they are expected to be, as highlighted by Matsushima et al. (2014). A common feature for the HVNSR directional results is that, for sites close to tectonic structures, the
<figure><img src="image_135.png" /><figcaption>FIGURE 12
Directional HVNSR of Section 2.</figcaption></figure>
<figure><img src="image_136.png" /></figure>

spectral ratio curves display minimum amplitude around the f0 for the azimuth almost parallel to the strike of the fault. This seems to prove that the presence of lateral subsoil heterogeneities strongly influences the noise wavefield. This statement becomes more difficult to justify for TF12-TF13-TF14-TF15 sites (Section 2), where the amplitude values change for azimuth for the two HVNSR peaks which seem to be uncorrelated to the strike of the fault, maybe for a combined effect with other subsoil structures unknown so far. Concerning the azimuthal dependence of f0 peak values, while Section 1 does not show important f0 variations close to the Trasacco fault, and Section 2 is complex to interpret because of the disagreement in the structural reconstruction from commercial seismic lines which raise uncertainties in the existence and location of possible fault segments, Section 3 (HVNSR contour-plot in Figure 5 and

directional HVNSR in Figure 13) leads itself to comment on that topic because its subsoil structure can be reasonably interpreted as a 2D structure. As discussed in the Results section, directional HVNSR curves for TF05 and TF08 sites (see location in Figure 3), among others, show a slight variation of the peak frequency related to the azimuth of computation. TF05, which we assume to be located on a structural high respect to the TF08 site (see Figure 9 for the numerical 1D simulation and comments in the related section), has an HVNSR curve with a lower f0 value (1 Hz) for the azimuth perpendicular to the fault strike and, conversely reaching a maximum of 1.3 Hz for the azimuth parallel to the fault. TF08 shows azimuthal dependence but with an opposite behavior in terms of f0 values changing from 0.8 to 1 Hz from azimuth parallel to perpendicular, respectively. We
propose the following: these sites can be affected by the proximity
of a lateral variation in terms of frequency of the resonant layer;
this means that the trend of subsoil geometries for the area drives
the energy of the noise wavefield. To demonstrate directional
effects on HVNSR measurements, Matsushima et al. (2014)
showed that performing numerical modeling in these cases
can help understand the contribution of 1D versus 2D and
3D structures on the wavefield and, therefore, on the HVNSR
curves.

Many of the noise measurements were recorded simultaneously, and gave us the opportunity to use these data to estimate the dispersion characteristics of the noise wavefield in the 1.5–10 Hz frequency range. Finally, the dispersion curve obtained using both the cross-correlation and constant velocity stack (CVS) analyses was inverted to retrieve a reliable velocity profile for the shallower layers of the western edge of Section 1. Moreover, the dispersion curve (Figure 9, right) was used, together with the HVNSR curves, as a constraint for the calculation of 1D velocity models (Figure 9) for some test sites located close to boreholes where the stratigraphy was pretty well known. Of course, this inversion process assumed that the dispersion curve obtained for the western part of Section 1 is representative of the shallow velocities of the subsoil for the entire area. We believe that this assumption is fairly realistic because it is supported by the stratigraphic logs available for the area which reflect the general homogeneity of the lithologies representative of the subsoil. However, we assigned variable depth and thickness for all the layers included the inversion process as free parameters to avoid over-constraining of the initial model. The results of these inversions in terms of velocity profiles show a variability of the depth of the main geological interface.

Finally, taking into account all the limitations of the HVNSR
technique, our study has demonstrated that if this method is used
as a preliminary investigation tool to identify and locate the
hypothetical hidden heterogeneity, it can give some important
contributions for seismic risk assessment studies, with the
advantage of being an easy application and cheap when
compared with other geophysical techniques. Remarkably, the
HVNSR technique can give significant advancements on the
geometrical reconstruction and the geological evolution of an
area similar to the FB, for which extensional quaternary tectonics
played an important role for its general setting. The results of this
study highlight the presence of a structural high hidden by
Quaternary deposits in the northern part of the Trasacco
territory within the basin area. Its shape was very well
constrained by the dense mesh of HVNSR measures. Another
important finding is that the directional analysis for the sites of
Section 2 located over this structure (from TF12 to TF15),
combined with the interpretation of seismic lines and a good
stratigraphic knowledge for the area, suggests that the double
peak in HVNSR curves is not related to a stratigraphic effect but
probably connected to the lateral heterogeneity of the subsoil.

Nevertheless, the FB has been the object of several
geophysical studies in the past, and apart from some
uncertainties in the precise localization, most of the hidden
structures are known.

The aim of our noise measurement survey was to verify if,
even in a complex context as the FB, the HVNSR technique was
able to give useful information.

Our study in the Trasacco area has demonstrated that most
of the HVNSR curves show a unique frequency peak (as in the 1D
conditions), but the more the measurement point is close to
hidden tectonic elements, the more complex the HVNSR curves
become. This major complexity also affects the azimuthal
variation of the HVNSR frequency peaks. The loss of
continuity of the resonance peak as well as the increase of
azimuthal variations are directly related to the hidden faults.
Therefore, the HVNSR technique can be used as a preliminary
investigation tool to reveal the presence of these hidden tectonic
elements.

Therefore, the HVNSR technique has a double potential:
assessing the thickness and then the geometry of the sedimentary
layers, and being an indicator of potential hidden structures.

Data availability statement

The original contributions presented in the study are
included in the article/Supplementary Material; data presented
in this study are available on request from the corresponding
author. Further inquiries can be directed to the corresponding
author.

Author contributions

DF: substantial contributions to the conception or design of the work, drafting the work or revising it critically for important intellectual content, acquisition, and analysis or interpretation of data for the work. FC: substantial contributions to the conception or design of the work, drafting the work or revising it critically for important intellectual content, acquisition, and analysis or interpretation of data for the work. GDG: drafting the work or revising it critically for important intellectual content, acquisition and analysis or interpretation of data for the work. MV: acquisition and analysis or interpretation of data for the work. GM: drafting the work or revising it critically for important intellectual content.

Funding

This work has been funded by MZS Avezzano livello 3
project. The authors DF, FC, GDG, and MV during part of
this research work held temporary positions at INGV supported
by FIRB Abruzzo - UR07 (n. RBAP10ZC8K).

Acknowledgments

The authors thank the Laboratorio Effetti di Sito (https://www.ingv.it/it/monitoraggio-e-infrastrutture-per-la-ricerca/laboratori/laboratorio-effetti-di-sito) for providing seismic equipment used in the field measurements. The authors are grateful for the insightful comments and suggestions provided by Francesco Maesano about the study area and Alessia Mercuri for providing useful suggestions with Geopsy software tools and her comments in different stages of this research project. Useful contribution to improve the paper was given by the reviewers PB and SM.

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.

References

Accardi, G. (1975). Nuovi affioramenti di lacustre sollevato a terrazzi al bordo settentrionale del fucino (abruzzi). *Boll. Soc. Nat. Napoli* 84, 1

Amoroso, S., Bernardini, F., Blumetti, A. M., Civico, R., Doglioni, C., Galadini, F., et al. (2016). Quaternary geology and paleoseismology in the fucino and L'aquila basins. *Geol. Field* 8 (1,2), 1–88. doi:10.3301/GFT.2016.02

Bard, P.-Y., and Bouchon, M. (1980a). The seismic response of sediment-filled valleys. Part 1. the case of incident SH waves. *Bull. Seismol. Soc. Am.* 70, 1263–1286. doi:10.1785/bssa0700041263

Bard, P.-Y., and Bouchon, M. (1980b). The seismic response of sediment-filled valleys. part 2. the case of incident P and SV waves. *Bull. Seismol. Soc. Am.* 70 (5), 1921–1941. doi:10.1785/bssa0700051921

Bard, P.-Y., and Bouchon, M. (1985). The two-dimensional resonance of sediment-
filled valleys. *Bull. Seismol. Soc. Am.* 75, 519–541. doi:10.1785/bssa0750020519

Bensen, G. D., Ritzwoller, M. H., Barmin, M. P., Levshin, A. L., Lin, F., Moschetti, M. P., et al. (2007). Processing seismic ambient noise data to obtain reliable broadband surface wave dispersion measurements. *Geophys. J. Int.* 169 (3), 1239–1260. doi:10.1111/j.1365-246x.2007.03374.x

Blumetti, A. M., Dramis, F., and Michetti, A. M. (1993). Fault-generated mountain fronts in the central apennines (Central Italy): geomorphological features and seismotectonic implications. *Earth Surf. Process. Landf.* 18, 203–223. doi:10.1002/esp.3290180304

Bosi, C., Galadini, F., and Messina, P. (1995). Stratigrafia Plio-Pleistocenica della conca del fucino. *Il Quat.* 8 (1), 83

Cara, F., Di Giulio, G., Cavinato, G. P., Famiani, D., and Milana, G. (2011). Seismic characterization and monitoring of fucino basin (central Italy). *Bull. Earthq. Eng.* 9, 1961–1985. doi:10.1007/s10518-011-9282-2

Castellaro, S., and Mulargia, F. (2009). The effect of velocity inversions on H/V. *Pure appl. Geophys.* 166, 567–592. doi:10.1007/s00024-009-0474-5

Cavinato, G. P., Carusi, C., Dall'Asta, M., Miccadei, E., and Piacentini, T. (2002). Sedimentary and tectonic evolution of plio-pleistocene alluvial and lacustrine deposits of Fucino Basin (central Italy). *Sediment. Geol.* 148, 29–59. doi:10.1016/
s0037-0738(01)00209-3

Cella, F., Nappi, R., Paoletti, V., and Florio, G. (2021). Basement mapping of the fucino basin in central italy by ITRESC modeling of gravity data. *Geosci. (Basel)*. 11, 398. doi:10.3390/geosciences11100398

Di Giulio, G., Ercoli, M., Vassallo, M., and Porreca, M. (2020). Investigation of the Norcia basin (Central Italy) through ambient vibration measurements

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.

Supplementary material

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

SUPPLEMENTARY MATERIAL S1

Average HVNSR for all the single station ambient noise measurements
performed in the study area.

SUPPLEMENTARY TABLE S1

Main parameters of the HVNSR recordings: Site Name, Lat, Lon, f0n.

and geological surveys. *Eng. Geol.* 267, 105501. doi:10.1016/j.enggeo.2020.105501

Doglioni, C., Mongelli, F., and Pialli, G. P. (1998). Boudinage of the alpine belt in the apenninic back-arc. *Mem. Soc. Geol. It.* 52, 457

Famiani, D., Brunori, C. A., Pizzimenti, L., Cara, F., Caciagli, M., Melelli, L., et al. (2020). Geophysical reconstruction of buried geological features and site effects estimation of the middle valle umbra basin (central Italy). *Eng. Geol.* 269, 105543. ISSN 0013-7952. doi:10.1016/j.enggeo.2020.105543

Frezzotti, M., and Giraudi, C. (1992). Evoluzione geologica tardo-pleistocenica ed olocenica del conoide complesso della Valle Majelama (Massiccio del Velino - Abruzzo). *Il Quat.* 5 (1), 33

Galadini, F., Galli, P., and Giraudi, C. (1997). Geological investigations of Italian earthquakes: new paleoseismological data from the fucino plain (central Italy). *J. Geodyn.* 24, 87–103. doi:10.1016/s0264-3707(96)00034-8

Galadini, F., Galli, P., Giraudi, C., and Molin, D. (1995). Il terremoto del 1915 e la sismicità della piana del Fucino (Italia Centrale). *Boll. Soc. Geol. It.* 114, 635–663. doi:10.12910/EA12015-075

Galadini, F., and Galli, P. (1996). Paleoseismology related to deformed archaeological remains in the fucino plain. implications for subrecent seismicity in Central Italy. *Ann. Geofis.* 39, 925–940. doi:10.4401/ag-4025

Galadini, F., and Galli, P. (1999). The Holocene paleoearthquakes on the 1915 Avezzano earthquake faults (central Italy): implications for active tectonics in the central apennines. *Tectonophysics* 308, 143–170. doi:10.1016/s0040-1951(99)00091-8

Galli, P., Messina, P., Giaccio, B., Peronace, E., and Quadrio, B. (2012). Early Pleistocene to late Holocene activity of the magnola fault (fucino fault system, central Italy). *Boll. Di Geofis. Teor. Ed. Appl.* 53 (4), 435–458. doi:10.4430/bgta0054

Giraudi, C. (1988). Evoluzione geologica della Piana del Fucino (Abruzzo) negli ultimi 30.000 anni. *Il Quat.* 1 (2), 131.

Giraudi, C. (1986). Faglie ad attività olocenica nella piana del fucino. *Mem. Soc. Geol. It.* 35, 875

ISPRA (1881) Geological survey of Italy. Available at: http://sgi2.isprambiente.it/ithacaweb/Mappatura.aspx

ITHACA Working Group (2019). *ITHACA (ITaly HAzard from CApable faulting), A database of active capable faults of the Italian territory.* Version December 2019. ISPRA Geological Survey of Italy. Available at: http://sgi2.isprambiente.it/ithacaweb/Mappatura.aspx