Near-source ground motion estimation for assessing the seismic hazard of critical facilities in central Italy

We apply the Probabilistic Seismic Hazard Analysis (PSHA) and compute Physics-Based Simulations (PBS) of ground motion for three dams in the Campotosto area (Central Italy). The dams, which confine an artificial water reservoir feeding hydroelectric power plants, are located in an active seismic zone between the areas that experienced the 2009 L’Aquila and 2016–2017 Central Italy seismic sequences. The probabilistic disaggregation estimated for a return period of 2475 years, corresponding to the collapse limit state for critical facilities, indicates that the most dangerous fault is associated with a maximum magnitude of 6.75 ± 0.25 within a distance of 10 km. This fault is used in PBS to emulate the Maximum Credible Earthquake scenario. To capture the ground motion variability, we input a pseudo-dynamic source model to encompass spatial and temporal variations in the slip, rise time and rupture propagation, heavily affecting the near-source ground motion. Indeed, the ground motion above the rupture volume is mainly influenced by the epistemic uncertainties of rupture nucleation and slip distribution. The computed broadband seismograms are consistent with the near-source shaking recorded during the 2016 MW 6.6 Norcia earthquake and constrain the upper bound of the simulated ground motion at specific sites. Our modelling reinforces the importance of considering vertical ground motion near the source in seismic design. It could reach shaking values comparable to or larger than those of the horizontal components. This approach can be applied in other areas with high seismic hazard to evaluate the seismic safety of existing critical facilities.


Introduction
Seismic hazard quantifies the ground motion expected at a particular site and provides the seismic parameters needed to mitigate earthquake damage. Detailed hazard estimates are routinely performed for critical facilities such as nuclear power plants, dams, gas storage, and pipelines using probabilistic or deterministic approaches or integrating both (e.g., Renault 2014;Dello Russo et al. 2017;Tamaro et al. 2018;Moratto et al. 2021).
The Probabilistic Seismic Hazard Analysis (PSHA) (Cornell 1968) is currently the standard approach required by the Italian seismic code and is widely used for spatial and punctual planning strategies. The PSHA simulates the ground motion by systematically integrating all possible seismic source models obtained from joint seismological and geological information for a given area. It estimates the probability that various levels of earthquake-caused ground motions will be exceeded at a given location in a given time. Ground Motion Models (GMM) are widely used in seismic hazard, but they are not well constrained in the near-source (e.g. Mai 2009;Sgobba et al. 2021) because the representativeness of the data is not statistically significant at short distances from the sources, particularly for normal faults (e.g. Valentová et al. 2021). The finite fault effects are essential to simulate ground motion in near-source for seismic hazard assessment and infrastructure design; recent GMM (e.g. Spudich et al. 2014) account for finite-fault effects (e.g. focal mechanism, directivity), but the ergodic assumption, the poor integration of spatial correlation, and the lack of data in the near-source severely limit the accuracy of near-source ground motion prediction (e.g. Pischiutta et al. 2021;Xin and Zhang 2021). Oversimplified models based on isotropic assumptions may in fact be gross oversimplification in near-source domain and may be unreliable for assessing the seismic risk of critical facilities (Schiappapietra and Smerzini 2021). Furthermore, modern earthquake engineering requires realistic simulations of near-source broadband ground motions, also because the time series are needed as input for more realistic dynamic analyses of the critical facilities (Mai et al. 2010).
Therefore, it is preferable to adopt Physics-Based Simulations (PBS) of ground motion, which is defined as the seismic ground motion computed using numerical methods and models that take into account the physics of the earthquake source and the resulting propagation of seismic waves (Taborda and Roten 2015). Validation results show that the records can be satisfactorily reproduced by PBS, especially at low frequencies where both the source and the material models are better constrained. This method can be integrated at higher frequencies by the stochastic approach, which can reliably reproduce the physics of earthquakes, even if it does not necessarily solve the mathematical problem describing the physics of the source and wave propagation. PBS is a fundamental tool for assessing the potential impact of seismic events on critical facilities near-source, because it can enrich information about the shaking near the source, where the source effects may dominate, and strong motion records may be lacking. Typically PBS evaluates the strong ground motion at selected sites for a given earthquake scenario, usually the Maximum Credible Earthquake (MCE); the disaggregation of PSHA is an essential tool to identify the earthquake scenario for PBS computation. Therefore, PSHA can be applied for evaluations in the broader area, while PBS can be adopted for the detailed design of critical structures based on the controlling earthquake (Krinitzsky 2003).
This study estimates the seismic input for the hazard assessment in an active seismic zone of Campotosto (Central Italy), an area hosting the second-largest artificial water reservoir in Europe bordered by three dams and hydroelectric power plants (Fig. 1). Strong earthquakes struck the zone in 1619 and 1672, and other related earthquakes occurred in 1639 and 1703 (Rovida et al. 2020). The area is located between the 2009 L'Aquila sequence (M W 6.3), related to the Paganica fault (e.g. Chiarabba et al. 2009) and the 2016-2017 Central Italy seismic sequence (e.g. Chiaraluce et al. 2017;Vuan et al. 2017;2018), dominated by three mainshocks (M W 6.2, on August 24, 2016: M W 6.1, on October 1 3 26, 2016; M W 6.6, on October 30, 2016) and activating the Mt. Vettore-Mt. Bove fault system (e.g. Falcucci et al. 2018). These sequences occurred on NW-SE striking faults (ITHACA, 2019), SW dipping and aligned to the system of normal faults along the Apennine chain (Falcucci et al. 2018).
Paleoseismology evidenced significant surface-rupturing events capable of triggering an M > 6.5 earthquake (Galadini and Galli 2003) in the Campotosto area, where Falcucci et al. (2018) estimated a possible earthquake in the range M W 6.4-6.6, making this area Fig. 1 Map of the study area with the capable active faults (red lines) taken from the Fault2SHA CAD (Faure Walker et al. 2020); the historical (circles) and major instrumental earthquakes (stars) are also plotted one of the most problematic seismic gaps of Central Italy. Tondi et al. (2021) suggested that the seismogenic potential of this zone is in the order of M W = 6.0. The Campotosto area was affected by the migration of the seismicity in 2009 toward North from L'Aquila and in 2017 toward South from Amatrice; these seismic activities caused multiple moderate earthquakes (M W ≤ 5.5) nucleated in the Campotosto area. Several studies estimated the PSHA in the area of Campotosto (e.g. Meletti et al. 2006;Stucchi et al. 2011;Woessner et al. 2015;Martelli et al. 2017;Santulin et al. 2017); although the zonations used in these studies are similar, different Ground Motion Models (GMM) led to a generalized increase in the PSHA values.
In the following, after describing the PSHA conducted in this study, we estimate the seismic potential of active faults and related epistemic uncertainties of PSHA parameters, identifying the most dangerous seismic source for the dams. Such a source is used as input in PBS to model the MCE scenario by computing the ground motion for three near-source receivers located at the three dams along the edges of Campotosto Lake. We provide synthetic broadband waveforms and spectral parameters that can be used as input motion for the seismic design and retrofitting in the area.

Probabilistic seismic hazard
We evaluate the PSHA of the Campotosto region by modelling fault sources. The code used is Openquake (Pagani et al. 2014), which provides hazard curves for the desired Ground Motion Model (GMM), and uniform hazard response spectra (UHRS) for selected sites. The fault model is based on the 3-D fault system given in the Fault2SHA Central Apennines Database (CAD) (Faure Walker et al. 2020); rupture rates are parameterized by using the SHERIFS code (Chartier et al. 2019) that allows seismic ruptures involving multiple sections of a single fault or interactions between multiple faults.
The Fault2SHA code CAD (Faure Walker et al. 2020) brings together the knowledge of field geologists and seismic hazard modellers, and facilitates the definition of PSHA uncertainties.
The SHERIFS code requires the faults' geometries and the slip rates to calculate rupture rates. Given a fixed shape of the magnitude-frequency distribution (MFD) set for the entire fault system (target MFD), its iterative procedure establishes all possible earthquake rupture scenarios in the area. At each iteration, a magnitude is randomly selected based on the target MFD, and an earthquake rupture scenario that generates that magnitude is selected. For this scenario, an increase in the slip-rate budget of the affected faults is converted into an earthquake rupture rate for the magnitude under consideration. At the beginning of the iterative process, only the shape of the target MFD is set without knowing its absolute value. The iterative computation continues until the fault slip-rate budget is exhausted. The code compares the modelled earthquake rates with the local earthquake catalogue. We use for this comparison the rates obtained from a Gutenberg and Richter (GR;1944) double truncated exponential MFD model, based on the CPTI15 Parametric Catalogue of Italian Earthquakes (Rovida et al. , 2020 and the completeness periods of Meletti et al. (2019). The catalogue was declustered by applying a space and time window (Gardner and Knopoff 1974), based on information provided by the CPTI15 compilers during the project for the new Italian hazard map (Meletti et al. 2019). We assumed that the seismicity of the complete fault system also follows the form of a GR double truncated exponential MFD 1 3 model, to reproduce better the MFD shape of the CPTI15 earthquake catalogue (Rovida et al. , 2020. We started by defining the limits of the fault system in an area of about 100 km around Campotosto Lake (Fig. 1). Since the Fault2SHA CAD (Faure Walker et al. 2020) database does not cover the northern part of the study area, we integrated the system with the individual seismogenic sources of the DISS 3.2.1 catalogue (DISS Working Group 2018), collecting all the information needed by SHERIFS: the geometry (trace, dip and upper and lower seismogenic depth), the slip rate with uncertainties and the kinematics (normal, strike slip, reverse or rake) of the fault sections. We defined the length of the fault sections in order of the dimensions of the seismogenic depth, which was set at 10 km (Valoroso et al. 2013). Faults with a longer length are divided into smaller sections (however, smaller sections can break together).
We have analysed the uncertainties associated with the fault slip rate, the magnitude scale parameters, and the b-value. In SHERIFS, the uncertainty limits are defined and explored based on the number of user-defined random samples, which in our case were set to 20, each representing a model for each branch of the logic tree (see the electronic supplement Fig. ES1). For each model, a slip-rate value is chosen uniformly within the uncertainty limits associated with each fault in the database; the scale law parameters are independently selected according to a Gaussian distribution within their error limits, and a b-value is set within the range considered by the user, in our case 0.9-1.1. The uncertainties of the two magnitude scale relationships (Wells and Coppersmith 1994;Leonard 2010) in estimating the M max are accounted for by random sampling of the two-scale relations.
In summary, to account for the epistemic uncertainties, we implemented a logic tree of 120 branches, as shown in Fig. 2: the fault model is unique, as described previously; two magnitude-rupture scaling relationships (W&C94, Wells andCoppersmith 1994 andLeo10, Leonard 2010) to calculate the maximum magnitude each fault can accommodate, based on the rupture area for normal faults; twenty random samples within the lower and The logic tree adopted in our analysis consists of one seismotectonic model (Fault2SHA CAD), two scaling laws for the maximum magnitude definition (W&C94 and Leo10), twenty random samples to explore the epistemic uncertainties of the faults slip rate, the scaling law parameters, and the b-value adopted for the target MFD. In addition, the uncertainty about the GMM is explored by including three branches (Bindi11, Bindi14, Cauzzi15). The numbers indicate the branches at each node of the logic tree upper bounds of a centred triangular distribution to explore the epistemic uncertainty of (1) the b-value of the target MFD for the fault network, (2) the maximum magnitude (M max ) due to the uncertainty in the scaling law and (3) the slip rates of the faults. Also, we use the three GMM with the highest weights selected by Lanzano et al. (2020) to compute the new Italian seismic hazard map: the Cauzzi et al. (2015), the Bindi et al. (2014), and the Bindi et al. (2011) calibrated on a global, European and Italian database, respectively.
To assess the seismic hazard of the Campotosto Lake area, we use the synthetic seismicity rates resulting from the SHERIFS' computations as input for the Openquake code (Pagani et al. 2014), using the specific codes of Scotti et al. (2021). In this study, we did not analyze the influence of the different soils and considered the whole area as bedrock (Vs 30 > 800 m/s).
The most common result of a hazard analysis is the hazard curve, which represents the annual probability of exceeding certain shaking levels (Kramer 1996) in a fixed period. From the hazard curves obtained for PGA and response spectra computed at different spectral ordinates (e.g. 0.1 s, 1 s…), it is possible to obtain the uniform hazard response spectra, which have the same probability of exceeding, at all frequencies, for a fixed return period. The Campotosto hazard curves for the horizontal and vertical components are shown in Fig. 3a and b: the dispersion of each branch of hazard curves is associated with the different epistemic uncertainties analyzed. For the horizontal components, the branches of the logic tree show the largest differences when the different GMM used are highlighted, and this is more evident for long return periods (Fig. 3a).
Figures 3c and e illustrate the horizontal UHRS for the 2475 and 475 years return periods, respectively, obtained from the different branches, while Figs. 3d and f refer to the vertical components, for which only the Bindi et al. (2011) GMM is available. The obtained mean and different quantiles spectra are also represented as the weighted average of the various individual spectra of the logic tree. When analysing the logic tree branches for the horizontal components, a strong dispersion is noted in PGA, maximum spectral acceleration and its related spectral ordinate value. As can be seen both in the hazard curves and in the spectra, the strong dispersion is mainly linked to the different GMM used. The results and the shape of the mean curve are strongly influenced by the branches of the Cauzzi et al. (2015) model, especially within the PGA and 0.2 s range, while the PGA values related to Bindi et al. (2011) GMM are always the lowest for all return periods (Fig. 3a).
In a classical PSHA, all earthquakes contribute to building the seismic hazard of a site. In our analysis, we obtain the scenario that most affects the hazard of a site, expressed in terms of magnitude-distance pair, by disaggregating the probabilistic calculation. The disaggregation allows us to identify, among all the earthquakes considered to compute the seismic hazard of a region, the one that contributes most (McGuire 1995).
For a 2475 years return period (see the electronic supplement Fig. ES2a), the magnitude-distance pair that most affects the hazard of the study site is a magnitude M W 6.75 ± 0.25 within the first 10 km of distance, whereas for a 475 years return period (see the electronic supplement Fig. ES2b), it is M W 5.25 ± 0.25, within the same distance. To evaluate the contribution of each fault section to the hazard, we computed PGA curves for each rupture scenario at the site. The sum of all these hazard curves gives the total hazard at the site (Fig. 4a). The slip values of the fault sections included in the seismic hazard computation are shown in Fig. 4b.
The contribution of each section to the hazard of Campotosto is obtained by summing its hazard curves for each rupture scenario in which it is involved, normalising the sum to the total hazard curve. Figure 5a shows, for each PGA level, the sections contributing to the probability of exceedance (POE) of that PGA level. In Fig. 5b, we map the contribution of the sections to hazard at the Campotosto site (PGA 0.816 g) and surrounding area, obtained for the 2475 years return period: the hazard is mainly influenced by the nearest sections around our study site (Laga, Assergi and Barisciano sections) with similar contribution, while, for lower PGA levels, the Laga sections become prevalent (Fig. 5a), even if they are not in the highest slip rate classes (Fig. 4b). This is due to the close distance of these sections to the site. For example, comparing Sects. 2 and 3 of the Laga fault, we see (Fig. 5b) that section Laga3, which is closer to the Campotosto site, contributes the same as section Laga2, although the latter has a higher slip rate (Fig. 4b). Similarly, the   The PGA map for the 2475 years return period and the sections (thicker lines) of the faults that contribute most to the hazard of the Campotosto site (PGA 0.816 g) and surroundings for that return period compared to seismotectonic zoning approaches that smooth out the hazard in large zones, ignoring local information on active faults.

Broadband simulations using the pseudo-dynamic approach
The PBS ground motion is obtained by applying COMPEX (Moratto et al. 2015;, a hybrid approach to computing broadband seismograms combining the time series calculated in the low-frequency range by COMPSYN deterministic code (Spudich and Xu 2003) with the stochastic signals calculated in the high-frequency by EXSIM (Motazedian and Atkinson 2005;Boore 2009). The Fourier amplitude spectra of the computed velocity seismograms are merged by applying the matching filters at the intersection frequency; after that, the stochastic spectra are scaled on the deterministic spectra to reduce artificial spectral steps as described in Moratto et al. (2015). This approach has proven effective in computing realistic broadband simulations for horizontal and vertical motion components near the source (Moratto et al. 2015Parolai et al. 2020).
The wave propagation model and the finite fault parameters are the input data needed for simulating the strong ground motion. The source is modelled using the pseudo-dynamic approach (Guatteri et al. 2004;Mena et al. 2012) to emulate the temporal behavior of spontaneously propagating dynamic rupture models. Guatteri et al. (2004) demonstrated that the spatial and temporal variations in the slip, rise time and rupture propagation heavily affect the near-source ground motion. Mena et al. (2012) adopted the modified Yoffe function (Tinti et al. 2005) and upgraded the pseudo-dynamic model with a new set of relations that account for earthquake size, buried and surface ruptures and includes local rise-time variations and possible super-shear rupture speed. To increase the accuracy in shaping the asperities and the rise time of the sub-faults outside the asperities, we implemented a weighted distribution computed by the slip distribution normalized to the maximum slip of the rupture area.
To validate the correct application of the pseudo-dynamic model in our approach, we computed the source parameters and compared our values with those obtained by the waveform inversions of two earthquakes: the 2008 Iwate-Miyagi (Cultrera et al. 2013) earthquake and the 2009 L'Aquila event (Cirella et al. 2012), whose source parameter distributions were available in the SRCMOD database (Mai and Thingbaijam 2014). For the 2008 Iwate-Miyagi event, our computed rise time has a median difference of -0.05 from the inversion one; the differences of the peak of the slip-velocity in the source time function (V peaks ) have a median of -0.20. In the case of the L'Aquila earthquake, we obtain more significant differences (median of 0.14) with respect to the inversion results of Cirella et al. (2012), probably because of the heterogeneity in the velocity model; the differences in V peaks values also have a median of 0.14. In both cases, the rupture time distribution reproduces satisfactorily the results retrieved from the rupture inversions.

The Campotosto scenario input parameters
The selection of the finite fault input parameters for the PBS at Campotosto is driven by considering the MCE scenario case that, in our case, is provided by the probabilistic disaggregation analysis and geological evidence. The DISS database (DISS Working Group 2018) reports the Campotosto fault as the composite normal seismogenic source of Colfiorito-Campotosto (M max = 6.5); further, DISS (2018) and the ITHACA (2019) database report a debated seismogenic source of Mt. Laga based on geological evidence (Blumetti el al. 1993) and associated with an M max = 6.6, a striking angle of 150° and a steep dip. Boncio et al. (2004) linked the Mt. Laga system to the more detailed Mt. Gorzano fault, considered an active WSW dipping normal fault; in particular, the field observations map the Mt. Gorzano structure as a continuous fault surface 28 km long. In some stretches, the fault scarp is visible and well exposed in the central part of the structure; it strikes N140-150° and dips 60-70° to the SW. Boncio et al. (2004) estimated an M max = 6.7 for this fault. The same magnitude was assumed by Pace et al. (2006) in designing the seismogenic boxes to compute seismic hazard. Therefore, we model a 28 km long and 13 km wide rupture area capable of generating an M W 6.7 (M 0 = 1.26×10 19 N*m) earthquake (Fig. 6). The fault strikes N150° and dips 65° to the SW (Fig. 6) in agreement with ITHACA (2019) and Boncio et al. (2004); the top of the rupture area is fixed at a depth of 1 km, while the simulated earthquakes nucleate at 11-12 km and propagate toward the surface. The slip distribution is a random field (Mai and Beroza 2002) stochastically consistent with slip distributions of past earthquakes. We generate five different slip distributions with seven different nucleations placed at the bottom of the fault from NW to SE edge (Fig. 6) to simulate unilateral and bilateral rupture, resulting in 35 input source models (see the electronic supplement Fig. ES3). Our source models encompass the variability observed on the maximum slip (280-450 cm) within the predicted range (Thingbaijam and Mai 2016). Such variance in the maximum slip leads to the generation of different asperities and rupture areas, producing variability in the effective stress drop. In contrast, the static stress drop remains constant, being fixed by M W and the entire fault area.
The crustal wave velocity model employed for the deterministic computation is taken from Poiata et al. (2012). The high-frequency stochastic simulations require tuning the wave propagation parameters; the high-frequency propagation in Central Italy has been investigated in various recent studies (e.g. Pacor et al. 2016;Bindi and Kotha 2020). We use the parameters proposed by Malagnini et al. (2011), which have been validated on the moderate earthquakes of the Campotosto 2017 sequence (see next section). We do not include site effects in the time series computation.

Validation of the input parameters
To validate the source parametrization and the propagation parameters to be used for the Campotosto scenario, we compute the near-source synthetic seismograms for three moderate earthquakes that occurred in the study area on January 18, 2017, and we compare them with the corresponding records downloaded from the ESM-Engineering Strong-Motion database . In particular, the events were recorded in the nearsource by the Italian Accelerometric Network-RAN (Zambonelli et al. 2012) stations Poggio Cancelli (PCB) and Sella Pedicate (SPD). The finite fault parameters used in the simulations are taken from Vicic et al. (2020) and listed in Table 1.
We show the computed accelerogram and compare the recorded signals in Fig. 7 and the electronic supplement ES4, and ES5. We reproduce the recorded low-frequency displacement satisfactorily at PCB and SPD in amplitude and phase for the three events (Fig. 7a, and in the electronic supplement ES4a, ES5a). Some differences observed at the high frequency of the horizontal component amplitude spectra of PCB (see the electronic Fig. 6 The rupture area of the Campotosto fault (green rectangle), the seven nucleations of the rupture (green stars) and the three receivers (yellow triangles) considered in the scenario modelling. The associated focal mechanism (Boncio et al. 2004) is plotted in the top right inset. We also show the capable active faults (red lines)  Fig. ES5b), placed in the centre of the filling valley (Antonielli et al. 2021), can be due to site effects, as shown by the horizontal to vertical spectral ratio available from the ESM database . Other minor discrepancies can be related to unmodelled local propagation effects (e.g. topographic and valley edges). Finally, the accelerograms obtained from the combination of deterministic and stochastic simulations reproduce all the components of the recorded signals (Fig. 7c, and in the electronic supplement ES4c, ES5c).

Campotosto scenario
We compute the scenario for an M W = 6.7 earthquake in the Campotosto area at PCB and SPD located above the rupture area and to Rio Fucino (RFC), another station belonging to Fig. 7 Comparison between recorded (black lines) and simulated waveforms (red lines) for the earthquake (M W = 5.2) that occurred on January 18, 2017 at 09:25. a Recorded vs deterministic displacement signals filtered between 0.08 and 0.5 Hz, b recorded vs stochastic accelerometric spectra, c recorded vs broadband simulated waveforms the RAN network, which is very close to the rupture area (Fig. 6); these three sites represent three dams, which confine the Campotosto artificial water reservoir. The peak ground motions (PGMs) mean values obtained using the 35 source models are reported in Table 2. The results are given for each motion component of the receivers (see the electronic supplement Fig. ES6). We also provide the geometrical mean of the horizontal components.
The computed PGA, within the related standard deviation, are compatible with Akkar et al. (2014a, b); the PGV (Peak Ground Velocity) values are more significant than GMM by Akkar et al. (2014a), also considering one standard deviation and the vertical PGV value at SPD and PCB are larger than horizontal ones.
The horizontal PGA and PGV values computed for the three receivers (Table 2) are comparable with the strong motion recorded during the 2016 Central Italy sequence when PGA reached 0.87 g while PGV exceeded 80 cm/s at some near-source stations (Suzuki and Iervolino 2017). Pitarka et al. (2021) simulated ground motion for the 2016 Norcia earthquake, with PGV reaching values of up to 200 cm/s close to a large slip patch. On average, the most considerable shaking is estimated at PCB, located above the rupture area and strongly influenced by near-source effects (Fig. 6). The simulated vertical ground motion is comparable to the shaking recorded at the near-source station located in Castelluccio di Norcia (CLO) during the M W = 6.6 2016 Norcia earthquake (D'Amico et al. 2019); similarly to PCB, CLO was situated above the rupture area of a strong distensive earthquake.
The horizontal Peak Ground Displacement (PGD) values are consistent with the GMM prediction by Goldberg et al. (2021); the vertical PGD at PCB and SPD is higher than the horizontal values.
The Static Offset (SO) is negligible at RFC, located at the border of the rupture area, whilst it is relevant in the vertical components of PCB and SPD. However, the values are comparable with those obtained by the empirical relation of Kamai et al. (2014). The SO and the vertical response spectra simulated at PCB are also similar to those computed at CLO for the 2016 Norcia earthquake (D'Amico et al. 2019).
The horizontal SA are generally comparable with Akkar et al. (2014a) within the standard deviation (Fig. 8a, c, e). In comparison, the values of the vertical component (Fig. 8b, d, f) are higher than Akkar et al. (2014b) for all the receivers. To compare the simulated spectra with recorded ones, we retrieved the spectra of the 2016 M W = 6.6 Norcia earthquake (stations CLO and T1214) and the 2016 M W = 7.0 Kumamoto earthquake (stations KMM10 and KMM16) from the Ness2.0 database (Sgobba et al. 2021) at stations above the rupture area. In addition to a similar magnitude, these earthquakes were selected for their focal mechanism and vertical SO < -40 cm. Considering the standard deviations, our results are larger than recorded vertical spectra at periods of 0.8 s at PCB and SPD, stations located above the rupture area and strongly influenced by single asperities in the slip distributions and heterogeneities in rise-time and rupture velocity distributions within the source models. The simulated horizontal displacement spectra (SD) are compared ( Fig. 9) with the values estimated by GMM (Cauzzi et al. 2015) adjusted for fling-step factor at longer periods (Schiappapietra et al. 2022). Our mean displacement spectra are comparable with adjusted GMM within the standard deviation for RFC and SPD, while SD simulated at PCB overestimates the empirical values at 2.0 < T < 6.0 s; this difference is due to the position of the source that is shallow and close to the receivers. Also, the mean simulated SD are comparable with the recorded spectra in the whole period range for horizontal components; the mean SD simulated on the vertical component is comparable with the recorded spectra at SPD, while it is moderately larger at PCB (located above the central part of the rupture area) and lower for T > 4.0 s at RFC, which is located outside from the rupture area.

Discussion
The PSHA ground motion estimated for the study area can be extracted from the hazard curves for a few return periods. The Italian seismic hazard map (Meletti et al. 2006;Stucchi et al. 2011) defines a PGA of 0.260 and 0.448 g expected at return periods of 475 and 2475 years, respectively; more recent studies estimated higher values of PGA due to the different updated GMM used. Our PGA values (0.414 and 0.816 g, respectively) are comparable with these recent studies, although the complete logic tree is very different. The similarity of the results is probably linked to the choice of the same GMM, or at least of the same generation, despite a completely different approach to seismological sources. Although comparable acceleration values are obtained, the inclusion of seismogenic sources results in a more defined areal distribution of accelerations than when areal sources are used. This is a significant result in the near-source, obtained mapping hazard based on faults.
In general, the MCE-case PBS scenario ground shaking is comparable with the empirical GMM estimations within one standard deviation. The variability of input source parameters is limited, being M W , rupture area and the focal mechanism fixed on the basis of known tectonic settings and quite good knowledge of various seismogenic sources; these fixed parameters are satisfactorily constrained because Central Italy is a well monitored and studied area (see Sects. 1 and 3.2 for references). However, these limitations are reliable constraints for our results. The simulated shaking is very sensitive to source parameters such as spatial heterogeneity of slip, rise time and rupture velocity on the fault volume; the seismic moment distributions and asperities patches can generate very high strong motion in near-source, particularly on the longer period parameters and with significantly associated uncertainties (Moratto et al. 2017;Causse et al. 2021;Pischiutta et al. 2021;Pitarka et al. 2021). For these reasons, the scenario generated utilising 35 rupture models evidences a substantial variability; the standard deviations associated with the mean values are larger than those related to GMM, and they can exceed 50% of the mean values. We evidence that our simulated PGA values (Table 2) cover a range comparable with the PGA variability observed in near-source during the 2016 Norcia earthquake . The highest values are generated by the bilateral ruptures associated with the seismic moment distributions and the largest slip in our scenario. The shaking computed at various receivers is also influenced by the source approximation, the distance and the directivity effects of the rupture propagation; the updip propagation generates strong directivity effects around the surface exposure or its up-dip projection (Somerville 2003). In our scenario, the rupture area is modelled as a planar fault, but some authors (Valoroso et al. 2013;Vicic et al. 2020) suggest a listric fault; the simplification of the fault geometry could increase PGV values on the hanging wall and a decrement on the footwall (Passone and Mai 2017). The amplitude ratio between the three motion components is also linked to different radiation patterns generated by various source models. Therefore a translation of the fault plane with respect to the original model could generate results different from those computed in this study.
The increase of the recording stations and the available seismograms evidences the relevance of the vertical motion for the shaking observed in the near-source (e.g., Di Sarno et al. 2010;Suzuki and Iervolino 2017;D'Amico et al. 2019); at the same time, the development of the computation models allow to simulate accurately these effects giving a good agreement between the theoretical models and the observations. In our scenario, the most significant shaking values in the vertical component are observed at longer periods for receivers located above the fault (PCB and SPD). In this case, the normal-type source generates a strong subsidence effect. The ratio between the vertical and horizontal ground acceleration peak exceeds the unity in many other cases (e.g., Elgamal and He 2004;Shrestha 2009;Kim et al. 2011).
The validation process is fundamental to obtaining a reliable estimation of the seismic hazard and to assessing the seismic safety of large structures (e.g., dams), which may be strongly influenced by spatial variability in one or both the horizontal directions of the input motion (Dello Russo et al. 2017). Grimaz and Malisan (2014) evidenced a deficiency of seismic safety in the near-source if the design does not consider the seismic effects.

Conclusions
Reliable estimates of the ground motion in the near-source are crucial to evaluating the seismic hazard, mainly if some critical facilities operate there. Ground motion models are weakly constrained in the near-source because of the lack of recorded waveforms. Numerical simulations can fill the gap and provide the input design for specific sites.
We apply the PSHA and, later, we compute the PBS scenario to estimate the nearsource seismic hazard for three dams closing the water impoundment in the Campotosto area (Central Italy). The PSHA is performed following the approach of seismotectonic probabilism by modelling fault sources. The disaggregation of the probabilistic seismic hazard, calculated for the return period of 2475 years, provides the magnitude-distance pair of M W = 6.75 ± 0.25 within the first 10 km of distance. This magnitude-distance pair 1 3 represents the most influencing earthquake for the hazard estimate in the study area. The PBS computation assumes as the MCE case an M W = 6.7 earthquake, consistently with field geology studies (e.g., Boncio et al. 2004) and provides the ground motion for three receivers at the three dams.
We reproduce the variability due to the source rupture parameters that affect the simulated ground motion and that are not reproduced by the GMM in the near-source. This variability is related to the rupture propagation model heterogeneities and reliable finite fault parameters input in the modelling. The most significant ground motion values are obtained above the rupture area, at the receiver PCB, particularly at long periods. At the same time, high-frequency shaking characterizes the receiver RFC, placed at the edge of the main fault. On average, the associated variability can reproduce the strong motion characteristics recorded in near-source for similar earthquakes.
The Campotosto scenario shows that the vertical component might experience shaking values comparable with the horizontal, and, therefore, it cannot be ignored in near-source studies. Our vertical peak parameters and spectra are consistent with the values recorded for the Norcia 2016 mainshock, and the observed low frequency contribution (PGD and SO). GMM probably underestimates the vertical component's ground motion, especially at longer periods and above the rupture area. Numerical simulations can exceed the limits of the GMM in the near-source to fit observations and define reliable seismic input. However, ground motion simulations require a constrained definition for the seismic source, the propagation and local site effects; adequate knowledge of these parameters' bounds is a primary key to defining reliable, strong motion and uncertainties of the results. Future developments will be focused on the site-specific hazard analysis, in particular for the embankment dam at PCB.