Hybrid GMPEs for Region-Specific PSHA in Southern Italy

Abstract

This paper describes the main findings of the project HYPSTHER (HYbrid ground motion prediction equations for PSha purposes: the study case of souTHERn Italy; supported by the Italian Institute of Geophysics and Volcanology). The goal of the project is to develop a methodological approach to retrieve hybrid Ground Motion Prediction Equations (GMPEs) based on integration of recorded and synthetic data. This methodology was applied to the study area of southern Italy, focusing on the southern Calabria and Sicily regions. The target area was chosen due to the expected high seismic hazard levels, despite the low seismic activity in recent decades. In addition, along the coast of the study area, there are many critical infrastructures, such as chemical plants, refineries, and large ports, which strongly increase the risk of technological accidents induced by earthquakes. Through the synthetic data, the predictions of the hybrid GMPEs have been improved under near-field conditions, with respect to empirical models for moderate to large earthquakes. Attenuation at distances greater than 50 km is instead controlled by the empirical data, because attenuation is faster with distance. The aleatory variability of the hybrid models has strong impact on probabilistic seismic hazard assessment, as it is lower than the sigma of the empirical GMPEs. The use of the hybrid GMPEs specific for the study area can produce remarkable reductions in hazard levels for long-return periods, mainly due to changes in median predictions and reduction of the aleatory variability.

1. Introduction

In recent decades, the calibration of reliable Ground Motion Prediction Equations (GMPEs) in regions of interest have become a critical issue in Probabilistic Seismic Hazard Assessment (PSHA). GMPEs generally describe the amplitude of the ground shaking as a function of the magnitude, distance, and site condition, and they are generally derived from empirical data of past events, when available. The most known and used GMPEs were derived from strong-motion data of worldwide areas where there are dense recording networks, such as California, Japan, Turkey, and Italy [1,2,3].

In areas characterized by high hazard levels where few significant earthquakes have occurred in recent years or where the recording networks are not dense, empirical data can be scarce, and the use of synthetic ground-motion parameters can provide a way to represent ground motion in regions of interest. Indeed, when a site is close to a fault, the rupture processes are predominant, and the finite-source effects dominate the ground motion, such as for directivity, hanging-wall/foot-wall, radiation-pattern and slip distribution. Therefore, the empirical GMPEs can generally not be used to capture such features, because there are few strong-motion recorded data under near-source conditions.

In recent years, some studies have included the use of synthetic data to fill this lack of observations, especially for moderate to large earthquakes under near-field conditions. In particular, the Next Generation Attenuation (NGA)-East project [4] provides a set of new GMPEs for median and standard deviation of ground motions and their associated weights in the logic trees for use in PSHA for the central and eastern North-American regions, based on synthetic data. The calibration of physics-based GMPEs has also been attempted, as these are particularly effective in showing important ground-motion features in near-source regions (Network of European Research Infrastructures for Earthquake Risk Assessment [NERA EU] project, [5]).

A common approach to obtain a hybrid description of ground motion in a given area is the host-to-target method [6]. This calibrates an empirically well-constrained GMPE in a data-rich host region for use in a data-poor target region, based on adjustment factors that are obtained from response–spectral ratios of stochastic simulations in the host and target regions. Another method is the empirical approach that was introduced by Atkinson [7], for which the adjustment factors are determined using spectral ratios of observed motions in the target region adjusted to predictions of an empirical GMPE in the host region. Yenier and Atkinson [8] used this method to develop a GMPE for north-eastern USA through adjustments of the stress and anelastic attenuation, and calibrated the model using the NGA-East database. In Italy, stochastic simulations have been performed in southern Italy for the town of Cosenza, to provide synthetic ground-motion models for earthquakes in a magnitude range and at a distance for which the finite-fault effects are predominant [9]. The analysis of the synthetic ground-motion variability for the implementation of hybrid GMPEs was integrated into the PSHA procedures, as this allows evaluation of the probability of exceedance through statistical analysis of a huge number of possible ruptures that occur on specific composite faults.

In this study, we propose a methodological approach to retrieve ground-motion prediction models that is based on integration between recorded and synthetic data. This methodology is applied to the study case of southern Italy, focusing on the Calabria and Sicily regions. The target area was chosen based on the expected high hazard level, despite the low seismic activity in recent decades. A simple logic tree with only six branches was considered to perform PSHA in the study area: two branches account for epistemic uncertainty in the seismic zonation model, and three branches are related to alternative ground-motion prediction equations. Sensitivity analysis is performed to evaluate the impact on the seismic hazard assessment of the study area, using different GMPEs derived from empirical or hybrid (i.e., recorded and simulated) data. We selected three sites where critical facilities are present: Milazzo (refinery), Priolo Gargallo (chemical plants), and Gioia Tauro (large port). The hazard variability is evaluated considering two different return periods (T_R) for strategic infrastructures, according to the European (EC8, [10]) and Italian (NTC08, [11]) seismic codes, which correspond to damage and collapse-limit states.

2. Study Area

2.1. Historical Seismicity

Southern Calabria and Sicily are areas with long recorded seismic histories, as shown in Figure 1 and

Table 1. Events that have occurred in southern Calabria and Sicily with I₀ ≥ VIII.

Year	Month	Day	Epicentral Area	Latitude [°]	Longitude [°]	I0 Mercalli-Cancani-Sieberg	Moment Magnitude
1169	2	4	Southeastern Sicily	37.215	14.949	X	6.50
1184	5	24	Crati Valley	39.395	16.193	IX	6.75
1542	12	10	Southeastern Sicily	37.215	14.944	X	6.68
1626	4	4	Central Calabrian	38.851	16.456	IX	6.07
1638	3	27	Central Calabrian	39.048	16.289	XI	7.09
1638	6	8	Crotone	39.279	16.812	X	6.76
1659	11	5	Central Calabrian	38.694	16.249	X	6.57
1693	1	9	Southeastern Sicily	37.141	15.035	VIII-IX	6.07
1693	1	11	Southeastern Sicily	37.140	15.013	XI	7.32
1743	2	20	Northern Ionian	39.847	18.774	IX	6.68
1783	2	5	Southern Calabria	38.297	15.970	XI	7.10
1783	2	7	Central Calabrian	38.580	16.201	X-XI	6.74
1783	3	28	Central Calabrian	38.785	16.464	XI	7.03
1786	3	10	Patti Gulf	38.102	15.021	IX	6.14
1791	10	13	Central Calabrian	38.636	16.268	IX	6.14
1818	2	20	Catania	37.603	15.140	IX-X	6.28
1832	3	8	Crotone	39.079	16.919	X	6.65
1836	4	25	Northern Calabria	39.567	16.737	IX	6.18
1854	2	12	Cosenza	39.256	16.295	X	6.34
1870	10	4	Cosenza	39.220	16.331	IX-X	6.24
1894	11	16	Southern Calabria	38.288	15.870	IX	6.12
1905	9	8	Central Calabrian	38.811	16.000	X-XI	6.95
1908	12	28	Messina Straits	38.146	15.687	XI	7.10
1968	1	15	Belice Valley	37.756	12.981	X	6.41
1978	4	15	Patti Gulf	38.385	15.086	VIII	6.03

. Figure 2 shows the seismic history of the study area in terms of epicentral intensity (I₀; hereinafter, Mercalli-Cancani-Sieberg (MCS) scale) and magnitude. Over the last four centuries, the study area was struck by several crustal earthquakes with Mw ≥ 6.0 and I₀ ≥ VIII, as derived from the relevant Italian earthquake catalog (Italian Parametric Earthquake Catalogue, 2015 version [CPTI15], [12]).

The largest events in Calabria include: the main earthquakes of the seismic sequence of February–March 1783, which struck a large region between the southern end of the Gioia Tauro plain and the region south of Catanzaro (I₀ = XI; Mw = 7.03); the 16 November 1894, earthquake near the northern end of the Messina Straits (I₀ = IX; Mw = 6.12); the 8 September 1905, earthquake in the Gulf of Sant’Eufemia (I₀ = X-XI; Mw = 6.95); and the December 28, 1908, earthquake in the Messina Straits (I₀ = XI; Mw = 7.10), which ranks among the most catastrophic events in Italian history (Figure 2a). Other Mw > 6.0 earthquakes affected central Calabria on 4 April 1626 (I₀ = IX; M = 6.07), November 5, 1659, (I₀ = X; Mw = 6.57), and October 13, 1791 (I₀ = IX; Mw 6.14). For Sicily, the largest earthquake that hit the island was the 11 January 1693, eastern Sicily event (I₀ = XI; Mw = 7.32). This earthquake was preceded by a large foreshock on 9 January 1693 (I₀ = VIII-IX; Mw = 6.07). Additional significant earthquakes for Sicily include the 10 March 1786 (I₀ = IX; Mw = 6.14) and the 15 April 1978 (I₀ = VIII; Mw = 6.03), events. Both hit the northern coast of the island between Patti, Milazzo, and the Aeolian Islands, while the northern offshore saw two earthquakes of magnitude >5.0, as the 5 March 1823 (Mw = 5.81) and the 6 September 2002 (Mw = 5.92), events. Over recent decades, several earthquakes of moderate magnitude occurred in the study area, even though the shallow seismicity recorded in the upper crust of the Ionian offshore was limited. In some cases, the moment magnitude was larger than 5.0: 11 March 1978, for Aspromonte (Mw = 5.22), 11 May 1947(Mw = 5.7) for the Ionian coast of Calabria, and 26 May 2001 (Mw = 4.5; [13]), for the Ionian offshore. As for Sicily, Calabria was affected by the 1978 Patti Gulf (Mw = 6.03), the 1979 Sicily Channel (Mw = 5.34), the 1980 Aeolian Islands (Mw = 4.57), and the 2002 Palermo (Mw = 5.92) earthquakes.

2.2. Seismotectonic Framework

Calabria and Sicily are among the most seismically active regions of Italy, and they are characterized by the highest seismic hazard (according to MPS04, the Italian official seismic-hazard map [14]). The Calabrian Arc lies on top of a subduction zone and has a complex geological structure that is largely inherited from the early stages of the collision between the Africa and Eurasia plates (Figure 3).

The subduction of oceanic crust occurs at present only along a relatively small portion of the arc ([19], and references therein), with a rate of a few mm/y (e.g., [20,21]). Contraction has led to the formation of a large accretionary wedge that is currently in the Ionian offshore, which arose by progressive folding of the thick sedimentary cover and includes several thrust faults [16,17,22]. No large earthquake is known for these sources; nevertheless, the size and geometry of the accretionary wedge thrust faults are fully consistent with an expected maximum magnitude of ≥7.0. The Calabrian Arc is characterized by a series of blind, low-angle, normal faults that alternate with strike-slip structures ([18], and references therein). The main fault system of Calabria runs onshore along the Tyrrhenian coast for its entire length and consists of large normal faults that are believed to be responsible for the largest earthquakes in Calabria (and the whole of Italy), such as the M > 7. 1908 and 1905 earthquakes, and part of the 1783 earthquake sequence ([15,23,24,25], and references therein). Strike-slip faults are responsible for the along-strike segmentation of the Calabrian Arc, which divided it into portions that are characterized by a diverse range of geomorphological features and by different uplift rates ([18,26], and references therein). To the southwest, a NW-SE series of faults reflects the presence of a major deep-seated shear zone that bounds the active portion of the Ionian subduction (e.g., [19,27,28]). In the upper plate, the Tindari-Letojanni Line is an important NNW-trending dextral—or oblique—strike-slip system that was associated with the Mw 6.1, April 15, 1978, earthquake. This separates the extension along the Calabria western coast from the contraction for the southern Tyrrhenian coast ([29,30], and references therein).

Tectonic processes in Sicily appear to be more directly affected by the Africa and Eurasia convergence that according to current kinematic plate models acts along a N-NW to NW direction, at 3 mm/y to 8 mm/y (e.g., [31]). Sicily is characterized by thrust faults both in the northern offshore (southern Tyrrhenian) and in the southern part (e.g., [32,33]). In the southern Tyrrhenian, an EW narrow contraction belt runs from the Sicily Channel to the Aeolian Islands, which are about 50 km off the northern coast of Sicily. In the past 30 years, several events with Mw > 5 originated in this area, including the Mw 5.9 Palermo earthquake of 2002 [13]. Also, south-eastern Sicily has large thrust faults. The reverse kinematics of the sources bordering to the north and south of the Catania Plain agrees with geological observations and present-day stress-field and GPS data, which show NNW-directed active shortening between the Hyblean Plateau and Mt. Etna (e.g., [34]). Strike-slip tectonics are seen for the Hyblean foreland. These faults belong to a high-angle, N-S-oriented shear zone that extends both onshore and offshore from the Hyblean foreland area to the frontal part of the Apennines–Maghrebide orogenic wedge. This is a long-lasting structure that has a strong morphological and structural overprint from an older tectonic phase, which indicates right-lateral displacement along a pre-existing crustal inhomogeneity. Although the geological evidence for the present activity is faint, it indicates the present-day left-lateral sense of motion along the shear zone ([35], and references therein).

2.3. Critical Infrastructure

The term critical infrastructure is used to define the essential services to society that represent the backbone of the economy, security, and health. The safety and resiliency of the critical infrastructure is a top priority, and advanced and standardized tools for hazard assessment are necessary, including the use of low-probability events in the application of these tools [36,37]. In the present study, three aspects of the critical infrastructure in the study area were selected (Figure 4) to perform site-specific PSHA: (i) the port of Gioia Tauro; (ii) the refinery of Milazzo; and (iii) the petrochemical complex of Priolo Gargallo.

3. Empirical Dataset

3.1. Dataset Description

An empirical dataset (HYPST_emp_db) has been assembled related to southern Calabria and Sicily that includes to date about 3200 three-component waveforms that were generated by 174 earthquakes that occurred in the period of 1978 to 2016. These were acquired by about 230 recording sites (Figure 5).

The selection criteria are: (a) events and stations included in the geographic area with latitude = 35.6°–40.2° N and longitude = 12°–18° E; (b) waveforms recorded by accelerometric (~1200) and velocimetric (~2000) sensors; (c) epicentral distance R_epi ≤ 200 km; and (d) magnitude >3.5. As the magnitude is tied to the radiated seismic energy, the moment magnitude M_w was adopted as the preferred measure. When it was not available from authoritative sources (e.g., for some events with M_w < 4.5), the local magnitude, M_L, was used.

The bulk of the data (Figure 5b) comes from accelerometric records of the network managed by the Italian Department of Civil Protection (Rete Accelerometrica Nazionale [RAN], code IT, doi:10.7914/SN/IT) and by the Istituto Nazionale di Geofisica e Vulcanologia (INGV; Italian National Seismic Network [INSN], code IV, doi:10.13127/SD/X0FXnH7QfY). Other contributions were derived from the networks MN (Mediterranean Network [MEDNET] project, doi:10.13127/SD/fBBBtDtd6q), managed by the INGV, and YD (Calabria Appennine Tyrrhenian/Subduction Collision Accretion Network, doi:10.7914/SN/YD_2003), managed by the Incorporated Research Institutions for Seismology program for Array Seismic Studies of the Continental Lithosphere, and operated over the time interval of 2003 to 2006. Accelerometric records have been integrated by velocimetric waveforms by networks managed by the INGV or other institutions, which are usually denser than accelerometric waveforms.

All the waveforms were uniformly processed using the strong-motion processing tool of the Engineering Strong Motion database (ESM, http://esm.mi.ingv.it/processing/) [38]. The metadata of both the events and stations were revised: the Italian Seismological Instrumental and Parametric Database catalog (http://iside.rm.ingv.it) was used to update the locations and magnitudes (i.e., Mw, M_L), as well as the moment tensor solutions (i.e., strike, dip, rake) and information on the geometries of the seismic source. Recording sites are classified according to the EC8 seismic code [10] using, when available, the V_S,30 value; otherwise the site category was defined based on the surface geological information [39].

The spatial distribution of the events is denser in the Tyrrhenian Sea and in the eastern part of Sicily, with a cluster of low-magnitude events in proximity to Mt. Etna (Figure 5a). The Tyrrhenian Sea events are generally characterized by deep hypocenters. The recording stations are almost uniformly distributed over the study area, except for the central and western parts of Sicily (Figure 5b).

The dataset is well sampled in terms of recordings for each distance bin (Figure 6a) but is not satisfactory in terms of magnitudes: more than 90% of the recordings are in the magnitude range 3.5 to 5.0 (Figure 6b). Figure 6c shows the focal mechanisms distribution. The prevalent style of faulting within the dataset is normal (NF, ~30%); however, a significant number of records are associated to undefined mechanisms (UN, ~30%), as information on the geometries of the seismic source are not available for low-magnitude events. The EC8 site categories (Figure 6d) are mainly rock and stiff soils, inferred from surface geology (A*, B*, respectively). In particular, few waveforms were recorded by stations identified as EC8 A class by means of V_S,30 measurements.

3.2. Reference Rock Site Identification

Due to the paucity of waveforms recorded by stations classified as EC8 A by the V_S,30 proxy, an alternative approach is used to identify rock sites that are not affected by amplification, hereinafter called reference rock (RR) sites. According to Puglia et al. [40], the method relies on the natural aggregation of the empirical amplification function of the sites, which is obtained by normalizing the acceleration response spectra of the recorded motions according to the EC8-A rock spectrum derived from the most updated ground-motion model for Italy (ITA10, [41]). The mean values of the normalized spectra are thus aggregated based on different amplification trends. Stations that belong to the cluster characterized by the lowest amplification are considered as RR sites. The cluster analysis was performed for 5% damped response spectral acceleration in the period range 0.04 s to 4 s for all the stations included in HYPST_emp_db (about 230). The 22 stations with the lowest misfit between median recorded spectra and ITA10 class A response spectra were selected as RR sites and are reported in Figure 7 (see Puglia et al. [40], for further details). They generally show a lower amplification level with respect to those used to calibrate the ITA10 EC8-A class.

It can be seen that a significant number of sites in Sicily and classified as EC8-A (Figure 8a) cannot be considered as RR sites (Figure 8b). Conversely, several sites located in Calabria classified as EC8-B or EC8-C (Figure 8a) are included in the reference site cluster (Figure 8b).

4. Synthetic Dataset

4.1. Ground-Motion Simulation

HYPST_emp_db is the most complete collection of both weak-motion and strong-motion data for Calabria and Sicily to date, although it is not fully representative of the seismogenic potential of the area under study. We note a lack of information due to either the insular/peninsular geomorphology of the area or the low occurrence rate of moderate-to-strong earthquakes (M > 5.0) in the observation period. In general, the azimuthal coverage of the recording stations is lacking. Several earthquakes with magnitudes between 4.0 and 5.0 are localized offshore with epicentral distances >50 km with respect to the nearest recording station. Most of all, intensity measures of strong earthquakes (M > 6.0) are lacking (Figure 6b).

To fill the gaps in the empirical flat-file, a complementary dataset was assembled (HPST_syn_db) using synthetic ground-motion data. To this end, several region-specific shaking scenarios at the bedrock were performed (for peak ground acceleration [PGA] and pseudospectral acceleration [PSA] at 0.3, 1, 3 s), varying: (i) the magnitude of the simulated events; (ii) the location and kinematic parameters of the individual ruptures; (iii) the stress parameter; (iv) the dip of the simulated fault; and (v) the style of faulting.

We simulate the ground motion using two different simulation methods, as follows:

(1) EXtended fault SIMulation (EXSIM) [42,43], to model a set of Generic Sources (GSs) embedded in the crust structure beneath southern Italy. GSs are not defined by geophysical or geological data; they are constructed for several moment magnitudes (from 5.0 to 7.0, with steps of 0.5), top fault depths, fault mechanisms, and dip angles. The geometry of the simulation sites varies as a function of the moment magnitude and dip, so that the location of the phantom receivers is denser over and in proximity to the Earth surface projection of the fault.

(2) Stochastic-Method SIMulation (SMSIM) [44,45], to model Point-Like Sources (PLSs) characterized by lower magnitudes (3.5 ≤ M < 5.0) and source-to-site distances up to 200 km.

The regional parameters related to the propagation medium properties are fixed for all simulations. A one-dimensional multilayered model representative of the study area was used, in agreement with seismic imaging studies for southern Italy [46,47]. To simulate the ground motion, we considered the quality factor and the geometric spreading valid for Calabria and north-eastern Sicily [48].

The simulated events are normal (NF), reverse (RF), and strike-slip (SS) earthquakes at six magnitudes (5.0, 5.5, 6.0, 6.5, 7.0, 7.5). The dip and top depth of the simulated faults vary as a function of magnitude, while the strike is fixed along the North direction. Point-like sources are modeled for the lower magnitudes (4.5, 5.0, 5.5). For GSs modeling, rupture fronts are used that radially propagate with three different constant velocities. Rupture velocities (RV) are defined in terms of the percentages of the shear-wave velocity (V_S) of the embedding fault medium (70%; 80%; 85%). We consider only one nucleation point that is randomly located over the fault, and random distributions of the slip. For GSs simulations, five different values of the stress parameters are used (50, 100, 150, 200, 250 bars), while for PLSs, the stress-drop values are randomly sampled from a normal distribution with a mean of 150 bars and a standard deviation of 30 bars. We only simulated scenario events for bedrock conditions. To this end, we considered a kappa of 0.035 s, and site amplification factors for hard rock sites (National Earthquake Hazards Reduction Program A, vs. ≥2000 m/s) [4].

4.2. Validation Exercises

The common approach to address the problem of the reliability of simulations is to validate the modeling against records of past earthquakes [49]. A goodness-of-fit between recorded and simulated events is typically given by comparisons of response spectra from simulations with those of selected earthquakes. If the resulting goodness-of-fit is unbiased, simulations can be extended to virtual receivers to complete the ground-motion field of occurred seismic events. However, this is a sound strategy only when a sufficient number of ground-motion recordings is available. In our case, a classical goodness-of-fit is very hard to obtain because strong events are still lacking in the study area. For this reason, we conducted three types of validation exercises:

• Comparisons between synthetics and GMPEs: to verify whether simulated ground motion relatively centered on empirical models, and to determine whether the associated parametric variability can reproduce reasonable ground motion as a whole.
• Comparisons between synthetics and observed data: to evaluate the reliability of simulations against past recorded earthquakes. This exercise is driven by one of the major earthquakes in southern Italy over the past 40 years (1978/04/15, Patti Gulf Earthquake, Mw = 6.0). Some records of the event considered are compared to simulation GSs outcomes characterized by similar magnitude, focal mechanism, epicentral distance, site class, and PGA. The reason here is to determine whether there are simulated waveforms that are well matching to recorded ones in HPST_syn_db.
• Comparisons between synthetics and macroseismic data: to determine whether the simulation method can reproduce the intensity field of past earthquakes. This exercise is driven by the same earthquake of #2 above, for which a consistent number of macroseismic intensity points is available (DBMI15, [50]).

An example of the comparisons between simulations and pre-existing empirical GMPEs among global (CZ14, [51]; BSSA14, [52]), pan-European (BI14, [53]), and regional (ITA10, [41]) GMPEs is shown in Figure 9. All the GMPEs have applicability ranges that are consistent with the simulated ones. The judgement of “very confident” with empirical models defined the scenarios for which at least 80% of the samples (simulated ground motion) falls within the boundary of the propagated standard deviation (Figure 9, gray lines) for all the selected ground-motion parameters.

For the comparison against the observed data, the benchmarks chosen are the records of the 1978/04/15, Patti Gulf earthquake Mw = 6.0 related to the stations IT.GRR (EC8-A*), IT.NAS (EC8-C), and IT.MSS1 (EC8-B*). For the last two stations, amplification factors from Boore and Joyner [54] are accounted for, as specific for generic soil (V_S,30 = 310 m/s) and generic rock (GR) (V_S,30 = 620 m/s) sites. For the three stations, we were able to retrieve at least a pair of synthetics that reproduce the observed ground motion in terms of amplitude, duration, and frequency content of the recorded ground motion. Figure 10 shows that the plots of the simulated waveforms on the S-wave phase of the records provide first-order verification that the synthetics are reliable.

Finally, we simulate the 1978 event considering the source geometry proposed by DISS (https://goo.gl/hhrxgx). We consider the same quality factor and geometric spreading already used to generate scenarios for GSs and PLs, together with a random slip distribution. Each site of the 1978 earthquake macroseismic intensity field (Figure 11a) was previously classified as generic very hard rock site (V_S,30 = 2900 m/s), a GR site (V_S,30 = 620 m/s), or a generic soil site (V_S,30 = 310 m/s), and hence we simulated synthetic acceleration time series by applying amplification factors for the Fourier spectrum that are specific for the soil conditions [12]. After this, the intensity points (MCS) are converted into PGA using three different empirical relationships (GOMAL, [55]; FAEMI10, [56]; MAR992, [57]). We observe that the PGA distribution obtained by EXSIM for the different soil conditions (hard rock, stiff and soft soil) agrees with those provided by GOMAL and MAR992 only for epicentral distances <50 km. Instead, for greater distances, the simulated values are lower (Figure 11b). The best agreement is obtained between synthetic data and PGAs obtained by FAEMI10, within the whole distance range (Figure 11b).

4.3. Dataset Description

The synthetic dataset HYPST_synt_db encompasses more than 180,000 strong-motion data. Figure 12 reports the magnitude–distance distribution of both simulated (GS red circles; PLS, blue circles) and empirical data from HYPST_emp_db (gray circles), together with histograms of the style of faulting, magnitude, and distance for the GS synthetic dataset. About a third of the data are relative to magnitudes >7.0; more than half of the dataset is composed of waveforms of receivers located in the distance range 0 to 25 km. Finally, the dataset is mainly composed of normal and reverse events, rather than strike-slip.

The values of the total standard deviations (in decimal logarithm units) of the synthetic data are reported in

Table 2. Standard deviations of GSs synthetic dataset compared to that of the reference GMPEs for Italy (ITA10, [41]).

Dataset	PGA	SA 0.3s	SA 1s	SA 3s
MW 7.5	0.4073	0.3983	0.3791	0.3603
MW 7.0	0.4871	0.4733	0.4404	0.4123
MW 6.5	0.4718	0.4525	0.4188	0.3869
MW 6.0	0.5167	0.4879	0.4491	0.4148
MW 5.5	0.5505	0.5091	0.4675	0.4466
MW 5.0	0.5812	0.5297	0.4873	0.4725
ITA10	0.3370	0.3630	0.3600	--

and are compared to those proposed by the most recent ground motion for Italy (ITA10, [41]). The standard deviations obtained by the simulations are larger than the empirical ones, especially at lower magnitudes.

To investigate how the parametric variability due to different input simulations contributes to the ground-motion variability, we analyze the synthetic cumulative distribution functions of the PGA computed for all the simulated magnitudes and related to the variation in the style of faulting, depth of the fault, and stress parameter. No significant changes from the overall distributions are observed for the different styles of faulting for all the magnitudes. For the depth of the fault, as the magnitude increases the variability decreases, which shows that larger uncertainties in the source position for lower magnitudes produces larger variability of the ground motion. The most significant contribution is derived from the variability on the kinematic process over the fault (rupture velocity combined with random sample of the nucleation points) and stress parameter (Figure 13).

5. Hybrid Ground Motion Prediction Equations

5.1. Empirical GMPEs

As the first step, we calibrate a set of GMPEs for shallow active earthquakes in southern Italy based on the recorded data collected in HYPST_emp_db. To this end, we select a subset of events with magnitude >4.0 and focal depth <25 km, excluding the inslab events of the Calabrian Arc. We also disregard the volcanic events of Mount Etna and the Aeolian Islands, with depths <5 km, as the ground motion is significantly different from the shallow active one [58,59].

The calibration dataset includes 840 records of 48 events that are recorded by 194 stations. We have adopted the functional form of the ITA10, with some modifications:

$$lo{g}_{10}Y = a+F_D\left(R,M\right)+F_M\left(M\right)+F_S+F_{sof}$$

(1)

where a is the offset, F_D is the distance term:

$$F_D\left(R,M\right) = \left[c_1+c_2\left(M-M_{ref}\right)\right] lo{g}_{10}\left(\frac{\sqrt{R^2+h^2}}{R_{ref}}\right)$$

(2)

where M is the moment magnitude (or the local magnitude without conversion, if the moment magnitude is not available), and R is the Joyner-Boore distance (or the epicentral distance when the fault geometry is unknown). The reference magnitude (M_ref = 5.0) and distance (R_ref = 1 km) are the same as ITA10. The pseudo-depth h and the coefficients c₁ and c₂ are derived from data regression. In contrast to ITA10, the anelastic attenuation in the distance term has been neglected, as the regression coefficients are unstable.

F_M is the magnitude term:

$$F_M\left(M\right) = b_1\left(M-M_h\right)+b_2{\left(M-M_h\right)}^2$$

(3)

where b₁ and b₂ are the calibration coefficients, and M_h is the hinge magnitude. After checking the ground-motion scaling with magnitude, we adopt the hinge magnitudes as the same as for ITA10 (M_h = 6.75). Since the hinge magnitude is outside the range of the empirical dataset, we set F_M to zero for M > M_h. F_sof = f_jE_j (with j = NF, SS, UN) is the style of the faulting term, where E_j are dummy variables and f_j are the regression coefficients. We set the coefficients of class UN to 0 and neglect the reverse fault records (TF) because they are not well represented in the dataset (see Figure 6c). F_s = s_iS_i (with i = RR, GR, ST, SO) is the term that describes the site effects. S_i are dummy variables, and s_i are the coefficients for the four soil classes: (a) RR; (b) GR, i.e., sites of A class, according to EC8 (CEN, 2003), which cannot be considered as RR sites (e.g., affected by amplification); (c) stiff soil (ST), i.e., sites of the EC8-B and EC8-E classes; and (d) soft soil (SO), i.e., sites of the EC8-C and EC8-D classes. In the regression, the coefficient of class RR is set to 0. Figure 14 shows the distribution of the numbers of records for each site class: RR sites represent the 19.4% of the dataset.

The GMPEs were calibrated for the geometric means of the horizontal components of the PGA and three ordinates of the 5% damping acceleration response spectra (0.3, 1, 3 s) that are commonly used for the shake-map calculations [60]. As the spectral ordinates included in the regression are consistent with the waveform-filtering interval, only 832 records were used for the GMPE calibration at T = 3 s. The regressions are performed by applying the random-effects approach [61], with separation of the total residual into between-event and within-event terms [62]. The calibration coefficients of the empirical GMPEs (named as SI17ref) are given in

Table 3. Coefficients of SI17ref. PGA, peak ground acceleration; SA, 5% damping acceleration response spectra.

Intensity	Coefficients of SI17ref
Measures	a	b1	b2	c1	c2	H [km]	fNF	fSS	sGR	sST	sSO	φ	τ	σ
PGA	3.863	0.004	−0.070	−2.039	0.222	11.91	0.036	−0.036	0.479	0.475	0.617	0.107	0.322	0.339
SA 0.3s	4.526	0.375	−0.041	−1.884	0.099	12.24	0.037	−0.037	0.488	0.448	0.607	0.112	0.335	0.353
SA 1s	3.929	0.467	−0.122	−1.548	−0.010	11.16	0.018	−0.018	0.370	0.295	0.484	0.107	0.320	0.337
SA 3s	2.383	0.369	−0.083	−1.301	0.241	10.54	0.058	−0.058	0.335	0.228	0.416	0.135	0.320	0.348

, where σ, τ, and φ are the standard deviations of the total, between-event, and within-event residuals, respectively. Although not shown, the residual trends as a function of distance and magnitude are unbiased (see Supplementary Material).

In Figure 15, the PGA median predictions of SI17ref for soil classes RR and GR are plotted as a function of the distance and are compared to ITA10 EC8-A. We also report the empirical data points used for the regression, and as expected, the GMPEs are not well constrained under near-fault conditions, as very few recordings are available at distances <10 km.

The predictions for the RR class are systematically lower than those obtained for the GR class. The average reduction over all the intensity measures is about 60%, and it is even greater than that observed for Italy, at about 35% [63]. The predictions of ITA10 for EC8-A are similar to those obtained for the GR sites.

Figure 16 shows the comparison between the standard deviations of SI17ref and ITA10 for the intensity measure considered in the GMPE calibrations. The total variabilities of the two models are similar, with a slight reduction of σ at longer periods, by about 6%. The between-event standard deviation τ of SI17ref is significantly lower than for ITA10, with a reduction in the range of 40% to 50%. This is expected, as a local dataset is used for the model calibration (see also [64]), and the event metadata in ESM were revised and improved with respect to those used for ITA10. Finally, the within-event variability of SI17ref is larger than for ITA10.

5.2. Hybrid GMPEs

In this section, we propose a set of hybrid GMPEs calibrated on a mixed dataset of empirical and simulated ground-motion records. We initially compare the different datasets (i.e., HYPST_emp_dtb, HYPST_syn_dtb) we want to use to calibrate the hybrid model. The magnitude–distance distributions of the three datasets (i.e., EXSIM, SMSIM, recorded) are shown in Figure 17.

As the simulations are performed for outcropping rock without topographic amplifications, we considered only the recorded data relative to the RR stations (Figure 17, blue circles), which corresponds to about 180 records. We decided to use the SMSIM data to constrain the near-fault ground motion at magnitudes 4.0 and 4.5 (Figure 17, red circles). For this reason, only the records up to 50 km are included, for 3648 usable records. No selection criteria are applied to the EXSIM dataset, which encompasses more the 187,000 records (Figure 17, green circles). The empirical data are 0.09% of the mixed dataset, while the SMSIM and EXSIM datasets correspond to 1.91% and 98%, respectively. This record imbalance makes it impossible to obtain a hybrid GMPE from the combined data, because of the strong predominance of simulated data with respect to empirical data.

To use a more balanced dataset, we impose that all the empirical data correspond to 15% of the calibration dataset, while the SMSIM and EXSIM contributions are set to 20% and 65%, respectively. We randomly sample the synthetic datasets, respecting the fixed proportions, to obtain a test dataset composed of 1200 records (i.e., 180 empirical, 240 SMSIM, 780 EXSIM).

The basic concept behind the random sampling of synthetic waveforms is to have an almost constant density of data in the magnitude-distance space of interest for the PSHA in this study area; i.e., 0 km to 150 km, and magnitude 4.0 to 7.5. As the empirical data cover only a limited portion of this space (the complete range of distances is only for the magnitude interval 4.0–5.0), most of the data of the hybrid dataset remain synthetic. The record distributions of a generic calibration dataset (named SET1) as a function of magnitude, distance, and style of faulting are shown in Figure 18.

SET1 is still not completely balanced, as some distance–magnitude ranges are not covered by data; e.g., distances >150 km and magnitudes >5.0. On the other hand, SET1 has a consistent amount of data (more than 400) for the largest magnitudes (7.0, 7.5) and for distances <50 km. Similarly, there is a lack of data for the undefined style of faulting (UN) and for magnitudes >5.0. This evidence is related to the characteristics of the simulations: the EXSIM data (M > 5.0) are derived from finite fault simulations and the style of faulting is always known; the SMSIM data are point-source simulations and no focal mechanism is associated to the records (they are flagged as undefined).

Figure 19 shows the distributions of the PGA and the spectral acceleration at 3 s as a function of distance for several classes of magnitudes. SET1 can reproduce the ground-motion attenuation with distance and the scaling with magnitude. Moreover, especially at high frequencies, a magnitude-dependent attenuation can also be noted. The figures of ground-motion parameters as a function of magnitude, for three classes of distance, are available in the electronic supplement: the hinge magnitude M_h is still in the range 6.5–7.0 and the ground motion amplitudes are constant after M_h.

We use the same functional form of Equation [1] without the site effect term (F_S), because only the model for RR (F_S = 0) is derived. In contrast to the empirical GMPEs, we also calibrate the model for reverse faulting. As the coefficients are dependent on the random selection of the simulated data, we consider 50 different replications of the original dataset, with estimation of the median value and the variability of each calibration coefficient. Each replication has the same number of records, to respect the proportions among the three different datasets.

Table 4. Median values and standard deviations of the coefficients of SI17hyb.

Coefficients	Intensity Measures (IMs)
		PGA	SA 0.3 s	SA 1 s	SA 3 s
a	median	3.713	4.031	3.446	2.665
	st. dev.	0.046	0.048	0.043	0.048
b1	median	−0.384	−0.073	0.069	0.363
	st. dev.	0.021	0.027	0.022	0.022
b2	median	−0.056	−0.131	−0.155	−0.117
	st. dev.	0.007	0.008	0.007	0.007
c1	median	−2.242	−1.776	−1.519	−1.38
	st. dev.	0.02	0.016	0.012	0.015
c2	median	0.526	0.276	0.231	0.219
	st. dev.	0.013	0.015	0.014	0.014
H [km]	median	10.33	11.42	10.5	10.23
	st. dev.	0.173	0.13	0.137	0.149
fNF	median	0.052	0.065	0.082	0.098
	st. dev.	0.014	0.014	0.012	0.015
fSS	median	0.105	0.098	0.066	0.055
	st. dev.	0.02	0.021	0.015	0.017
fTF	median	−0.056	−0.055	−0.048	−0.052
	st. dev.	0.016	0.017	0.013	0.013
σ		0.299	0.289	0.278	0.271

shows the median coefficients and the associated uncertainties of the hybrid GMPEs (named SI17hyb). The aleatory variability (σ) is estimated through the error propagation of the standard deviations of the residuals of the different replications. The residuals distributions with distance show unbiased trends (see Supplementary Material Figures).

The uncertainties of the calibration coefficients are significantly lower with respect to those calculated by [41] for ITA10, where 40 bootstrap replications of the starting dataset were considered. The PGA median predictions (±standard deviations) of the SI17hyb (Figure 20, black curves) and SI17ref (Figure 20, gray curves), and the data points of the original dataset (Figure 17) are plotted as a function of distance in Figure 20 for different magnitudes and focal mechanisms. SI17ref predictions reported in Figure 20e,f are calculated outside of the validity range of the GMPEs, as they are calibrated up to magnitude 6.0.

The SI17hyb predictions are significantly different from SI17ref under near-field conditions, as the hybrid GMPEs are better constrained by the simulated data, with larger median values. However, the synthetic and empirical data do not perfectly match from 10 km to 30 km for smaller magnitudes, as the SMSIM amplitudes are greater than the empirical ones. The SMSIM data might need to be adjusted using corrective factors, to tune the ground-motion amplitudes to match the recorded data.

Table 5. Percentage increments of the predictions of SI17hyb with respect to SI17ref.

Intensity Measures	Percentage Increment (%)
	M = 4.0	M = 5.0	M = 6.0	M = 7.0
PGA	394.0	242.3	153.6	100.8
SA T = 0.3 s	25.6	77.5	63.6	13.1
SA T = 1 s	65.2	64.9	42.1	9.8
SA T = 3 s	13.0	48.7	69.4	66.4

reports the percentage increments of the predictions of SI17hyb with respect to SI17ref, under near-fault conditions (R = 0.1 km), for different magnitudes. The results are averaged over the predictions for the different styles of faulting (i.e., normal, strike-slip, undefined focal mechanisms).

Significant increments of the predictions are found for PGA and lower magnitudes, with a maximum value (394%) that approximately corresponds to a factor of five. At longer distances, the predictions of the hybrid GMPEs are more influenced by the empirical data, as they attenuate faster with distance, with respect to the simulated data.

Figure 21 compares the total variabilities of SI17hyb, SI17ref, and ITA10. The total variabilities of the hybrid model are significantly lower than those obtained for the empirical GMPEs, with a reduction of about 17% to 20%. This is because the point-source (SMSIM) simulations cannot reproduce the variability of a real dataset [65].

6. Hazard Assessment

We perform sensitivity analysis to evaluate the impact of the use of different seismotectonics and attenuation models on the seismic hazard assessment of the study area. In particular, the impact at regional scales of hybrid GMPEs is investigated. We also evaluate the hazard variability that corresponds to the Serviceability Limit State for damage control (SLD) and the Ultimate Limit State for collapse prevention (SLC) at three test sites in southern Italy that have critical infrastructures (Figure 4): a large port (Gioia Tauro), refineries (Milazzo), and chemical plants (Priolo Gargallo). Following the Italian building code (Norme Tecniche per le Costruzioni, NTC08 [11]), in case of critical infrastructure located in regions characterized by high seismic hazard, respecting the SLD level should ensure sufficient resistance and stiffness to guarantee immediate re-use, thus reducing the risk of technological accidents induced by natural hazards. Indeed, respecting the SLC should ensure that the construction still has a margin of safety for vertical actions, and a small margin of safety against collapse for horizontal actions. SLD and SLC correspond to specific return periods that depend on the relevance of the infrastructure considered here.

Table 6. Reference period (V_R) associated to the test sites (V_R is the product of nominal life V_N and coefficient of use C_U), exceedance probability for the limit states considered (P_VR SLx), and the related return period T_R SLx.

Test Site	VN	CU	VR	PVR SLD (%)	TR SLD	PVR SLC (%)	TR SLC
Gioia Tauro	≥50	2.0	100	63	100	5	1950
Milazzo	≥100	2.0	200	63	200	5	3900
Priolo Gargallo	≥100	2.0	200	63	200	5	3900

summarizes the reference period (V_R), nominal life (V_N), and coefficient of use (C_U) associated to each test site. The highest V_R is attributed to Milazzo and Priolo Gargallo because chemical plants are industrial installations that can be particularly dangerous for the environment and human health, while ports, such as that of Gioia Tauro, can be treated as an infrastructure that is subjected to overcrowding (i.e., lower V_N).

Logic Tree

A simple logic tree with only six branches is considered to perform PSHA in the study area (Figure 22): two branches account for epistemic uncertainty in the seismic zonation model, and three branches are related to the alternative GMPEs, as previously discussed.

In the present study, the seismic sources are wide seismogenic zones (SZs), in agreement with the previous zonation (ZS4 by [66], ZS9 by [67]) used for Italian seismic hazard maps [68,69]. The first zonation (Figure 22, SHARE; [70]) was proposed and applied for the computation of the European seismic hazard map. It is considered because the SZs and associated seismicity rates are very similar to those used for the present official seismic hazard map of Italy, MPS04 [71]. The second zonation (Figure 22, A1MPS16; [72]) is to be used in the next release of the Italian seismic hazard map (MPS16; https://ingvcps.wordpress.com/category/mps16/). With respect to ZS9 [67], used for MPS04 [69], this SZ considers narrower sources and it is based on new and updated seismological data of the CPTI15 catalog [12]. To compute the values of a and b of the Gutenberg–Richter distribution, the maximum likelihood method is adopted, according to the formulation proposed by [73]. The M_max values are identified based on the maximum observed or estimated earthquakes in each SZ and increasing those estimates by the related standard deviation. The seismic characterization of each SHARE and A1MPS16 zone (Figure 23a,b) is reported in

Table 7. Seismic characterization of each SHARE zone [70]. a and b are the parameters of the Gutenberg–Richter distribution. M_min and M_max are the minimum and maximum moment magnitudes for PSHA. FM: Focal Mechanism; NF: Normal Faulting; RF: Reverse Faulting. h: Reference depth. ASC: Active Shallow Crust; SSC: Stable Shallow Crust.

Id	a	b	Mmin	Mmax	FM	h [km]	Tectonics
ITAS296	4.00	1.00	4.7	6.6	RF	13.2	ASC
ITAS318	3.90	1.00	4.7	8	RF	13.2	ASC
ITAS319	3.90	0.90	4.7	7.7	NF	11.4	ASC
ITAS320	4.00	1.00	4.7	7.7	NF	11.4	ASC
ITAS321	4.30	1.10	4.7	7.7	NF	13.2	ASC
ITAS322	3.90	1.00	4.7	7.7	RF	13.2	ASC
ITAS323	4.30	1.10	4.7	6.6	RF	9.1	ASC
ITAS324	3.90	1.09	4.7	6.6	RF	9.1	ASC
ITAS326	2.80	0.80	4.7	7.7	RF	13.2	ASC
MTAS329	3.90	1.00	4.7	6.6	RF	13.2	ASC
ITAS331	0.18	0.99	4.7	6.6	RF	13.2	ASC
ITAS333	4.20	1.10	4.7	6.6	RF	13.2	ASC
ITAS309	3.60	1.00	4.7	6.5	RF	10	SSC

and

Table 8. Seismic characterization of each A1MPS16 zone [72]. a and b are the parameters of the Gutenberg–Richter distribution. M_min and M_max are the minimum and maximum moment magnitudes for PSHA. FM: Focal Mechanism; NF: Normal Faulting; SS: Strike-slip; RF: Reverse Faulting. h: Reference depth. ASC: Active Shallow Crust; SSC: Stable Shallow Crust.

Id	a	b	Mmin	Mmax	FM	h [km]	Tectonics
617	4.50	1.25	4.5	6.4	NF	17	ASC
602	4.57	1.25	4.5	5.8	SS	15	ASC
604	2.30	0.85	4.5	5.5	SS	25	ASC
606	3.34	0.96	4.5	6.4	RF	10	ASC
610	2.81	0.96	4.5	6.7	RF	15	ASC
611	3.32	0.96	4.5	5.8	RF	5	ASC
613	2.98	0.85	4.5	5.2	SS	10	ASC
605	2.85	0.96	4.5	5.8	RF	15	ASC
603	4.42	1.25	4.5	6.4	SS	10	ASC
600	4.85	1.25	4.5	5.5	SS	17	ASC
608	2.30	0.85	4.5	5.2	SS	10	ASC
679	2.28	0.85	4.5	7.6	RF	12	ASC
219	1.71	0.70	4.5	7.3	NF	7	ASC
200	2.67	0.73	4.5	7.6	NF	10	ASC
204	3.14	0.97	4.5	6.4	NF	9	ASC
205	2.17	0.70	4.5	7.3	NF	9	ASC
206	2.89	0.97	4.5	7	RF	8	ASC
208	1.98	0.70	4.5	7	NF	8	ASC
203	1.89	0.70	4.5	7.3	NF	8	ASC
214	2.39	0.97	4.5	5.8	RF	10	ASC
213	2.67	0.97	4.5	6.4	SS	12	ASC
212	2.52	0.85	4.5	5.5	RF	10	ASC
207	2.52	0.85	4.5	5.5	RF	30	ASC
265	3.60	1.00	4.5	6.5	RF	10	SSC
291	3.99	1.00	4.5	6.6	RF	13.2	ASC
13	3.52	0.87	4.5	5.8	RF	10.2	ASC
15	3.89	1.00	4.5	8	RF	10.2	ASC

, respectively.

For the node relative to the GMPEs, one branch is related to ITA10 [41], calibrated with Italian accelerometric data, while the other two branches are related to the empirical (SI17ref) and hybrid (SI17hyb) region-specific attenuation models calibrated for the study area. All the branches are equally weighted to highlight the influence of the different SZs and GMPEs on the hazard assessment. Figure 24 shows an example of the magnitude scaling of the PGAs considering the three GMPEs used for PSHA, defined for normal and reverse mechanisms.

As SI17ref is not calibrated for reverse faulting, due to the paucity of the empirical data, the median predictions for PSHA (Figure 24b) are relative to those calibrated for an undefined focal mechanism. The reverse faulting predictions of ITA10 are greater than the others, considering that the calibration dataset of ITA10 is represented by little data and is mainly composed of thrust fault events (e.g., M_W = 6.4, 1976 Friuli earthquake). The synthetic dataset, instead, is dominated by reverse faults with steeper dip angles.

The computer program CRISIS 2015 [74] is used here for computation of the expected ground motion in terms of the maps, curves, and uniform hazard response spectra. Figure 25 illustrates the hazard maps at various T_R (100–3900 years) for the PGA and spectral acceleration at 0.3 s and 1 s, respectively. As expected, the hazard level is enhanced, with an increase in the return period T_R, with similar behavior at all periods. In particular, the Hyblean Mountains area (ITAS326 in SHARE, SZ679 in A1MPS16; see Figure 23) goes from very low hazard at lower T_R to values comparable to the maximum values obtained in Calabria at higher T_R. This behavior might be related to the low b-value of the Gutenberg-Richter distribution (see Table 7 and Table 8), due to the occurrence of several historical earthquakes of high magnitude.

The seismic parametrization also influences the difference in PSHA between Milazzo (i.e., SZ291 in A1MPS16, b = 1) and Priolo Gargallo (i.e., SZ679 in A1MPS16, b = 0.85). Figure 26 shows increasing seismic hazard for the latter site, highlighted by the crossing of the hazard curves of Milazzo and Priolo Gargallo (Figure 26). It should be noted that there is a progressive increase in the return period at which the crossing of the two hazard curves occurs, moving from the PGA to higher spectral ordinates. The same inversion in the hazard behavior between Milazzo and Priolo Gargallo can be observed by comparing the uniform hazard response spectra in Figure 27. In general, the highest levels of seismic hazard are for Gioia Tauro (Figure 26 and Figure 27), with a maximum of about 1 g for T_R = 3900 years and spectral acceleration = 0.3 s.

As the final aim of this study is to explore the possibility of using hybrid GMPEs in region-specific PSHA, we analyze the influence of different typologies of attenuation models (among empirical and hybrid ones) by comparing the PGA hazard curves of the three branches of Figure 21 that correspond to the A1MPS16 node (Figure 28). The behavior of the other spectral periods is also shown in Santulin et al. [75].

Focusing our attention on differences between ITA10 (BRANCH#: 4) and the region-specific hybrid model (SI17hyb, BRANCH#: 6), we observe a reduction in the hazard levels of ≤10% for Gioia Tauro, ≤20% for Milazzo, and ≤50% for Priolo Gargallo at T_R = 3900 years. Hazard levels for SI17ref (BRANCH#: 5) are instead remarkably lower. For Gioia Tauro, the differences between BRANCH#: 4 and BRANCH#: 6 are small, both for SLD (T_R = 100 years) and SLC (T_R = 1950). This is because the main contribution to the hazard derives from normal earthquakes 10 km from the site and characterized by moment magnitudes of around 5.0 and 6.0, respectively (see the de-aggregation data in Santulin et al. [75]). In both cases, the median values of ITA10 and SI17hyb are almost identical (Figure 24a), but the standard deviation associated to the hybrid model is a little lower (Figure 21). In contrast to Gioia Tauro, the SZs of Milazzo and Priolo Gargallo are characterized by reverse faulting (Table 8). As a result, because the ITA10 model is systematically higher than SI17hyb (Figure 24b), the BRANCH#: 4 hazard levels for SLD (T_R = 200 years) and SLC (T_R = 3900 years) are always higher than those related to BRANCH#: 6.

7. Conclusions

This study is the first attempt to derive a hybrid region-specific GMPE to be injected into PSHA. This might be a sound strategy for areas such as southern Italy that are characterized by a lack of the empirical data needed to fully capture the peculiarities of the ground motion. To this end, we assembled two different datasets of the ground-motion parameters (PGA, spectral acceleration at 0.3, 1, 3 s) derived from empirical (HYPST_emp_db) and simulated data (HYPST_syn_db). The bulk of the data allows us to calibrate a set of empirical (SI17ref) and hybrid (SI17hyb) GMPEs for rock site conditions, which is useful for southern Calabria and Sicily.

SI17ref predictions are significantly lower than those predicted for Italian rock sites by Bindi et al. [41], by about 60%. The standard deviation is instead slightly reduced with respect to the Italian GMPEs (on average, 6%).

The validity ranges of SI17hyb span from 4.0 to 7.5 for the moment magnitude, and from 0 to 200 km for the R_JB distance. The use of hybrid data for GMPE calibration is efficient, as:

• SI17hyb predictions are larger with respect to SI17ref under near-fault conditions for PGA and lower magnitudes, with a maximum increasing value that corresponds to a factor of 5;
• SI17hyb predictions at distances >50 km are controlled by the empirical data, because they attenuate faster with distance, with respect to the simulated data;
• SI17hyb aleatory variability (σ) is significantly lower than that obtained for the empirical GMPEs (SI17ref), with an average reduction of about 17% to 20%.

The empirical (SI17ref) and hybrid (SI17hyb) prediction models, and the reference GMPEs for Italy (ITA10) were used for hazard calculation with a simple logic tree with only six branches (two branches for the epistemic uncertainty in the zonation model). In testing the performance of region-specific GMPEs for magnitude–distance pairs poorly sampled by recorded data, we highlight the variability of the seismic hazard in terms of probability of exceedance of the ground motion parameter (i.e., PGA) due to the use of different attenuation models.

We also assessed the impact on the seismic design of some of the critical infrastructure, including a port (Gioia Tauro), refineries (Milazzo), and chemical plants (Priolo Gargallo) located along the shoreline of the investigated area. For the region-specific hybrid model, we observe a reduction in the hazard levels of ≤50% for Priolo Gargallo for SLC, with respect to ITA10. Hazard levels for SI17ref are instead lower than SI17hyb, due to the change in the median predictions under near-fault conditions. The highest levels of seismic hazard are observed for the site of Gioia Tauro, with 1 g at spectral acceleration of 0.3 s for SLC.

Further developments in the use of GMPEs for PSHA will include: (i) generation of a new synthetic dataset using regional attenuation parameters (e.g., quality factor and geometrical spreading) inferred by generalized inversion techniques; (ii) a more complex functional form for the attenuation model that includes anelastic attenuation and stress parameter corrective terms; and (iii) a more detailed description of the ground-motion variability, introducing heteroscedastic models for sigma, dependent on magnitude.

This approach is timely for updating the next generation seismic-hazard maps (i.e., the new release of the Italian seismic hazard map MPS16; http://tinyurl.com/jg99xsc), including hybrid GMPEs, in areas where recordings are few and where ground-motion models are not well constrained, especially under near-field conditions. The region-specific hybrid GMPEs can have a strong impact on seismic design and retrofitting of critical infrastructure for both damage and collapse limit states.