Dual-proxy estimation of Vs30: the case study of the Marche Region (central Italy)

ABSTRACT

This study focuses on generating a shear-wave velocity averaged within the uppermost 30 m of the ground surface (Vs30) map for the Marche region (central Italy) using two commonly acknowledged proxies: topographic slope and lithological classification. The analysis is based on a comprehensive dataset of geophysical tests from the Italian seismic microzonation dataset, employed as a training set. Through regression analysis, Vs30 values are modelled as a function of lithology and topographic slope, with a random effect accounting for the combination of these variables. The resulting Vs30 raster map illustrates the spatial distribution of shear-wave velocities across the region, offering a representation of the subsurface seismic characteristics essential for various applications, including local seismic hazard assessment, prediction of seismic ground motion parameters, microzonation mapping, real-time shakemap generation, and seismic design of engineering structures.

KEYWORDS:

• Seismic hazard
• shear-wave velocity
• Vs30 map
• Marche region

1. Introduction

The determination of the ground shaking induced by destructive earthquakes is becoming a crucial aspect of seismic hazard assessment. The complexity of the soil conditions at a site might induce seismic wave amplification which influences the local seismic response. Thus, mapping the geophysical characteristics of soil is essential to predict the effects of an earthquake, due to site amplification. The shear-wave velocity averaged within the uppermost 30 m of the ground surface (Vs30) is the most used parameter for the site classification and the evaluation of the soil amplification for the seismic design (Parker et al., 2017; Stewart et al., 2014; Wald & Allen, 2007). In addition, most of the ground motion models (GMMs) use Vs30 as data input representing the local site conditions (Abrahamson & Silva, 2008; Boore & Atkinson, 2008; Campbell & Bozorgnia, 2008; Chiou & Youngs, 2014; Lanzano et al., 2019).

Since Vs30 is considered a key factor for ground motion prediction, it is fundamental to evaluate its spatial distribution to calculate shaking scenarios. Vs30 is a measure directly calculated from the shear-wave velocity (Vs) profiles obtained by geophysical investigations. Borehole methods (i.e. Down-Hole, DH; Cross-Hole, CH), seismic refraction, Multichannel Analysis of Surface Waves (MASW), Spectral Analysis of Surface Waves (SAWS), and REfraction MIcrotremors (REMI) are well-known invasive and non-invasive techniques for estimating Vs (Comina et al., 2020; Kramer, 1996; Wang et al., 2022; Xia et al., 2002). Since detailed geophysical investigations are expensive and time-consuming, it is unrealistic to perform very dense geophysical characterisation of an area, therefore, various proxy-based methods have been proposed for mapping Vs30 over a large region. These methods are typically based on data available for the area, such as the topographic slope and geomorphology-based terrain categories obtained from a digital elevation model (Allen & Wald, 2009; Wald & Allen, 2007; Yong et al., 2012), or geotechnical categories and surface geology (Chiou & Youngs, 2014; Wills & Clahan, 2006; Wills et al., 2015) obtained from geological maps at various scales.

In Italy, systematic seismic microzonation (MS) studies have been carried out after the 2009 L’Aquila earthquake as part of the National Programme for seismic risk prevention and mitigation, an initiative led by the Department of Civil Protection (DPC). These studies provide a considerable amount of data, and therefore, alternative Vs30 maps of the whole national territory are present in the literature (e.g. Forte et al., 2019; Mori et al., 2020). These maps are mainly based on surface geology at 1:100.000 scale, and geological and geomorphological information. Forte et al. (2019) calibrated two maps of seismic shallow soil classification based on soil class (e.g. EC8) and Vs30 from Vs measurements, and a large-scale geological map of Italy. Mori et al. (2020) proposed a Vs30 map for Italy using the 1:100,000 geological map and geomorphological classes as proxies of Vs30 obtained by multiple in-situ geophysical tests from the Italian seismic microzonation dataset.

In this study, we present a Vs30 map (Main Map) tailored for the Marche region, Central Italy, according to its distinct geological and topographic characteristics. The adopted approach differs from previous studies for the implementation of a reclassification of geological formations and a higher-resolution Digital Elevation Model (DEM) to ensure a better representation of the lithology and topography of the region. Moreover, the public microzonation data are integrated for the Marche region with additional data (curtesy of local freelance geologists), and the dataset of the Engineering Strong motion (ESM) database (Luzi et al., 2016, 2020) to ensure a comprehensive representation of the seismic site characteristics. However, the Vs30 measurements are limited by a nonuniform spatial distribution, with a higher amount of data in populated areas, possibly leading to an underrepresentation of Vs30 values for certain lithologies and/or slope classes. The final Main Map is a spatial mapping of the results of regression analysis between Vs30, lithology, and topographic slope.

2. Geological and geomorphological setting of the Marche region

The Marche region belong to the Umbria-Marche Apennines, a fold-and-thrust belt, that resulted from the convergence between the continental of Corsica-Sardinia European margin to the west, and the Adria margin of African origin to the east (e.g. Carmignani & Kligfield, 1990; Malinverno & Ryan, 1986; Mazzoli et al., 2002; Pierantoni et al., 2013). The Umbria-Marche Apennines constitute the external sector of the Northern Apennines and is characterised by an arcuate shape and by the presence of asymmetrical anticlines, that are mostly faulted, verging mainly to the north-east. The study area is characterised by a Mesozoic-Tertiary sedimentary succession represented by the so-called Umbria-Marche Sedimentary Succession. This continuous stratigraphic sequence can be divided into three main geological units: i., Meso-Cenozoic carbonates and hemipelagic limestones; ii., Oligo-Miocene siliciclastic turbiditic deposits; and iii., Plio-Pleistocene peri-adriatic succession (Bigi et al., 1997; Mazzoli et al., 2002; Pierantoni et al., 2013). The Meso-Cenozoic portion of the sedimentary succession crops out extensively in the western part of the Marche region, in the so-called Umbria-Marche and Marche ridges, respectively, which merge to the south to form the Sibillini Mountains (e.g. Calamita & Deiana, 1988). This portion is characterised by a prevalence of carbonate lithotypes (e.g. Calcare Massiccio, Corniola, Calcari Diasprigni, and Maiolica Fms.) and marly lithotypes (e.g. Marne a Fucoidi, Scaglia Bianca, Scaglia Rossa, Scaglia Variegata, and Scaglia Cinerea Fms.) indicating the transition from shallow water to deep marine and hemipelagic environments. Siliciclastic turbiditic marine deposits which constitute the upper portion of the sedimentary succession are distinguishable in the long intra-Apennine synclinal structure represented by the Camerino basin and in the southernmost part of the Marche region. This portion of the sequence involves the Bisciaro, and Schlier Fms., and the Miocene turbidites of the Camerino and Laga Fms. In the eastern Adriatic sector of the Marche region, the Plio-Pleistocene peri-Adriatic succession characterised by the Argille Azzurre Fm., crops out (Pierantoni et al., 2013). Finally, the late Quaternary continental deposits are widely distributed in the foothill zone along the alluvial plains of the Marche region. However, for the purpose of the study, the geological formations are grouped into lithological complexes according with the Vs30 values as in Table 1.

Table 1. List of the lithological complexes identified in the study area.

Name of Lithology	Lithological subclass	Description	Number of observations	Mean Vs30 (m/s)	Standard deviation of Vs30 (σ)
Carbonate	CC	Limestones, dolomitic limestones, micritic limestones	8	492.6	264.7
Marls	MC	Marlstones	2	537	67.9
Calcareous marl	CmC	Marly limestones, calcareous marlstones	55	600.7	212
Argillaceous marl	McC	Marlstones and silty marlstones	63	474.8	174.3
Arenaceous	AC1	Arenaceous-marly clayey turbidites (well-cemented sandstones)	171	471.9	147
	AC2	Arenaceous-silty turbidites (poorly cemented sandstones)	30	373	132.5
Arenaceous flysch	AfC	Resedimented gypsum, gypsum arenites, clays with interbedded sandstones and conglomerates	44	430	107.8
Clays	ClC	Mudstones, silty clays interbedded with sandstones and conglomerates	328	349.8	93.8
Quaternary deposits	QC1	Continental deposits, alluvial and terraced deposits, colluvium, alluvial fan	603	427.2	141.4
	QC2	Shallow water marine deposits, travertine, and soft limestone	22	430.6	101

From a topographic point of view, the landscape of the Marche region is characterised by the presence of higher relief with steeper slopes on the western side, in correspondence of the Sibillini Mountains ridge which is characterised by stiffer lithologies (i.e. calcareous formations). On the Adriatic side, the region is characterised by a predominantly high-hilly territory consisting of soft lithologies and dominated by relatively wide valley floors (Bisci & Dramis, 1991). The main geomorphological features are those connected with fluvial processes. The Marche valleys are characterised by several orders of alluvial terraces of Quaternary age and located at elevation ranging from a few metres (fourth order terrace) up to more than 100 m (first order terrace) above the present-day valley floor.

3. Data collection and analysis

The data collection consists of the geological, geotechnical, and geophysical tests carried out in the framework of microzonation studies acquired by publicly accessible datasets, available at (https://sisma2016data.it/microzonazione/) and (https://qmap-protciv.regione.marche.it/cs/), in-situ investigations, and the shear-wave velocity profiles available in the ESM database (http://esm.mi.ingv.it; Luzi et al., 2016, 2020). The dataset consists in 1806 shear-wave velocity (Vs) profiles. The dataset is validated and processed to eliminate the Vs profiles that are incomplete, without coordinates, with marked velocity inversions and unreliable Vs values (e.g. Vs < 80 m/s). The Vs30 values are calculated and converted into a point map with the following attributes: point ID that is a combination of the Italian Census Codes and the geophysical test ID, coordinate of the survey location, type of geophysical tests, and Vs30 values (Figure 1).

Figure 1. Distribution of the Vs30 values obtained from the geophysical tests carried out in the Marche region as part of the seismic microzonation studies (DPC, 2019). The geophysical tests are categorised into two types: red triangles with the tip pointing downwards for in-hole tests (DH and CH), while black triangles with the tip pointing upwards for surface tests (MASW, REMI, and SAWS), (projected coordinate system UTM zone N33 expressed in metres).

Map of Marche region, Central Italy, showing Vs30 values obtained from geophysical tests conducted for seismic microzonation studies. Different symbols represent the specific survey types.

The geo-lithologic information is obtained from the 1:100.000 geo-lithological map of Italy (courtesy of the Italian Institute for Environmental Protection and Research, ISPRA – Italian Geological Survey; Servizio Geologico d’Italia, 2004). For the sake of simplicity, the original categories described by ISPRA for the whole Italian territory were grouped into a set of homogeneous lithological complexes consistently constraining the lithological framework of the study area. The original categories are reclassified into 10 subclasses, including two ‘Quaternary deposits’ and eight different types of ‘geologic bedrock’ (Jurassic to Pleistocene) encompassing the prevailing lithotypes and similar expected geotechnical behaviour (Table 1). The thematic map representing the lithological complexes outlined in Table 1 is obtained and presented in Figure 2. Approximately 30% of the total area is represented by sandy-clayish lithotypes, whereas marly-arenaceous lithotypes and Quaternary deposits cover 24% and 22% of the total area, respectively. The complexes corresponding to prevailing calcareous lithotypes cover 18% of the total area, while the remaining 4% corresponds to the AfC complex, which represents evaporitic rocks. The lithological complexes with a percentage of area lower than 1% are not considered.

Figure 2. Map of the Marche region showing the identified lithological subclasses of Table 1 (projected coordinate system UTM zone N33 expressed in metres).

Lithological map of the Marche region. A color map showing the lithological subclasses of the Marche region from Table1.

The first step of the procedure consisted of the overlay of the Vs30 data with the lithological units cropping out in the Marche region. Since microzonation studies aim to characterise lithologies that can be affected by seismic amplification, Quaternary deposits and soft rocks have a large number of observations. On the contrary, very few geophysical tests are performed in correspondence of stiff rocks such as limestones (CC class, Table 1), where amplification is not expected.

The mean and standard deviation of Vs30 values derived for each lithology subclass are reported in Table 1. An analysis of variance (ANOVA test) is carried out to verify whether the lithological complexes listed in Table 1 significantly differ from each other in terms of Vs30 values and the relative box-plot is shown in Figure 3. The Anova test revealed a statistically significant difference of Vs30 between the complexes (overall P-value = 3.99e-40), although some lithological subclasses (e.g. McC and AC1; QC1 and QC2) exhibit similar Vs30 values.

Figure 3. Box-plot showing the distributions of Vs30 for the lithological subclasses listed in Table 1. The first quartile, the median value, and the third quartile are reported with the minimum and maximum values of the distribution and the outliers.

Box-plot showing the distribution of Vs30 values for each lithological subclass of Table 1. The box plot also includes quartiles, median, minimum, and maximum values, and outliers.

An additional analysis is performed to verify the correlation between Vs30 and topographic slope.

The topographic data are derived from an ASTER Global Digital Elevation Model (GDEM) with a pixel size of 30 m, from NASA, METI, AIST, Japan Spacesystems and US/Japan ASTER Science Team (2019), which is compatible with the regional scale of the work. The topographic slopes are calculated using the ‘Slope’ tool of the ArcGIS software and are expressed in degrees. This tool employs an algorithm that uses the eight adjacent pixels associated with a specific weighting matrix, as proposed by Horn (1981). The different classes of slope in degrees and their description are shown in the Main Map. In the Marche region, the slope ranges from 0° to more than 35° (i.e. very gentle to very steep).

Given the limited number of observations and the similarity in Vs30 values for certain lithological classes (i.e. McC and AC1; QC1 and QC2), we further investigated whether these classes could be distinguished based on morphological characteristics. The zonal statistics are performed using the ArcGIS Zonal Function to calculate the slope variability inside each lithological subclass. The mean, median, and standard deviation values obtained from the zonal statistics for each lithological subclass are listed in Table 2.

Table 2. Results of the zonal statistics.

Lithological subclass	Slope Mean [°]	Slope Median [°]	Slope Standard Deviation (σ)
CC	31.5	32	14.1
MC	28.3	28	13.6
CmC	24.3	23	13.4
McC	16.6	15	10.2
AC1	20.5	19	11.4
AC2	14.3	13	8.9
AfC	15.4	14	8.6
ClC	13	12	8.1
QC1	8	5	9.4
QC2	15.3	14	9.8

The zonal statistics revealed that the abovementioned subclasses (i.e. McC and AC1; QC1 and QC2), which share similar Vs30, appear to be dissimilar in terms of topographic slope. Therefore, in the subsequent analysis they are kept as separate complexes. On the contrary, we could have joined the subclasses AC2 and ClC, as the distribution of Vs30 values and slope proves to be comparable. However, we kept them separated because they have different lithological composition.

The Vs30 data points are overlain to slope to obtain the final dataset. A linear regression using the mixed-effects model (Abrahamson & Youngs, 1992) is used to obtain an empirical equation to predict Vs30 values as a function of lithology and slope, and using the random effect on the combination lithology and slope.

We obtain the following equation: (1) $\mathrm{Vs}{30}_n={\alpha }_0+S_{\mathrm{Ln}}\ast \text{β}+L_n+ϵ$(1) where α₀ is a constant equal to 457; S_Ln is the slope coefficient for each geo-lithological complex (obtained by the mixed-effect regression); L_n is the coefficient of the geo-lithological complex; β is the topographic slope, in degrees, and n is the number of complexes; ϵ is the error. The total root mean square error is 134.79, while the root mean square error of each lithological subclass is listed in Table 3. Following the approach proposed by Karimzadeh et al. (2019), we also tested the log-normal distribution variable of the Vs30 in relation with the Log10(slope), without obtaining a significant improvement in prediction accuracy. The Main Map is processed with a spatial resolution of 30 m using the ArcGIS Map Algebra tool.

Table 3. Regression coefficients of Equation (1) for the geo-lithological complexes (L_n) and the slope coefficients for each geo-lithological complex (S_Ln), and the root mean square error (RMSE) obtained for each geo-lithological complex.

Lithological subclass	Ln	SLn	Number of observations	RMSE
CC	−16.48	2.40	8	248.8
MC	27.19	−0.11	2	66.1
CmC	132.38	0.007	55	212.1
McC	−2.16	1.38	63	173.1
AC1	11.77	0.15	171	146.9
AC2	−72.22	−0.17	30	132.3
AfC	−29.85	0.34	44	107.7
ClC	−122.17	1.49	328	92.9
QC1	−52.74	3.15	603	138.9
QC2	−26.12	0.25	22	100.9

4. Discussion

The final Main Map that represents the spatial distribution of Vs30 can be subdivided into three main zones that share similar shear-wave velocities. The Vs30 distribution reflects the gradual transition from the Apennine mountains, with the predominance of Jurassic-Cretaceous calcareous geological formations, to a hilly area with the prevalence of terrigenous deposits of Miocene-Pleistocene age. The highest Vs30 values are associated with calcareous lithotypes which constitute the Umbria-Marche ridge. This zone shows average shear-wave velocities of 620–800 m/s. The intermediate Vs30 values in the 440–520 m/s interval are representative of the main arenaceous complex (AC1) in the north-western and southern parts of the region. On the Adriatic side, the silty-clayish complex (ClC) is characterised by lower Vs30 values (≤ 320 m/s). The Quaternary complex represented by the QC1 subclass formed by conglomerates and finer alluvial sediments shows slightly lower Vs30 values (320–400 m/s), while the subclass (QC2) represented by shallow water marine of sands and conglomerates shows a Vs30 range 420–470 m/s.

The analysis of the Vs30 residuals as a function of the lithological complexes and topographic slope is illustrated in Figure 4(a,b). The distribution of the residuals as a function of the lithological subclasses, in Figure 4(a), shows that the complex MC represented by Marls is the class with large positive residuals. Nevertheless, since this complex is characterised by only two observations, and represents less than 1% of the total area, the error could be acceptable. The distribution of the Vs30 residuals as a function of the topographic slope has zero mean and no trend with the slope, as illustrated in Figure 4(b).

Figure 4. (a) Box-plot illustrating the distribution of Vs30 residuals across the lithological subclasses of Table 1; (b) Distribution of the Vs30 residuals across the topographic slope. The red dashed line is the residuals mean.

Figure showing a pair of plots. The first graph is a box-plot that shows the distribution of Vs30 residuals across the lithological subclasses listed in Table 1. The second graph illustrates the distribution of Vs30 residuals across the topographic slope, with a red dashed line indicating the mean of the residuals.

Finally, to investigate if the 30 m resolution of the Main Map is significant in terms of Vs30 variability, the focal statistic using the ArcGIS Focal Function tool is performed on the final map. This function calculates a statistic for the input cells within a neighbourhood. The Vs30 variability of the Main Map raster is calculated in terms of standard deviation (σ) using different neighbourhood shapes (i.e. 10 × 10, 5 × 5, and 3 × 3). Small neighbourhoods (e.g. 3 × 3 or 5 × 5) share similar variability, in terms of standard deviation (with maximum observed standard deviation of 172). Increasing the size of the neighbourhood, the focal statistic reveals a reduction in variability in terms of standard deviation (with maximum observed variability equal to 158). Therefore, a Vs30 resampling from 30 m to 150 m does not alter the ground motion variability, with the advantage of a reduction of size of the map. A resampling of the Vs30 map to 300 m could reduce the variability of ground motion.

5. Conclusions

In this work, the variability of the local site conditions that might induce a seismic wave amplification is explored in terms of the spatial distribution of the average shear-wave velocity in the uppermost 30 m, a proxy commonly used for soil amplification in ground motion models. A set of shear-wave velocity profiles is used as a training dataset to calibrate an empirical relation to predict Vs30 based on lithology and topographic slope. The Vs30 Main Map is obtained by linear regression, with the random effect on the lithology-slope combination. The distribution of Vs30 reflects the predominant lithotypes characteristic of the Marche region, as depicted in Figure 3, with the topographic slope contributing to the variation within each lithological complex. The map is distributed with 30 m resolution and a resampling to 150 m does not alter the original Vs30 variability. However, it is worth noting that since Vs30 are punctual values, the presented methodology, aimed at the spatialization of Vs30 covering unsampled areas is likely affected by possible inaccuracy due to other factors (e.g. the geological basemap scale, DEM inaccuracy, other proxies not considered in this study). Therefore, the final Vs30 map can be employed in the framework of modelling or even other studies at a regional scale (e.g. shakemaps), but it is not designed for site-local scale purposes, possibly leading to an inaccurate estimation of the site response if employed at this site-scale.

Software

The vector and raster data (geology, DEM, and geophysical investigations) used to develop the basemaps and the final Main Map were processed using ESRI ArcGIS 10.8 ® (Academic Licenses provided by University of Camerino, Unicam). The statistical analyses were performed using MATLAB (R2021a) (Academic Licenses provided by University of Camerino, Unicam) and the Spatial Analyst tools of ArcGIS 10.8.

Acknowledgements

This work was supported by the FAR Unicam project ‘Novel Approach for Seismic Hazard Analysis—NoHard’, coordinator prof. Emanuele Tondi. An additional thanks to Dr Giovanni Lanzano and the National Institute of Geophysics and Volcanology group of Milan Division for their assistance and support in this project.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The ASTER Global Digital Elevation Model (30 m cell-size) is free downloadable from NASA website (https://search.earthdata.nasa.gov/). The geo-lithological map of Italy (scale 1:100.000), courtesy of the Italian Institute for Environmental Protection and Research, ISPRA – Italian Geological Survey; Servizio Geologico d’Italia, 2004. The training set used for the analysis was collected from various sources, including the Italian Seismic Microzonation website: https://sisma2016data.it/microzonazione/, the web portal of Civil Protection of the Marche region: https://qmap-protciv.regione.marche.it/cs/, the ESM database: http://esm.mi.ingv.it, and additional data from professional studies.

The database of Vs30 adopted in this paper is available upon request to the corresponding author. However, to ensure the accessibility of the research results, the final map is provided in TIFF format.

Additional information

Funding

This work was supported by University of Camerino, Italy (No. STI000104).