1 Introduction

Seismic hazard assessment is a fundamental task to ensure public safety because it provides the seismic parameters needed by designers to mitigate earthquake damage to structures and the environment. Strategic facilities are usually the most critical infrastructures, and detailed hazard estimations are routinely performed for nuclear power plants, dams, gas transit lines, and pipelines (e.g., Toprak and Taskin 2007; Schmid and Slejko 2009; Slejko et al. 2011; Renault 2014; Grimaz and Slejko 2014a, b; Dello Russo et al. 2017; Tamaro et al. 2018).

Assessing the hazards posed by faults that cross pipelines is crucial both during the design and life cycle of a pipeline project. The consequences of potential pipeline failures are usually severe. Earthquake damage to buried pipelines can be caused by transient ground-shaking effects and permanent ground deformations, or a combination of both. The related permanent ground deformations occur due to surface faulting, liquefaction, landslide activity, and differential settlement from the consolidation of cohesionless soil (Toprak et al. 2015). They create large stress concentrations at structural joints (Wham and O'Rourke 2015). The methodologies for estimating potential pipeline damages are mainly based on empirical relationships obtained by past earthquakes that relate pipeline damage to ground motion parameters (e.g., Pineda Porras et al. 2010; O'Rourke et al. 2014). However, to account for different source-types and wave propagation effects and develop robust regressions for future fragility analyses, it is necessary to add new data as they become available after earthquakes (Toprak et al. 2015).

To effectively include earthquake effects in the design phase, a specific seismic hazard assessment of the area around a pipeline is needed by combining both probabilistic and deterministic approaches (Slejko et al. 2021).

Probabilistic Seismic Hazard Analysis (PSHA) considers all the seismological information for a given area (i.e. earthquakes and faults). It estimates the probability that various levels of earthquake-caused ground motion will be exceeded at a given location at a given future time. A detailed description of the method is given in Slejko et al. (2021). However, PSHA is rarely constrained in the near-source areas, where the observed variability in ground-motion intensity is larger than that at greater distances (e.g., Shakal et al. 2006; Mai 2009). The ground motion empirical relationships used by PSHA suffer from a lack of near-field recordings and often do not adequately represent the observed source effects (Dalguer and Mai 2011). Therefore, numerical simulations are a viable solution that account for the physical characteristics of sources, propagation paths, and site effects.

Deterministic seismic hazard analysis (DSHA) evaluates the strong ground motions at selected sites for a specific earthquake scenario (or a combination of earthquakes). In each scenario, an earthquake is generated by a fault that is modelled using its source parameters. In the latter approach, it is crucial to define a worst-case scenario (Kramer 1996) that is obtained by computing the ground motion peaks and integral parameters at the shortest source-to-site distance; the goodness of the worst-case scenario tied to the geological and seismotectonic knowledge of an area and in some case may not be at all possible. The input data needed for calculating strong ground motions are usually a structural model, soil (local geology) site responses, and finite fault parameters. As expected, accurate knowledge of the input parameters results in more reliable and robust estimates of the considered DSHA scenarios. DSHA can be helpful for designing critical structures, while PSHA can be applied for preliminary evaluations (Krinitzsky 2003). Overall, McGuire (2001) asserted that determinism versus probabilism is not a bivariate choice but rather a continuum in which both analyses are conducted (Bommer 2002).

In this paper, we present the DSHA that we performed for a compressor station along the Trans Adriatic Pipeline (TAP), which is an infrastructure that transport natural gas from Greece through Albania and the Adriatic Sea to southern Italy and further to western Europe; this is a complementary work of the PSHA study performed by Slejko et al. (2021). By employing DSHA analysis, we modelled the worst-case scenario, calculating the strong ground motion generated by two different earthquakes in the area of a compressor station located close to the town of Fier along the Albanian coastal region (Fig. 1). The PSHA disaggregation analysis (Slejko et al. 2021) identified the most critical (highest hazard) faults for the compressor station and provided input to formulate the corresponding physics-based scenarios for the DSHA analysis.

Fig. 1
figure 1

Map of the FCS area depicting the projection of the employed rupture areas onto the ground surface (rectangles), their expected surface traces (coloured dotted lines), considered nucleation points (stars) and receivers employed for the simulations (triangles); ALCS220 and ALCS230 are shown in red and green, respectively. The black dotted lines show the position of the upper edge of other faults referred to in the text (KSD = Karaburuni-Sazani detachment; NKT = North Kurveleshi thrust; AR = Ardenica backthrust; and LTFZ = Lushnje transfer fault zone)

From the broadband synthetic seismograms computed at the compressor station for broad frequency content, we retrieve the ground motion peak measures (PGM), i.e., Peak Ground Acceleration (PGA), Peak Ground Velocity (PGV), Peak Ground Displacement (PGD) and the spectral acceleration (SA), which are usually used for relating damage rates to particular earthquakes (Pineda Porras et al. 2010). We also compute the permanent displacement along the pipeline (static offset (SO)), which is usually considered the potentially most severe cause of damage to pipelines (Falcao Silva and Bento 2004). All the waveforms are first computed for stations on bedrock, and successively, site amplification effects are included in the calculations by applying specific transfer functions determined from the local geology.

2 Seismotectonic view of the Fier area

The TAP route crosses an area characterized by a complex tectonic setting that includes, from east to west, the North Anatolian strike-slip fault zone, the northern Aegean extensional realm, the Hellenides and Albanides fold-and-thrust belts, and the less deformed carbonate Apulia Foreland at the future pipeline's end. All of these tectonic structures are active and capable of generating strong earthquakes.

The Fier compressor station (FCS hereafter) is located on the Adriatic coast, 70 km SW of Tirana in the Albanian territory. The area is characterized by a convergent tectonic regime that generated an accretionary wedge in continuity with the southernmost Hellenides and northernmost Dinarides (e.g., Robertson and Shallo 2000; Nieuwland et al. 2001; Jouanne et al. 2012). Major east-dipping thrusts, minor back thrusts and some transfer zones (e.g., Roure et al. 2004; Lacombe et al. 2009; Jouanne et al. 2012; Caputo and Pavlides 2013) contribute to segment, and hence, reduce the overall size of the thrust system’s decreasing maximum magnitudes.

For DSHA, the Vlore backthrust (ALCS220 from now on) and the Fier thrust (ALCS230 from now on) are considered the most relevant seismic sources due to their proximity to the FCS, the associated maximum magnitude associated with these faults, and the recurrence time interval (Table 1), as also identified by the disaggregation analysis of the PSHA results (Slejko et al. 2021).

Table 1 Table 1

The ALCS220 fault is located on the eastern side of a topographic high and separates the coastal lowland from the Adriatic Sea (Fig. 1). This backthrust is highlighted in several interpreted seismic profiles (Nieuwland et al. 2001; Argnani 2013) and is consistent with the compressional focal mechanism of moderate events occurring in the crustal volume, including the Vlore backthrust. The width of this west-dipping fault is constrained downdip by the Karaburuni-Sazani detachment, which represents the lowermost boundary of the orogenic wedge. The fault width is possibly further limited in the depth direction by secondary fault branches splaying up from, and synthetic with, the primary detachment. Concerning the fault length, the lateral tips are limited by the Lushnje Transfer fault zone to the north and, due to its opposing strike, the Karaburuni-Sazani detachment to the south. As a consequence, the associated maximum magnitude is most likely constrained to 6.7.

The ALCS230 fault belongs to the compressional belt developed as a consequence of the Adria-Albanides collision, the activity of which has been well documented by the seismicity characterizing the area. This seismogenic source is mainly detectable in seismic profiles (Velaj 2011; Skrami 2001; Aliaj 2006), showing that the Ardenica backthrust cuts the Fieri thrust at a depth of ca. 5 km. The Ardenica backthrust is a shallow secondary, west-dipping tectonic structure. The fault should join the basal detachment at a depth of ca. 12–13 km (i.e., maximum depth), while its length is constrained by the Lushnje transfer fault zone to the north and is considered to terminate in correspondence with a large overstep zone to the south, possibly connecting this structure to the North Kurveleshi thrust (Fig. 1). The corresponding maximum magnitude is likely constrained to 6.6.

3 Method

To compute the broadband synthetic seismograms used in this study, we adopted the hybrid approach described in Moratto et al. (2015) by merging low-frequency waveforms with high-frequency stochastic signals. In particular, the low-frequency signals were computed using the discrete wavenumber technique COMPSYN (Spudich and Xu 2003) up to 2 Hz, with the high-frequency stochastic synthetics calculated by EXSIM (e.g., Motazedian and Atkinson 2005; Boore 2009) at a lower cut-off frequency of 0.1 Hz. We merged these seismograms in the Fourier amplitude domain utilizing the matching filters at the target intersection frequency, which is usually set between 0.5 and 1.5 Hz to optimize the procedure; further details are described in Moratto et al. (2015). The extended fault was modelled following a pseudodynamic approach, while the slip distribution is a random field (Mai and Beroza 2002) that is consistent with the slip distribution of past earthquakes. Given the slip on the rupture area, the pseudodynamic approach estimates the peak slip velocities, rupture velocities, and rise time distributions needed to compute the ground motion through the empirical relationship proposed by Guattieri et al. (2004). The rupture velocity is typically ranged between 0.7 and 0.8 of the S-wave velocities. The corresponding stress-drop and fracture energy distributions are derived from the previously generated slip; hereafter, the distributions of the rupture velocity, rupture time, peak slip velocity, rise time, and pulse width are computed.

We modelled different rupture nucleation points and seismic moment distributions since these parameters affect the rise time and rupture velocity in the pseudodynamic model (Guattieri et al. 2004) and thus the simulated ground motion (e.g., Moratto et al. 2017; Cultrera et al. 2010). The PGM mean values (with their related standard deviation) are then extracted, and the related uncertainties are estimated. We also compute the SO along the pipeline using a wavenumber integration method (Spudich and Xu 2003).

Seismic waveforms are computed for a grid of receivers to reproduce and oversample the ground motion spatial variability in the near field (Lin et al. 2005; Douglas 2007; Moratto et al. 2017). Multiple receivers allow us to account for uncertainties due to source mislocation that affects the source-receiver distance, evidencing possible azimuthal variations due to the source radiation pattern, directivity, and other parameters that could influence ground motion.

Since no experimental measurements of site response or resonant frequency are available for the FCS area, and because there is no knowledge of a possible sharp velocity contrast with depth, we estimate the site amplification using an equivalent linear analysis code (equivalent linear earthquake response analysis (EERA) of layered soil deposits, e.g., Bardet et al. 2000). The synthetic seismograms related to the bedrock scenarios are finally convolved using the average response obtained for the different S-wave profiles. To evaluate the vertical amplification, we use a procedure suggested by Beresnev et al. (2002). For practical purposes, the design motions are synthesized as nearly vertical (SV) waves propagating across the specified soil column. This procedure can be implemented through one-dimensional horizontal-component modelling using EERA and empirical depth-dependent V/H correction factors to account for inclined propagation paths.

4 Input parameters

We set up two different scenarios considering the ALCS220 and ALCS230 faults as the most critical fault for our case study (Table 1). All of the input parameters used to simulate the strong motion have been defined based on the field surveys, literature data, and geophysical investigations carried out in the framework of the TAP project, as described in detail by Slejko et al. (2021).

The strong ground motion is computed for a grid of receivers (Fig. 1) with an interstation distance of 3 km, and the PGMs are extracted for each receiver. The receiver grid, placed at a maximum Joyner–Boore distance of 30 km from the fault, follows the TAP route, and it is 27 km long in the W–E direction and 15 km wide in the N–S direction, covering an area of 30 × 18 km2. The grid orientation and position depend both on the fault strike and pipeline route direction. The FCS receiver configuration is fixed in the two scenarios. The anelastic attenuation at low frequencies is considered to be negligible for the involved source-receiver distances.

4.1 The seismic sources

Geological observations, integrated with geodetic, seismic, and seismological data, contributed to defining the main parameters used to model the ALCS220 and ALCS230 faults (Table 1). The maximum possible magnitude (Mmax) is obtained from the empirical relationship of the rupture area versus magnitude (Wells and Coppersmith 1994), while assuming the worst-case scenario in which the entire fault surface breaks and slips. A comprehensive description of the fault models for the entire TAP route is presented in the companion paper by Slejko et al. (2021). All source parameters of Table 1 are retrieved from Appendixes A1 and A2 in Slejko et al. (2021), except the maximum displacement (Dmax), which was calculated in agreement with the Thingbaijam and Mai (2016) relation. The average stress drop is assumed to be 2.5 MPa, according with values estimated worldwide (Allmann and Shearer 2009) and for the Mediterranean region (Kostantinou 2014). The minimum distances of the FCS from the ALCS220 and ALCS230 faults are approximately 7.5 km and 11 km (Fig. 1), respectively.

For each fault, five slip distributions are generated (Fig. 2) by adopting the hypothesis that Dmax varies within the range predicted by Thingbaijam and Mai (2016). For each slip model, we consider three different rupture nucleation points in the deeper part of the fault (Fig. 1). In this way, we have generated fifteen source models for each fault, which are used for both the low- and high-frequency stochastic simulations. The slip distributions with the largest Dmax value are related to well-delimited asperities, whereas those with the smallest Dmax value are related to the smoothed distribution of the slip. The peak slip velocity, rupture velocity, and rise time are computed starting from the slip distribution (Guatteri et al. 2004).

Fig. 2
figure 2

Slip distributions utilized in the simulations of the ALCS220 (on the left) and ALCS230 (on the right) faults; magenta stairs indicate the positions of nucleation points on the rupture area

4.2 The velocity model and propagation parameters

We used the average 1D velocity model retrieved from the 3D tomography of Papazachos et al. (2002), which is currently the best available source of information for the velocity structure of the broader Albania area (Fig. 3). The S-wave velocity at 0 km depth (2.5 km/s) is considered the bedrock velocity, which partly reflects the characteristics of deep consolidated sediments in the coastal areas of western Albania. The density model has been derived from the Nafe–Drake relationship (Ludwig et al. 1970).

Fig. 3
figure 3

1D average models of P- and S-wave velocities for the Fier area

To compute the high-frequency content of the stochastic waveforms, we utilized the parameters (Q value, geometrical spreading, duration function) proposed for the Hellenic area (Margaris and Boore 1998; Margaris and Hatzidimitriou 2002) since no high-frequency propagation models were available for the Albania region:

  • Q:Q = 275(f/0.1)−2.0 for f ≤ 0.2 Hz, Q = 88(f/1.0)0.9 for f ≥ 0.6 Hz, Q for 0.2 < f < 0.6 Hz determined from power-law fit to values of Q at f = 0.2 and f = 0.6 Hz.

  • Geometric spreading: 1/R.

  • Duration: (1/f0) + 0.05R.

4.3 The site effects

To estimate the FCS amplification factor, we generated different S-wave velocity models to be used with EERA. The linear equivalent response is preferred because no detailed geotechnical parameters or strong ground motion records are available. Duni and Kuka (2011) estimated an S-wave velocity value for the shallow 30 m depth interval (VS30) of approximately 170 m/s, thereby classifying the FCS site as a D ground type according to Eurocode 8 (CEN 2004) and an E soil class according to the International Building Code (IBC 2003). The VS30 profile is defined down to 50 m, and no sharp S-wave velocity contrast is found within the shallower sediments. Geotechnical investigations in the Semani area near Fier (see Allkja and Dhimitri 2012) show that medium dense sand with silty clay, silty sand, and occasional peat inclusions characterize the upper 50 m.

Site response estimation, due to the very low S-wave velocities at the surface, requires exploring all possible velocity models at various depths. As no velocity models were available for depths below 50 m, to obtain a fair estimate of the local site response, we created a set of models characterized by different velocity contrasts at variable depths and gradients while keeping the surface velocity structure fixed for the upper 50 m. For this purpose, we included various possible S-wave velocity contrasts at different depth interfaces (50, 100, 200, 400, and 600 m) and considered the possibility of a smooth increase in the S-wave velocity up to a depth of 1 km (Fig. 4a). The EERA considers the material shear modulus and damping constant within each layer and for the entire duration of the earthquake. In this study, the shear modulus and damping ratio were retrieved from Idriss and Sun (1992) and Vucetic and Dobry (1991) for sands and silty clays, respectively. The shallow stratigraphy at the FCS site was retrieved from Allkja and Dhimitri (2012).

Fig. 4
figure 4

a Amplification factors versus frequency from EERA. The blue and red curves represent the average response of the horizontal and vertical components, respectively. Dashed lines represent ± 1σ. The inset plot in the right part of the figure presents the velocity models considered. b Comparisons between the simulations computed for the bedrock and with the incorporation of site effects that consider the linear response and the linear equivalent models

To evaluate the vertical amplification, we used the empirical depth-dependent V/H correction factors developed by Darragh et al. (1999), specifically for Eurocode 8 D sites. The strain-dependent soil properties, serving as input to EERA, are generally assumed for S-waves. In Fig. 4a, the mean horizontal and vertical amplifications estimated by the suite of velocity models are shown. Such amplification curves enhance the ground motion of the horizontal and vertical components in the range from 0.5 to 3.5 Hz at the FCS site. Higher frequencies are strongly damped for all S-wave velocity models due to the D-type soil and significant thickness of the soft soil deposits.

For softer soils (D and E soil types), the equivalent linear solutions can have up to 30–50% larger amplitude peaks below 2 Hz than nonlinear formulations because they include no time dependence in the solution (Hartzell et al. 2004). Thus, in the case of the nonlinearity of soil under strong ground motions, 1D equivalent linear modelling cannot correctly reproduce hysteric soil behaviours. The possible underdamped seismic response below 2 Hz is handled by averaging the amplification between the set of considered S-wave velocity models. Figure 4b shows the comparison between the waveforms computed under bedrock and soil conditions by applying the amplification function shown in Fig. 4a.

5 Results

The mean values of the two scenarios for PGA, PGV, PGD, and SO computed at the FCS for bedrock and soil conditions are shown in Table 2.

Table 2 DSHA results for the FCS (peak mean values and standard deviation obtained by the 15 utilized

In general, we observe that the north–south component values are larger than the east–west component values. In contrast, vertical shaking reaches peak ground motions similar to horizontal shaking, as observed in past earthquakes (e.g., for central Italy 2016, see Luzi et al. 2017; D'Amico et al. 2019). The vertical component in the near field shows the largest values as the result of the combined effect of the radiation and directivity of the thrust fault type. The ground motion shows a strong dependence on fault location, and the PGM decreases as the distance from the source increases. As a consequence, the uncertainties in the fault position could play a crucial role in the variability of the results. This is shown in Figs. 5 and 6, where we consider the two faults (ALCS220 and ALCS230), respectively, and the maximum PGA decreases from 0.6 g close to the sources to 0.1 g at the more distant receivers (Fig. 5a), while the vertical component can reach 1 g above the rupture area (Fig. 6b). The PGV decreases from 60 cm/s to 20 cm/s leading away from the source (Fig. 5c), while the vertical component can reach 100 cm/s above the rupture area (Fig. 6d). PGD varies between 10 and 30 cm (Fig. 5e), with a vertical component of 60 cm in the near-source region for ALCS220 (Fig. 5f). For the case of ALCS230, the horizontal PGD is lower than 30 cm (Fig. 6e), but it can reach 90 cm for the vertical component (Fig. 6f). The SO is estimated to be approximately 5 cm at the FCS (Figs. 5g and 6g) and increases to 15 cm above the fault plane for ACLS220 (Fig. 5h) and to 50 cm for ACLS230 (Fig. 6h).

Fig. 5
figure 5

PGA, PGV, PGD, and SO computed for the rupture scenario of fault ACLS220; the largest of the horizontal components is plotted in the left column, while the vertical component is shown in the right column. The blue square depicts the FCS site; the black triangles are the virtual receivers used in the simulations, while the blue open rectangle shows the projection of the rupture area; the grey line displays the pipeline route

Fig. 6
figure 6

Same as Fig. 5 for the ACLS230 fault rupture scenario

A comparison of the PGA and PGV computed using the Ground Motion Model (GMM) is shown in Fig. 7 for the ACLS220 scenario and Fig. 8 for the ACLS230 scenario. We considered the GMM proposed by Campbell and Bozorgnia (2014; CB14 hereafter) for the horizontal component and that of Bozorgnia and Campbell (2016; BC16 hereafter) for the vertical component. For ACLS220, both the horizontal PGA and PGV mean values are comparable with those proposed by CB14 within the related standard deviations (Fig. 7a, b), whereas the vertical component of PGA, particularly the PGV, exceeds the associated + 1σ of BC16 (Fig. 7c, d), mainly in the near field. The simulated PGV values attenuate with the source distance, computed as the distance from the rupture area (RRUP > 20 km), more slowly than the empirical relationship. Considering the ACLS230 scenario, the horizontal PGA and PGV values (Fig. 8a, b) underestimate the same GMM (within -1σ), and at a fixed distance, their scattering can exceed one order of magnitude. The vertical PGA values fit those of BC16 (Fig. 8c), while the PGVs overestimate the GMM (Fig. 8d). The discrepancies observed in the comparison between the simulated peaks and GMMs of the two scenarios are mainly related to the diverse position of the receiver grid relative to the two faults considered (Figs. 5 and 6), as well as the different depths of the two ruptured areas (Table 1) and their directivity effects. Therefore, the different source models lead to a large scattering of the simulated peak parameters, and the standard deviation values in some simulations are similar to the mean values, while the standard GMMs have a standard deviation of approximately 30%.

Fig. 7
figure 7

Comparison of the ACLS220 fault rupture scenario between the GMM (blue lines) and the simulations: a horizontal PGA vs. CB14; b the horizontal PGV vs. CB14; c the vertical PGA vs. BC16; and d the vertical PGV vs. BC16. The geometrical means of the horizontal PGA and PGV components are orientation-independent computed (i.e., RotD50, Boore et al. 2006). The rupture distance (RRUP) is plotted along the x-axis. Black circles are the values at the receiver computed for each simulation, while red circles are the mean values for each receiver

Fig. 8
figure 8

Same as Fig. 7 for the ALCS230 fault rupture scenario

The comparison of the horizontal response spectra (SA) computed at the FCS is in good agreement with the CB14 prediction (Fig. 9a). ALCS230 has a deeper fault top edge than ALCS220, and therefore, the corresponding mean spectra are quite lower than that of CB14 (Fig. 9b). The mean SA computed for the vertical component is constantly larger than that of BC16 for T > 0.3 s, while the GMM experiences a peak at 0.1 s (Fig. 9c, d); the amplification of the simulated spectra in the range of 0.5 s < T < 1.5 s can also be related to the site effects (Fig. 4a) that affect the vertical ground motion (particularly PGV) in this period range. At periods T > 2.0 s, the ALCS220 simulations overestimate that of BC16 (Fig. 9c). The single simulated spectra (black lines in Fig. 9) show a large ground motion variability associated with the different source models utilized in the DSHA computation.

Fig. 9
figure 9

Comparison between GMMs (blue lines) and DSHA mean spectra (red lines) at the FCS for ACLS220 (a, c) and ACLS230 (b, d); a and b represent the horizontal component compared with CB14, while c and d represent the vertical component compared with BC16. The horizontal geometric mean is computed as RotD50 (Boore et al. 2006). The black lines show the SA extracted from each DSHA simulation

Our static offset values are compared with the empirical relationships proposed by Kamai et al. (2014). In the comparison, we consider only the receivers placed above the ALCS230 rupture area (RJB = 0; Fig. 6) where the highest vertical static offset value is observed. In the ALCS220 scenario, the static offset values are generally lower, and only three receivers are placed close to the projection of the rupture area on the surface. The empirical relationships by Kamai et al. (2014) predicted a vertical static offset of 25.9 cm (14.2–47.2 cm, if we consider one standard deviation) for the receiver placed above the rupture area; our simulations give comparable values with a mean vertical static offset of 30.5 ± 13.6 cm and overall values ranging between 5.3 and 54.1 cm. For the horizontal components, Kamai et al. (2014) predicted a static offset of 21.0 cm (11.5–38.2 cm with the related standard deviation), while we estimate a lower mean value (10.3 ± 4.5 cm), with our simulations varying overall between 4.4 and 14.8 cm. These comparisons show that our simulations give results comparable to the selected empirical predictions; further, our static offset values show a large scattering that reproduces the variability of the empirical relationships represented by the associated standard deviation.

6 Discussion

The results obtained for the ALCS220 and ALCS230 faults indicate that the ALCS220 fault is the most dangerous seismic source for the FCS site, not only because it has a slightly larger maximum magnitude (seismic moment and energy) but also due to the different geometry and setting of the fault. The upper edge of the fault is shallower than that of the ALCS230 fault, and if the rupture propagates from south to north, significant directivity effects could be observed, resulting in increased seismic motions. Thus, the horizontal ground motion estimated for ALCS220 is almost double, at least compared with the shaking related to the ALCS230 scenarios; the difference is less relevant on the vertical component, where the results are comparable within long periods for the displacement parameter (PGD and SO). For both scenarios, the PGM values on the east–west component are low, while the north–south and vertical components show the largest peak values.

The ALCS230 simulations show the strongest shaking in the vertical component for the receivers close to the source, while the ALCS220 fault generates its most significant ground motion values along the north component. This behaviour is related to both source and propagation effects.

The estimation of the near-field vertical component and its role in seismic hazard represents a challenging task for the seismological community. Shrestha (2009) observed that the 2/3 vertical/horizontal ratio (Newmark 1973) is not realistic in the near field, and this ratio can become larger than 1.0 within a 5 km source radius and larger than 2/3 within a 25 km radius. Thus, quite often, the vertical component of ground motion can be underestimated in the structural analysis of infrastructures and facilities located in the near-source area. Such behaviour could be deemed critical at the compressor station because resonance phenomena are generated when oscillation frequencies overlap with the building's characteristic content, and it may induce a non-conservative estimation of the seismic capacity (Di Michele et al. 2020). Nevertheless, the energy content of the vertical component of strong motion originating from a thrust fault can severely affect the maximum variation in the axial force in the building piers; indeed, if the vertical ground motion generates significant axial force variations in the piers, it influences the related shear and flexural capacities, which may drastically decrease (Di Michele et al. 2020).

In the near field, the modelled ground motion variability is larger than that obtained when considering GMMs because, within the near-source distance, very few data are recorded. For this reason, it is of critical importance to develop a DSHA, mainly of the near field, particularly for the vertical component, for the design of critical facilities. The uncertainties in the average ground motion parameters are high and comparable with the associated mean. Various input models generate very different shaking patterns, with the corresponding PGM values scattered within a broad range. This scatter confirms that the definition of the principal fault parameters, such as geometry and position close to the compressor station, is fundamental in estimating realistic seismic input. The shaking computed at the receiver grid depends on the different seismic radiation patterns generated by combining slip distributions and earthquake nucleation positions.

The source modelling adopted in this study and the resulting scenario are reasonable for the investigated area if we compare it with the effects of the MW = 6.4 earthquake that struck north-western Albania on November 26, 2019 (https://bbnet.gein.noa.gr) and caused the deaths of 51 people and injured approximately 3000. The maximum peak values were recorded at Durrës station (80 km north from Fier), located on soft soil approximately 15 km from the 2019 Albanian mainshock epicentre. The maximum PGA value recorded at Durrës was 192 cm/s2, the PGV was 38 cm/s, the PGD was 14 cm and the maximum SO was 5 cm (Duni and Theodoulidis 2020); it is worth mentioning that the signal stopped after 15 s due to electrical supply problems. Even if the magnitude of this event was slightly lower than those assumed in our scenarios, source distances and soil conditions are quite similar, and the simulated average north component values of PGV and PGD for the ALCS230 fault rupture scenario (Table 2) are consistent with the peaks recorded by the Durrës station.

The simulations produced ground motion parameters comparable to the shaking values observed worldwide in recent years; strong ground motion values were recorded in the vertical component, especially when thrust events (MW ≥ 6.0) occurred (Suzuki and Iervolino 2017), with a vertical PGA exceeding 1 g, while the PGV values exceeded 50 cm/s. Pacor et al. (2018) observed that these strong values do not seem to exhibit a clear dependence on magnitude and distance, but they may be related to the complexity of the rupture process, such as the localized failure of different portions along the fault. In this study, the rupture complexity on the finite fault is modelled using the pseudodynamic approach, which is strictly related to the slip distributions utilized in the computation; reproducing the observed near-field characteristics is a clear goodness marker for our simulation procedure.

Comparison with the recorded PGD and SO values extracted from available datasets can be crucial. Past near-field observations suggest that the vertical component is larger than the horizontal component at high frequencies, whereas the trend is reversed at lower frequencies (Bozorgnia et al. 1995). However, D’Amico et al. (2019) demonstrated that standard processing could underestimate PGD and PGV values in the case of strong motion accelerograms affected by static offset effects. For this purpose, they processed the waveforms recorded during the 2016 Norcia earthquake (MW = 6.5) by adopting an ad hoc tuned procedure; the new PGMs are comparable with our simulations, both for PGD and the static offset, on the overall ground motion components.

We also compared the DSHA at the compressor station with the ground motion obtained by applying the PSHA approach (Slejko et al. 2021). PSHA estimates the ground motion at different return periods by considering the ensemble of faults that could affect the site under investigation. In contrast, DSHA returns the shaking parameters from 15 different source models related to the most hazardous ALCS220 fault. Unlike GMMs (Fig. 9), PSHA is computed up to a cut-off period of 2 s, and the PSHA/DSHA match is limited to lower periods. In Fier, the average DSHA PGA is 514 cm/s2 (north component at soil conditions, Table 2), compared to a PSHA PGA of 525 cm/s2; this latter value refers to a 475-year return period. If we consider one standard deviation, the calculated PGA of 904 cm/s2 is similar to the PSHA value for a return period of 2475 years (893 cm/s2). The average DSHA PGV value of 73 cm/s (north component at soil conditions, Table 2) is also comparable with the PSHA PGV value of 70 cm/s computed for a return period of 2475 years, while the one-standard deviation peaks are comparable but for a more extended PSHA return period (9975 years). The PGV values are higher, mainly because the soil transfer function estimated for the Fier site amplifies the intermediate frequency content (approximately 1 Hz), which is usually the band in which PGV is more pronounced, and the directivity effect is enhanced. PGV is frequently adopted to quantify the dynamic effect of seismic action, and it is also used in damage functions (e.g., Pineda-Porras 2010). Therefore, the PGV values estimated in the DSHA analysis are crucial for the seismic design of the pipeline; high correlations have been observed between ground motion parameters and pipeline damages caused by earthquakes, as discussed by Bommer and Alarcon (2006). Furthermore, PGV is a critical parameter for evaluating the liquefaction potential assessment.

Regarding the acceleration response spectra (Fig. 10), the DSHA mean spectrum is comparable for T < 1.0 s with the PSHA spectrum for a return period of 475 years. One standard deviation DSHA spectrum values are within the PSHA spectra calculated for the 975- and 2475-year return periods. If we consider two standard deviations, the DSHA spectral values are within the PSHA spectra for return periods of 2475 and 9975 years. We also extract the five source models of the ALCS220 fault that generated the worst DSHA scenarios in terms of ground motion parameters, and we compare the mean spectral values with the PSHA results. As shown in Fig. 10, the DSHA spectral values are within the PSHA spectral values for the 975- and 2475-year return periods. If we add one standard deviation, our results are comparable with the spectra for a 9975-year return period. At longer periods (T > 1.5 s), the DSHA approach estimates higher values than those calculated with PSHA at very long return periods; these DSHA results can be of significant importance for the design of buried pipelines and the related infrastructure at the FCS site.

Fig. 10
figure 10

Comparison of PSHA spectra computed at various return periods for FCS with the mean ALCS220 spectra of the largest horizontal component (red line) and mean of 5 worst scenarios (blue line)

PSHA and DSHA are two different approaches that provide independent results. While PSHA gives an exceedance probability of a ground motion threshold within a selected return period, DSHA generates potential maximum values without considering any specific time interval. However, a qualitative comparison between PSHA and DSHA can help evaluate how much the ‘worst’ DSHA scenario is reproduced within the PSHA results and how much the DSHA simulations, including near-field effects, differ from the PSHA estimations.

7 Conclusions

Since the management of critical facilities requires detailed and reliable seismic hazard assessments, in this paper, we computed deterministic broadband scenarios for the TAP compressor station located in Fier (western Albanian coast). Previous analyses revealed that the most dangerous faults for the compressor station are the Vlore backthrust and Fier thrust. The ground motion is computed at bedrock conditions and on soft soils by applying specific transfer functions estimated from the available geotechnical data for the Fier area (Allkja and Dhimitri 2012). The site conditions at the compressor station amplify the ground motion, particularly in the intermediate frequency band generally used in estimating the dynamic effects and damage functions on pipelines. However, since the compressor station site is characterized by very soft sediments (sands and silty clay), further investigations on possible liquefaction phenomena are needed in the future to estimate their potential impact on facilities. Indeed, the synergistic role of extreme ground displacements and soil liquefaction can deeply affect gas pipeline operativity.

The position and geometry of the two considered faults and possible directivity effects play a crucial role in the estimation of the seismic input. For the near field, the vertical component is larger than the horizontal component, also in the low-frequency content, likely as an effect of the fault reverse kinematics. Accordingly, strong vertical ground motion can severely impact the seismic design of the compressor station, particularly the response to building structures, leading to a nonconservative estimation of the seismic capacity of such facilities. If these seismic effects are not considered in seismic design (Grimaz and Malisan 2014), seismic safety deficiencies can affect the near-field domain. Although it is known that the vertical component can be larger than the horizontal component in the high-frequency content (e.g., PGA), we reproduce strong vertical values at low frequencies (e.g., PGD).

When comparing, in a qualitative way, the present DSHA with the PSHA (Slejko et al. 2021) results for the FCS compressor station, we observe that the DSHA mean spectrum is comparable for T < 1.0 s with the PSHA spectrum for a return period of 475 years, while if we add one standard deviation, the DSHA average response spectra fall within the PSHA spectra calculated for 975- and 2475-year return periods. The highest values obtained with the DSHA analyses are comparable with the PSHA results at long periods (T > 1.5 s) that are critical to the seismic design of buried pipelines.

This study also describes a pioneering approach to estimating, in such detail, the seismic hazard assessment for a compressor station, and can be adopted for other pipelines crossing active faults.