Hybrid GMPEs for Region-Specific PSHA in Southern Italy

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.


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 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 0 = XI; Mw = 7.03); the 16 November 1894, earthquake near the northern end of the Messina Straits (I 0 = IX; Mw = 6.12); the 8 September 1905, earthquake in the Gulf of Sant'Eufemia (I 0 = X-XI; Mw = 6.95); and the December 28, 1908, earthquake in the Messina Straits (I 0 = 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 0 = IX; M = 6.07), November 5, 1659, (I 0 = X; Mw = 6.57), and October 13, 1791 (I 0 = IX; Mw 6.14). For Sicily, the largest earthquake that hit the island was the 11 January 1693, eastern Sicily event (I 0 = XI; Mw = 7.32). This earthquake was preceded by a large foreshock on 9 January 1693 (I 0 = VIII-IX; Mw = 6.07). Additional significant earthquakes for Sicily include the 10 March 1786 (I 0 = IX; Mw = 6.14) and the 15 April 1978 (I 0 = 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.

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).

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).  [15], Minelli and Faccenna [16], Polonia et al. [17], Tiberti et al. [18]. Subduction interface from Maesano et al. [19].
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 NNWtrending dextral-or oblique-strike-slip system that was associated with the Mw 6.1, April 15, 1978, Figure 3. Sketch map that illustrates the main tectonic features of the study area. Faults from Monaco and Tortorici [15], Minelli and Faccenna [16], Polonia et al. [17], Tiberti et al. [18]. Subduction interface from Maesano et al. [19].
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 6 of 33 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).

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. 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 NNWdirected 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).

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.

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).

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 Repi ≤ 200 km; and (d) magnitude >3.5. As the magnitude is tied to the radiated seismic energy, the moment magnitude Mw was adopted as the preferred measure. When it was not available from authoritative sources (e.g., for some events with Mw < 4.5), the local magnitude, ML, was used. 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 ( 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. 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, ML), 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 VS,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 VS,30 measurements.

Reference Rock Site Identification
Due to the paucity of waveforms recorded by stations classified as EC8 A by the VS,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.

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.

Reference Rock Site Identification
Due to the paucity of waveforms recorded by stations classified as EC8 A by the VS,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.

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.

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.
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].

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:

1.
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.

2.
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.

3.
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.  [52], CZ14 [51], BI14 [53], and ITA10 [41]. CZ14 was calibrated for rupture distances, and it is plotted against Joyner-Boore distance, considering a scale factor 1:10 for RJB < 20 km. Black line, average of the median prediction from the four GMPEs; dotted gray line, confidence interval calculated considering the upper and lower bounds of the GMPE predictions for the four models; grey circles, single combination of scenario parameters. The percentages of the simulated intensity measures that follow within the confidence interval are also reported.
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 (VS,30 = 310 m/s) and generic rock (GR) (VS,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.  [52], CZ14 [51], BI14 [53], and ITA10 [41]. CZ14 was calibrated for rupture distances, and it is plotted against Joyner-Boore distance, considering a scale factor 1:10 for R JB < 20 km. Black line, average of the median prediction from the four GMPEs; dotted gray line, confidence interval calculated considering the upper and lower bounds of the GMPE predictions for the four models; grey circles, single combination of scenario parameters. The percentages of the simulated intensity measures that follow within the confidence interval are also reported.
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 (VS,30 = 2900 m/s), a GR site (VS,30 = 620 m/s), or a generic soil site (VS,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). 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). 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 (VS,30 = 2900 m/s), a GR site (VS,30 = 620 m/s), or a generic soil site (VS,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). versus epicentral distance) obtained using a set of relations between macroseismic intensity and PGA (GOMAL, [55]; FAEMI10, [56]; MAR992, [57]) and simulated by EXSIM for different soil categories (hard rock, stiff and soft soil).

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. versus epicentral distance) obtained using a set of relations between macroseismic intensity and PGA (GOMAL, [55]; FAEMI10, [56]; MAR992, [57]) and simulated by EXSIM for different soil categories (hard rock, stiff and soft soil).

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 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. The values of the total standard deviations (in decimal logarithm units) of the synthetic data are reported in Table 2 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).  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).

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: , where a is the offset, FD is the distance term:

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: where a is the offset, F D is the distance term: 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 1 and c 2 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: where b 1 and b 2 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 j E 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 i S 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.
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 (Mref = 5.0) and distance (Rref = 1 km) are the same as ITA10. The pseudo-depth h and the coefficients c1 and c2 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. FM is the magnitude term: where b1 and b2 are the calibration coefficients, and Mh is the hinge magnitude. After checking the ground-motion scaling with magnitude, we adopt the hinge magnitudes as the same as for ITA10 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, 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). 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, 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). Table 3. Coefficients of SI17ref. PGA, peak ground acceleration; SA, 5% damping acceleration response spectra.

Intensity
Coefficients of SI17ref 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.  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 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. 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. Figure 16. Standard deviations of SI17ref and ITA10. σ, τ, and φ are the standard deviations of the total, between-event, and within-event residuals, respectively.

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, Figure 16. Standard deviations of SI17ref and ITA10. σ, τ, and φ are the standard deviations of the total, between-event, and within-event residuals, respectively.

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. 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.

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, 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. 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 magnitudedependent 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 Mh is still in the range 6.5-7.0 and the ground motion amplitudes are constant after Mh.
We use the same functional form of Equation [1] without the site effect term (FS), because only the model for RR (FS = 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 shows the median coefficients and the associated uncertainties of the hybrid GMPEs (named SI17hyb). The aleatory variability (σ) is 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 .   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. 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   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. 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 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). Table 5. Percentage increments of the predictions of SI17hyb with respect to SI17ref. 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].   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].

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

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 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 ). 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.

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.

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 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 Tables 7 and 8, 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 Tables 7 and 8), due to the occurrence of several historical earthquakes of high magnitude. 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 TR (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 TR, 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 TR to values comparable to the maximum values obtained in Calabria at higher TR. This behavior might be related to the low b-value of the Gutenberg-Richter distribution (see Tables 7 and 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 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 (Figures 26 and 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 regionspecific 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].   As the final aim of this study is to explore the possibility of using hybrid GMPEs in regionspecific 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]. 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].
Geosciences 2018, 8, x FOR PEER REVIEW 28 of 33 Figure 28. Probability of exceedance (50 years) for the three sites of interest, calculated for the PGA considering the logic tree branches related to the A1MPS16 node (see Figure 22).
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 TR = 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 (TR = 100 years) and SLC (TR = 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 (TR = 200 years) and SLC (TR = 3900 years) are always higher than those related to BRANCH#: 6.

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 RJB 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 regionspecific GMPEs for magnitude-distance pairs poorly sampled by recorded data, we highlight the Figure 28. Probability of exceedance (50 years) for the three sites of interest, calculated for the PGA considering the logic tree branches related to the A1MPS16 node (see Figure 22). 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.

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.