INTRODUCTION

Infrastructures that extend for kilometers and that have to last for some decades, as is the case of pipelines, can be affected in their integrity by differential movements of the ground (e.g., O’Rourke and Liu, 1999, and references therein). Earthquakes can determine sudden and relevant surface displacements so that the pipeline design in critical frames must consider the possible displacements in correspondence of major active or possibly active faults, where sudden permanent deformation could more easily occur. According to Caputo (2005), only large-magnitude earthquakes are capable of generating or modifying the surface morphology instantaneously and permanently (i.e., are linear morphogenic earthquakes). A classic example is the case of the 2002 Denali (Alaska) earthquake, moment magnitude (M_w) &equals; 7.9 (Haeussler et al., 2004). Among some recent events close to the pipeline route are the cases of the 1978 Thessaloniki and 1995 Kozani earthquakes (Caputo et al., 2008). Indeed, earthquake-induced shear, compression, or tension in the ground can strongly affect such buried infrastructures depending on the relative orientation of the pipeline and the fault (Ballantyne, 2008). The worst case occurs when the pipe is affected by pure compression, whereas in the case of tension, steel pipes with welded joints can distribute lengthening over hundreds of meters, thus minimizing localized stresses. According to O’Rourke and Liu (1999), pipeline breakage can occur for relatively small ground displacements (i.e., <12.7 cm [<5 in.]). Empirical formulas such as the ones of Wells and Coppersmith (1994), Ambraseys and Jackson (1998), or S. Pavlides and Caputo (2004) may be used to estimate the surface displacement for events of different magnitudes and fault mechanism.

Thanks to the global navigation satellite system (GNSS), we know that tectonic plates and crustal blocks are in an almost continuous relative movement, most pronounced in the narrow zones between plates. By repeatedly measuring the distances between particular points, we can determine the movement along faults or between plates, the rate, and the likely areas of strain accumulation (Lisowski et al., 1990; R. E. Reilinger et al., 1997; Mantovani et al., 2001; Wright et al., 2001; Cenni et al., 2012; Caporali et al., 2018). Plate boundaries move with variable velocities, from a few millimeters to some centimeters per year (from hundredth to tenth of an inch), implying significant total displacements in 50 yr (e.g., Altamimi et al., 2012). Even without the occurrence of strong earthquakes, a pipeline crossing active plate boundaries must cope with remarkable differential motions along its length because of the interseismic elastic strain accumulation within the upper crust. Such movements, leading to permanent ground deformation (PGD), can distress the pipe, whereas the anchor points can result in local stress concentrations.

The Southern Gas Corridor (SGC) is a significant infrastructure project for the construction of a pipeline that will transport natural gas coming from the fields of Shah Deniz II in Azerbaijan to the Italian Coast for distribution in central and southeastern Europe. The final step of the SGC project (in the following, SGC–Trans Adriatic Pipeline [TAP]), starting at the Turkey–Greece border, running through Greece and Albania, and, after passing beneath the Adriatic Sea, ending in Southern Italy, will be occupied by the TAP. The total length of the pipeline will be 870 km (540.6 mi), of which 545 km (338.6 mi) run in Greece, 211 km (131.1 mi) in Albania, 105 km (65.2 mi) across the Adriatic Sea, and 8 km (5 mi) in Italy.

The SGC–TAP crosses one of the most active zones of the Mediterranean, with many active faults (Caputo et al., 2012; Caputo and Pavlides, 2013; Ganas et al., 2013) and intense seismic activity (e.g., Papazachos and Papazachou, 1997). Several studies previously performed in the area crossed by the pipeline and, based on GNSS observations, evidenced a complex pattern of deformation velocities (e.g., Altiner et al., 2006; Pérouse et al., 2012; Devoti et al., 2017). By collecting and homogenizing the measurements throughout the Calabrian Arc, Balkans, Aegean, and Anatolia peninsula, the modeled velocity field shows a first-order distribution, characterized by two opposite toroidal crustal patterns, clockwise in northwestern Greece and counterclockwise in western Anatolia (Pérouse et al., 2012). On a smaller scale, however, differential motions between blocks separated by major faults or distributed shear zones can also be recognized (Jouanne et al., 2012).

The present study focuses on predicting the areas where critical differential motions could be expected along the SGC–TAP route, over the nominal pipeline lifespan of 50 yr. Based on recent tests on pipelines of 760 mm (30 in.), O’Rourke and Liu (2012) hypothesize that a total displacement of only 76–102 mm (3–4 in.) could cause the tearing of the pipe wall. The analysis of various cast-iron pipeline breakages evidences how the normalized pipe break rate is a nonlinear function of the PGD, with a faster rate for displacements smaller than 12.7 cm (<5 in.). Most of the breakages of pipelines occur for values approximately 50 mm (∼1.97 in.), apparently overall for compressional faulting and vertical displacements (Porter et al., 1991; O’Rourke and Liu, 1999). Pour-Ghaz et al. (2018) indicate that for segmental plain– and steel fiber–reinforced concrete pipelines, subjected to permanent ground displacement, the damage cracking in the immediate vicinity of the fault line begins with as little as 50 mm (1.97 in.) of displacement.

To identify the areas of the SGC–TAP where critical differential motions can most likely occur, we followed two approaches. The first one (continuum approach) considers a continuous crust not affected by faults; the second one (blocky approach) assumes a crust consisting of separate blocks, each one characterized by a relatively uniform tectonic style and behavior. We carefully analyzed the available GNSS data and compared the geodetic results to the pattern of the most significant mapped faults affecting the area (e.g., Caputo and Pavlides, 2013) for assessing the tectonic style and behavior along the route. Then, we processed the GNSS velocity data. In the continuum approach, the sparse, available strain rate information is interpolated along the pipeline route with a linear approach, considering the region as a continuum. In the blocky approach, we instead applied an original algorithm that allows the linear interpolation within individual blocks, characterized by a relatively homogenous strain velocity, or a specific tectonic setting, independently upon the neighboring ones. One of our aims is verifying the efficiency of an interpolation, driven by the knowledge of tectonic domains, and the local consistency of the velocity solutions. The results of the two methods are finally compared by calculating the maximum displacement that would cumulate in the next 50 yr of the pipeline lifespan and, for the blocky approach, the differential displacements that could cause possible bending phenomena to the pipeline structure.

REGIONAL FRAMEWORK

Geodynamics of the Area

From a geodynamic viewpoint, the broader investigated area spans the Apulian plate, crosses the Hellenides fold–thrust belt, and reaches, to the east, the southern Balkans–Northern Aegean extensional province. Apulia is often considered as a promontory of the African (Nubian) plate, and it is currently moving to the north-northeast (e.g., Devoti et al., 2011). The Hellenic orogen is the result of the polyphase convergence between Gondwana and Eurasia, locally oriented southwest-northeast, which has been active from the Mesozoic. The process of oceanic subduction of the Nubian plate beneath the Eurasian southern margin has led to a nappe stacking of the several interposed microcontinents and oceanic basins characterizing the Neo-Tethyan ocean (e.g., Biju-Duval et al., 1977; &Scedil;engör, 1985; Robertson et al., 1991; Papanikolaou, 2013). As a result, tectonic transport and imbrication of structures to the south-southwest has occurred since the Early Jurassic, with the progressive involvement of more external terranes and associated development and/or rejuvenation of contractional, northwest-southeast–striking structures (e.g., Jacobshagen, 1986; Mountrakis, 2006; Papanikolaou, 2013).

Along the Albanian and the northwestern Greek sectors of the convergence front (northwest of the Cephalonia transform fault; Figure 1), the subduction of the oceanic lithosphere ended in the late Eocene–early Oligocene (Royden and Papanikolaou, 2011). Slow continental collision, that is, 5 mm/yr (0.2 in./yr) as a long-term average, has been occurring until the present. However, in the back-arc region, widespread extension began since the late Eocene–early Oligocene for the Internal Hellenides (Circum Rhodope and Axios-Vardar zones) and since early Miocene for the External Hellenides (Pelagonian, Pindos, Ionian; Mercier et al., 1979a, 1989; Mountrakis, 2006). In the beginning, the crustal stretching was oriented east-northeast–west-southwest, therefore suggesting the occurrence of a postorogenic collapse (Mercier, 1976; Mercier et al., 1989; Caputo and Pavlides, 1993) or gravitational spreading (Le Pichon and Angelier, 1979). However, the onset of the current stress field for the Aegean Sea region and the central-eastern part of northern Greece characterized by an almost north-south maximum tensile direction, with a slight deflection of the trajectories toward north-northwest–south-southeast in northwestern Greece, can be traced back to the Pleistocene (late Pliocene).

Figure 1. Horizontal global positioning system velocities (gray arrows) and corresponding error ellipses (gray ellipses) used in the present research. The Southern Gas Corridor–Trans Adriatic Pipeline route is in green; primary seismogenic structures affecting the broader investigated area are simplified as red lines. Details on the Nubia–Eurasia rotation pole are in Pérouse et al. (2012). CTF &equals; Cephalonia transform fault; FYROM &equals; former Yugoslav Republic of Macedonia.

An additional element for the comprehension of the tectonics of the northern Aegean Region is the propagation of the dextral strike-slip North Anatolian fault zone (McKenzie, 1972) into the Aegean Sea, where it generates the North Aegean trough (see Koukouvelas and Aydin, 2002 for a review on the geometry and architecture of the North Aegean trough). Such a dextral shear motion across the Aegean Region likely started in the early Pliocene (Armijo et al., 1996), progressively propagating westward.

The current geodynamic setting of the area crossed by the SGC–TAP is characterized by tectonic elements belonging to either compressional, extensional, or transcurrent regimes (e.g., Caputo et al., 2012; Figure 1). More specifically, the western sector around the Albanian shoreline and the southern Adriatic Sea is characterized by both offshore and onshore east-northeast–dipping low-to-medium dip angle active thrusts (Figure 2). Along the highlands of the Albanides, a sudden change of the stress field occurs in terms of a likely swap between the stress tensor components (maximum [σ₁] and minimum [σ₃] principal stresses), and, as a result, dip-slip kinematics on north-south–striking normal faults is observed (Copley et al., 2009). Moving farther to the east into the western Macedonia region (Greece; Figure 1), the direction of the stress component σ₃ slightly rotates, becoming northwest-southeast (e.g., Mercier et al., 1989; Mountrakis, 2006), and the main tectonic structures are represented by southwest-northeast–striking normal faults (S. B. Pavlides and Mountrakis, 1987; Doutsos and Koukouvelas, 1998; Goldsworthy and Jackson, 2000; Doutsos et al., 2006; Caputo et al., 2012 and references therein). Farther to the east (Figure 3), the Thessaloniki and Mygdonian Graben region is subjected to a north-south extension on mainly east-west–striking dip-slip normal faults (e.g., Hatzfeld et al., 1987; Mountrakis et al., 1996a, b; Goldsworthy et al., 2002). Lastly, west-southwest–east-northeast– and east-west–striking normal– and oblique-slip faults characterize southern Thrace and the region close to the border with Turkey. The oblique-slip faults mainly exhibit dextral kinematics, thus suggesting the likely influence of the North Aegean trough on the regional stress field (e.g., Mountrakis et al., 2006; Caputo et al., 2012).

Figure 2. Zoom on the western part of the route showing the global positioning system velocities (black arrows) and relative errors (ellipses). The yellow polygons represent the horizontal projection of the seismogenic sources provided in Greek Database of Seismogenic Sources (Caputo and Pavlides, 2013), where the thick border indicates their upper tip projection (red, green, and blue correspond to a prevailing reverse, strike-slip, and normal kinematics, respectively). For the sake of graphical simplicity, the seismogenic sources cited in the text are labeled with only the numerical part, whereas the country code (AL and GR) and the source code (CS for composite source or IS for individual source) are omitted. Location map in Figure 1. The purple arrows indicate the Ohrid Graben and the Ptolemaida Graben, cited in the text. Details on the Nubia–Eurasia rotation pole are in Pérouse et al. (2012). Black labels are global navigation satellite system station names. AG &equals; from D’Agostino et al. (2008, 2011); FL &equals; from Floyd et al. (2010); FYROM &equals; former Yugoslav Republic of Macedonia; M &equals; from Kotzev et al. (2006); RE &equals; from Reilinger et al. (2006).

Figure 3. Zoom on the eastern part of the route showing the global positioning system velocities (black arrows) and relative errors (ellipses). Location map in Figure 1. The yellow polygons represent the horizontal projection of the seismogenic sources provided in Greek Database of Seismogenic Sources (Caputo and Pavlides, 2013), where the thick border indicates their upper tip projection (red, green, and blue correspond to a prevailing reverse, strike-slip, and normal kinematics, respectively). For the sake of graphical simplicity, the seismogenic sources cited in the text are labeled with only the numerical part, whereas the country code (AL and GR) and the source code (CS for composite source or IS for individual source) are omitted. Black labels are global navigation satellite system station names. FYROM &equals; former Yugoslav Republic of Macedonia; MA &equals; from Matev (2011); NAFZ &equals; North Anatolian fault zone; RE &equals; from Reilinger et al. (2006).

During the last decade, several sectors of the investigated area were monitored by both permanent GNSS networks and repeated campaigns on temporary stations. The resulting velocity data constitute the input for geodynamics and regional strain field models.

Along the Adriatic Sea, we record 3–5 ± 1 mm/yr (0.1–0.2 ± 0.03 in./yr), in the northeast direction in the Apulia area (Devoti et al., 2011). High velocity (8–10 mm/yr [0.3–0.4 in./yr]) is also recorded in some islands in central Adria between the Gargano zone and the central Dinarides, with northeast-oriented vectors (Altiner et al., 2006). Across the Adriatic Sea, the subduction of the Adria microplate under the Dinarides–Albanides causes a velocity drop of 8.7 ± 0.5 mm/yr (0.3 ± 0.02 in./yr) in 200 km (124.3 mi). On the eastern side of the Adriatic Sea, the velocity field is characterized by a south-westward motion of the external Albanides increasing from north to south and by a southward velocity of the inner Albanides, former Yugoslav Republic of Macedonia (FYROM), Bulgaria, and northern Greece increasing from north to south (Jouanne et al., 2012; Figure 1). The major dextral strike-slip faults in Albania are associated with an inversion of the velocity parallel to the coastline of 7.1 ± 0.5 mm/yr (0.3 ± 0.02 in./yr) occurring in 350 km (217.5 mi) (Caporali et al., 2009).

The subduction along the Hellenic Arc dominates the tectonic deformation in the study area (Papazachos, 1988; Wortel et al., 1990), with a 30.5 mm/yr (1.2 in./yr) mean southwestward motion roughly perpendicular to the west Hellenic trench. We record 20.6 ± 0.8 mm/yr (0.8 ± 0.03 in./yr) in eastern and central Turkey, 24.6 ± 1.0 mm/yr (9.7 ± 0.04 in./yr) in western Turkey, and up to 31.1 ± 0.9 mm/yr (1.2 ± 0.04 in./yr) in the central and southern Aegean (R. Reilinger et al., 2006; Devoti et al., 2017). Eastern Macedonia and Thrace in northern Greece are affected by the stretching of the upper crust, with a relative velocity up to 3.2 ± 0.9 mm/yr (0.1 ± 0.04 in./yr) in 440 km (273.4 mi) (Caporali et al., 2009).

Seismogenic Structures

In the following section, we briefly analyze the primary active fault zones close to and crossing the pipeline path, moving from the western sector offshore Albania to the eastern termination of the investigated area at the Greek–Turkish border (Figure 1). The principal seismotectonic characteristics, as well as the labels of the structures (Figures 2, 3), are from the Greek Database of Seismogenic Sources (GreDaSS) (Caputo and Pavlides, 2013; http://gredass.unife.it), to which the reader is referred to for the relevant references.

The first tectonic structure crossed by the pipeline path is the Karaburuni–Sezani detachment (Albania Composite Source [ALCS]201 in GreDaSS; Figure 2) described as an east-northeast–dipping, low-angle, reverse fault. Although the maximum expected earthquake has been estimated to have an M_w greater than 7.0 (Caputo and Pavlides, 2013), the scarce instrumental and historical seismicity available for this structure leaves some uncertainties on the possible occurrence of internal segmentation, as well as on the slip rate and recurrence time of each segment. Once on land, the route crosses other reverse structures such as the Ardenica back thrust (ALCS235) and the Osum thrust (ALCS240). The former is a secondary fault with a small seismogenic potential, whereas the possible activity of the latter as an out-of-sequence structure is still debated, although the faults could be capable of generating moderate events. Because of their overall dimensions, setting, kinematics, and location, it is doubtful that the associated earthquakes could be linear morphogenic ones (Caputo, 2005).

It is noteworthy that the Albanides convergence front is segmented along its roughly north-south–trending strike by some significant transverse faults, with prevailing transcurrent kinematics and transfer role, such as the Shkodra (ALCS390), Lushnje (ALCS290), and Qeparo Transfer faults (ALCS085) (Figure 2). Regarding the SGC–TAP, the closest among these structures is the Lushnje Transfer fault, striking east-northeast–west-southwest, steeply dipping SSE and characterized by a dextral sense of shear. The maximum magnitude M_w associated with this structure is approximately 6.3 (Caputo and Pavlides, 2013), and the minimum distance from the fault trace is approximately 10 km (∼6.2 mi).

No active faults have been recognized east of Berat (near the GNSS station, BERA, in Figure 2) along the pipeline route and for more than 50 km (>31.1 mi) across the high-topography internal region of southern Albania. Possible stability problems for the pipeline are only and basically of geotechnical concern. Toward the border with the FYROM and Greece, the most relevant active tectonic structures are represented by an 80-km-long (47.7-mi-long) fault system, forming the roughly north-south Ohrid Graben, which locally accommodates east-west crustal extension (Figure 2). The pipeline crosses its southern sector and, in particular, the normal faults named Korçë (ALCS515), Podgoritsa (ALCS510), Maliqi (ALCS505), and Bilisht (ALCS520) (Figure 2), whose strike ranges between N330°W and N30°W. The prevailing dip-slip kinematics is documented by structural data and focal mechanisms from moderate magnitude earthquakes (e.g., Muço, 1994; Kiratzi, 2011). A recent confirmation comes from an M_w &equals; 5.1 magnitude event that occurred in June 2019 (https://bbnet.gein.noa.gr/). According to their dimensions and geodynamic setting, all these structures could release moderate earthquakes; however, none of them could likely generate significant direct coseismic ground ruptures at intersection points with the pipeline.

Reaching the Greek territory, the tectonic regime changes once more as far as the Kastoria dip-slip normal fault (Greece Composite Source [GRCS]090) has an east-west orientation, thus documenting a roughly meridian crustal stretching direction. At the same time, eastward, the pipeline enters and crosses all along the northeast-southwest–trending Ptolemaida Graben (Figure 2). This structure is bordered by the Amyndeo (GRCS070), Ptolemaida (GRCS072), and Komanos (GRCS077) composite sources (Figure 2) as well as minor individual sources with faults dipping both northwest and southeast (Caputo and Pavlides, 2013). Their kinematics is mainly dip slip, confirming the change (relative to the western Albanian sector) in the Quaternary stress trajectories that characterizes the broader Aegean Region (e.g., Mercier et al., 1979a, 1989). Based on geological estimates and geodetic data, the slip rates of the Ptolemaida Graben faults are in the order of 0.1–0.3 mm/yr (3.9 × 10⁻³–1.2 × 10⁻² in./yr), suggesting prolonged recurrence times. In contrast, their estimated seismogenic potential in the case of multiple-segment ruptures could be as high as 6.6–6.9 for the M_w (Caputo and Pavlides, 2013).

Moving farther to the east (Figure 3), the stress trajectories and consequently the orientation of the major crustal faults slightly change (10°–20° clockwise rotation), as documented by the South Almopia fault (GRCS068), which is characterized by prevailing dip-slip kinematics with a minor dextral component of motion. The progressive change of the stress trajectories also continues in central Macedonia. The Anthemoundas (GRCS250), Asvestochori (GRCS245), and Mygdonia (GRCS100) composite sources (Figure 3) show dip-slip kinematics, north dipping, and east-west to east-southeast–west-northwest trending. The slip rates of these structures are in the order 0.1–0.4 and 0.1–0.7 mm/yr (3.9 × 10⁻³–1.6 × 10⁻² and 3.9 × 10⁻³–2.8 × 10⁻² in./yr) (e.g., Tranos et al., 2003), and the maximum possible earthquake could be as strong as M_w &equals; 7 (Caputo and Pavlides, 2013). In this regard, the Gerakarou fault (the central segment of the Mygdonia composite source) was responsible for the M_w 6.4 Thessaloniki earthquake of 1978 (Mercier et al., 1979b). It should be noted that north of Thessaloniki, the SGC–TAP turns northeastward; therefore, the infrastructure crosses only the westernmost part of the northern Chalkidiki fault system, where coseismic ground ruptures are likely to be limited, even in case of a future linear morphogenic earthquake.

North of the Chalkidiki peninsula, the pipeline turns eastward, running roughly parallel to the Serres fault (GRCS145; Figure 3) for approximately 50 km (∼31.1 mi) all along its hanging-wall block. This structure dips south- and southwestward as far as it is clearly segmented (Mountrakis et al., 2006). Estimated slip rates are very low (0.05–0.1 mm/yr [1.9 × 10⁻³–3.9 × 10⁻³in./yr]), and no historical earthquakes are reported (Papazachos and Papazachou, 1997). However, in the case of a reactivation of the entire source, the maximum magnitude predicted by empirical relationships should be approximately M_w &equals; 6.7 (Caputo and Pavlides, 2013).

More eastward, the pipeline runs just south of the Drama fault (GRCS140), and then along the hanging-wall block of a significant seismogenic source (Thrace fault zone; GRCS150; Figure 3). It consists of an approximately 120-km-long (∼74.6-mi) structure bordering the Rhodope Massif to the north and the Kavala–Xanthi–Komotini basin to the south (Mountrakis and Tranos, 2004). This composite source is divided into four major segments, but only the three eastern ones show evidence of recent activity (Xanthi [GRCS150a], Iasmos [GRCS150b], and Komotini [GRCS150c] faults; Figure 3). The fault barriers are of an angular type producing marked differences in strike and dip (basically southward) and generating a zigzag geometry. The seismicity is quite low, although long-term slip rates are in the range of 0.2–0.5 mm/yr (7.9 × 10⁻³–2.0 × 10⁻³in./yr). Although the maximum recorded earthquake was in 1784 (M_w &equals; 6.7) (Papazachos and Papazachou, 1997), if we assume the contemporaneous rupture of multiple segments, the estimated maximum magnitude could be even greater than 7 (Caputo and Pavlides, 2013).

In the very final sector of the investigated area, before crossing the border with Turkey, the SGC–TAP runs close to and is potentially affected by the Maronia fault (GRCS160). Such a fault zone strikes almost east-west dipping southward and exhibits slightly oblique-slip kinematics being possibly influenced by the transtensional stress field associated with the nearby North Anatolian fault zone–North Aegean trough system (Figure 3). The minor dextral strike-slip motion component (rake ∼250°–260°) is also confirmed by the analysis of the focal mechanisms (Kiratzi et al., 2005). The associated slip rate is in the range of 0.4–0.7 mm/yr (1.6 × 10⁻³–2.8 × 10⁻² in./yr), and the maximum M_w is approximately 7.0 (Chatzipetros and Pavlides, 2004).

DATA ANALYSIS

One of the primary needs for a study like the present one, involving the comparison, in modulus and directions, of velocities along the entire SGC–TAP, is to deal with a velocity field as homogeneous as possible, reducing the systematic errors arising from different reference frame adoptions. In this regard, Pérouse et al. (2012) performed a considerable work of revision and homogenization of the global positioning system (GPS) velocity data acquired and processed in different sectors of the central and eastern Mediterranean to present a new kinematic and strain model of this broad region. Thanks to their work, a catalog of the displacement velocities, with relative errors, is available and homogeneously calculated compared to a fixed Eurasia (Figure 1). Details on the Nubia–Eurasia rotation pole are in Pérouse et al. (2012). Devoti et al. (2017) made an analogous effort, but the coverage in the area of our interest is much lower.

Based on the catalog of Pérouse et al. (2012), Figure 4 shows a synthesis of the velocities observed in the broader area crossed by the SGC–TAP route, showing two distinct sectors: western and eastern. We averaged the values of the velocities with the same orientation: the color is proportional to the velocity. The western sector, across the southern Adriatic Sea–Albania (Figure 4, ellipse A), is characterized by interplate convergence between Apulia and the External Hellenides–Albanides, as clearly suggested by the opposite velocity vectors on the two coastal areas of the Adriatic Sea (Figures 2, 4). In contrast, the central-eastern sector is affected by a southward positive velocity gradient progressively increasing from 4 to 8 mm/yr (0.2–0.3 in./yr), thus suggesting the internal deformation of the crust, which is undergoing north-south stretching (Figure 4, ellipse B). Farther to the east, the velocities are higher and reflect the counterclockwise rotation of the Anatolian microplate (Figure 4).

Figure 4. Representative velocity vectors (colored arrows) relative to a fixed Eurasia frame (simplified from Pérouse et al., 2012) in the broader area crossed by the investigated sector of the Southern Gas Corridor–Trans Adriatic Pipeline (green line). The arrow color scale is indicative of the magnitude of the velocities. The gray ellipses mark the area characterized by interplate convergence between Apulia and the External Hellenides (A) and a part of the broader southern Balkan–Northern Aegean Region (B) affected by north-south crustal stretching. Details on the Nubia–Eurasia rotation pole are in Pérouse et al. (2012).

The Data

From the catalogs of the GPS measurements provided by Pérouse et al. (2012 and references therein), we selected the data in a buffer zone around the SGC–TAP (Figures 1–3). We gave priority to the data nearest to the pipeline (GPS sites within a distance of ∼120 km [∼74.6 mi]) and/or with the smallest error ellipses (major axis smaller than 6.8 mm [<0.3 in.]). For each station, the accuracy of the measurement was evaluated, considering whether it is a permanent or temporary station, the length of the time series for the permanent sites, or, in the case of a temporary station, the number of repeated measures.

Table 1 summarizes all relevant information used in our analyses. Regarding the station labels, in general, we used the official published ones; in other cases, for simplicity, stations are labeled based on the original papers and numbered sequentially. It is the case of the stations, that the labels indicate information comes from D’Agostino et al. (2008, 2011) AG01, AG02, and AG03 or from Matev (2011) MA01, MA02, and so on. Table 1 also reports the name in the original publication.

Quantitative Analyses

To give a quantitative estimate of the expected displacements along the whole route of the pipeline, we performed a three-dimensional interpolation on the data over equally spaced points (Figure 5). Our scattered data set is defined by locations, (x, y) or (longitude, latitude), and by the corresponding values d. The Delaunay triangulation function of (x, y), yielding the surface F(x, y) &equals; d, has been selected to interpolate the data. The function F has been used to estimate unknown values of d_i at locations (x_i, y_i) included in the original (x, y) domain (de Berg et al., 2008; Weisstein, 2010).

Figure 5. (A) Horizontal velocities (blue arrows) and corresponding estimated errors (gray ellipses) interpolated at closely and equally spaced points along the Southern Gas Corridor–Trans Adriatic Pipeline route (green line) based on the continuum approach. (B) The same map with the original data (gray arrows and ellipses) for reference. Black labels are global navigation satellite system station names. Additional labels are informational references: AG &equals; from D’Agostino et al. (2008, 2011); FL &equals; from Floyd et al. (2010); M &equals; Kotzev et al. (2006); MA &equals; Matev (2011); RE &equals; from Reilinger et al. (2006).

The five-dimensional F interpolating function has been inferred, using a linear method of interpolation, from the following values (or d, as mentioned above): velocity in the east-west direction, velocity in the north-south direction, east-west error, north-south error, and the correlation between the two components.

After that, we used the F function to evaluate the corresponding values d over 49 regularly spaced (20 km [12.4 mi]) points along the SGC–TAP route. We followed two approaches: (1) a continuum approach along the whole path and (2) a blocky approach, in which the F functions are calculated within distinct crustal blocks, showing an independent behavior. As mentioned above, such blocks are characterized by roughly uniform velocity trends markedly delimited and/or associated with different faults (Figure 6). With both approaches, we linearly interpolated all the available measurements within a distance of 120 km (74.6 mi) from the SGC–TAP route, weighted for the distance from the path.

Figure 6. Gray dashed lines represent the boundaries between the blocks assumed to have relatively homogenous geodetic behavior and used in the blocky approach. Each block is indicated with a black number in a circle. Gray arrows indicate the global positioning system velocities with their relative errors ( gray ellipses). Red lines represent the simplified upper tip of the seismogenic sources (Figure 1). Black labels are global navigation satellite system station names. Additional labels are informational references: AG &equals; from D’Agostino et al. (2008, 2011); FL &equals; from Floyd et al. (2010); M &equals; Kotzev et al. (2006); MA &equals; Matev (2011); RE &equals; from Reilinger et al. (2006).

Based on the velocities inferred employing the two approaches, we also calculated the cumulative displacements, with the relative uncertainties, expected in the 50 yr representing the average life of this kind of energy infrastructure.

Continuum Approach

Figure 5A shows the results of the data continuum interpolation. The horizontal velocities reconstructed by the interpolating algorithm and the relative errors have a good match with the observations (Figure 5B). The GNSS station, KAST, with a minimal error ellipse and very near to the SGC–TAP, ensures smaller error ellipses in the central part.

Shown in Table 2 are the velocity values obtained, with the relative errors and displacements expected after 50 yr, with the corresponding uncertainty. As mentioned above, all velocities are relative to fixed Eurasia (details on the Nubia–Eurasia rotation pole are in Pérouse et al., 2012). The minimum velocity and, hence, the expected minimum displacement, is recorded at the 10th interpolation point, in Albania, Berat district, near the city of Poliçan, with a velocity value of 0.6 mm/yr (0.02 in./yr), implying a predicted total displacement in 50 yr of only 30 mm (1.2 in.). However, the maximum velocity and displacement are recorded at the 28th and 29th points of the interpolation in Greece near the city of Chalkidona, 15–30 km (9.3–18.6 mi) northwest of Thessaloniki and to the southeast of the South Almopia fault (GRCS068), with values of 7.1 mm/yr (0.3 in./yr), implying a cumulative displacement of 355 mm (14 in.) in 50 yr, with an uncertainty of 55 mm (0.4 in.).

Blocky Approach

The blocky approach aims to highlight the existence of discontinuities in the deformation field. Distinct blocks and their relative boundaries (Figure 6) have been recognized by taking into account (1) the traces of the primary active faults (Caputo and Pavlides, 2013 and references therein) and the general geodynamic setting of the region; (2) in the absence of recognized tectonic structures, the relief characteristics that could justify a change in the observed horizontal velocities; and (3) the occurrence of homogeneous GPS clusters. The ensuing uncertainty in the position of the block boundaries, and, hence, of the differential displacements occurring, can be estimated to be in the order of 10–20 km (6.2–12.4 mi), comparable to the distance between two subsequent points.

Block 1 includes the Apulia and part of its offshore (Figure 2); block 2 incorporates the Karaburuni–Sezani detachment (ALCS201 in GreDaSS), the Ardenica back thrust (ALCS235), and the Shkodra (ALCS390), Lushnje (ALCS290), and Qeparo Transfer faults (ALCS085). At the same time, the Osum thrust (ALCS240) marks its boundary with block 3, which occupies the westernmost part of Albania (Figure 2). Block 4 includes the normal faults Korçë (ALCS515), Podgoritsa (ALCS510), Maliqi (ALCS505), and Bilisht (ALCS520) (Figure 2). The central part of the corridor falls almost entirely within block 5, characterized by a relatively homogeneous velocity field with roughly southward vectors. It is worth noting that the western sector of this block is affected by northeast-southwest–trending faults (Kastoria fault [GRCS090] and Ptolemaida Graben with the Amyndeo [GRCS070], Ptolemaida [GRCS072], and Komanos [GRCS077] composite sources) (Figure 2). In the remaining part of the same block, the major crustal faults are nearly east-west (i.e., South Almopia fault [GRCS068], Anthemoundas [GRCS250], Asvestochori [GRCS245], and Mygdonia [GRCS100] composite sources; Gerakarou fault [GRIS101]; and Serres fault [GRCS115]) (Figure 3). However, considering that all selected velocity vectors included in the block are almost parallel and the western part contains only two stations, we decided not to subdivide block 5 further. Block 6 includes the primary seismogenic source of the Thrace fault zone (GRCS150) as well as some other structures in southern Bulgaria. Finally, the easternmost block 7, including the region from Thassos to the Black Sea, is clearly influenced by the westward motion of the Anatolia microplate. It includes the Maronia fault (GRCS160) and the set of structures belonging to the North Aegean trough (Figure 3).

Based on the blocks shown in Figure 6, we group the stations into clusters, where we linearly interpolate the velocity. When the data coverage is not sufficient (generally toward the block boundaries), we fill the gaps with the last reliable value obtained. Hence, for the first two blocks, we use the nearest-neighbor interpolation algorithm so that we merely extend offshore the most reliable information we have on the coast. For the other blocks, the velocities are linearly interpolated within the block, weighting the data for their distance from the pipeline route. This procedure implies sudden changes at the boundaries of the blocks.

To safely interpolate the velocities and relative errors, we applied the same algorithm used in the continuum approach but using only the GNSS stations within every single block (indicated with different colors in Figure 7). In this way, where abrupt changes occur in the velocity distribution (Figure 7A), no intermediate orientations are artificially created in the interpolated points, and the pattern characterizing each separate block is better preserved. For facilitating a comparison between the two approaches, we applied the procedure to the same equally spaced points along the SGC–TAP route. Reported in Table 3 are the obtained velocity values (with the relative coordinates and uncertainties) as well as the estimated displacements after 50 yr (with their uncertainties). The minimum cumulative displacement is recorded at the 14th interpolated point in Albania in block 3 at the boundary between the Berat and the Korçë districts, near Lekas town. The maximum, instead, is recorded at the 28th interpolated point in Greece in block 5, 25 km (15.5 mi) to the northwest of Thessaloniki, to the southeast of the South Almopia fault (GRCS068) (i.e., nearly coinciding with the results of the continuum approach), with the same value of 7.1 mm/yr (0.3 in./yr), resulting in 355 mm (14 in.) in 50 yr, and uncertainty of 55 mm (0.4 in.).

Figure 7. (A) Horizontal velocities (colored arrows) and corresponding estimated errors (colored ellipses) interpolated at equally spaced points along the Southern Gas Corridor–Trans Adriatic Pipeline route (green line) based on the blocky approach (the points belonging to each of the blocks of Figure 6 are indicated with the same color: brown [block 1]; blue [block 2]; red [block 3]; dark green [block 4]; green [block 5]; magenta [block 6]; violet [block 7]). (B) The same map with the original data (gray arrows and ellipses) for reference. Black labels are global navigation satellite system station names. Additional labels are informational references: AG &equals; from D’Agostino et al. (2008, 2011); FL &equals; from Floyd et al. (2010); M &equals; Kotzev et al. (2006); MA &equals; Matev (2011); RE &equals; from Reilinger et al. (2006).

DISCUSSION

Comparing the information available on the major active faults affecting the region to the one obtained from the analysis of the GPS measurements could contribute to discriminating between faults that separate blocks with a nearly homogeneous behavior from other tectonic structures that do not affect the horizontal velocities among blocks. For this purpose, a geographic information system was utilized to analyze the velocity measurements in connection with the primary active faults affecting the investigated region (Caputo and Pavlides, 2013; Figure 1).

By examining the velocity distribution, it is possible to recognize areas characterized by a homogeneous pattern and other ones where abrupt velocity changes occur, both in module and direction, thus suggesting the kinematic effect of major active tectonic structures. For example, whereas Apulia is rigidly moving toward the north-northeast (AG01, AG02, AG03, AG14, and AG15 sites), the GPS stations on the Albanian shoreline (GNSS stations KRYE, APOL, and RRAD) point toward the north-northwest (Figure 2). The strain velocities in the internal region of southern Albania (GNSS stations LABO, DERV, and LESK) are perfectly parallel to the west-southwest direction of convergence, as documented by the slip vectors on the principal thrusts (Muço, 1994; Louvari et al., 2001; Serpetsidaki et al., 2012). This remarkable difference in the direction from the coastal to the internal sector is likely because of a strong mechanical coupling characterizing the Karaburuni–Sezani detachment (ALCS201) and the Durresi Offshore thrust (ALCS315). These structures represent the basal decollement of the accretionary wedge in this sector of the orogen. Therefore, the north-northeast motion of the underlying Adria plate is partly transmitted to the overlying wedge and, hence, up to the topographic surface (viz. geodetic stations). The seismogenic potential and the persisting activity of the accretionary wedge are also clearly documented by the damaging earthquake with an M_w of 6.4 (November 26, 2019; https://bbnet.gein.noa.gr/) whose causative fault is still under investigation and should correspond to the Durresi thrust (ALCS315; Figure 2).

In central Albania, in a more internal sector of the orogen, the crustal rheological conditions allow the formation of a sufficiently thick, weak, and ductile layer between the upper and lower plates (Maggini and Caputo, 2020). Here, the measured GNSS velocities are more consistent with the overall kinematics of this sector of the southern Balkan peninsula, reflecting a nearly south-southwest direction of motion (GNSS stations MUSH, LIXH, DAMG, and M124 in Figure 2).

Close to the Albanian–FYROM border, there is further change in orientation of the relative movements as far as all stations point south to south-southeastward (Figure 2). At a first approximation, this different pattern somehow reflects the crustal stretching pervading the broader Aegean Region (Figure 4). However, if we analyze in more detail the cluster of GNSS stations located in southwest FYROM and east Albania, it is possible to observe some divergence among the western (606_BIS and RE02) and eastern (POGR and KAPS) stations (Figure 2). Such behavior suggests, therefore, that some amount of east-west extension is presently occurring, associated with the activity of the normal faults belonging to the Ohrid Graben (Figure 2).

Regarding the other GNSS stations in continental northeastern Greece and Bulgaria, the southward direction of the velocity vector is more stable. In this sector of the investigated area, it is worth emphasizing that the modules markedly increase southward, being more than double in Chalkidiki (MA10 or MA08) compared to the stations PECH and MA03 or MOMC, close to the border between Greece and Bulgaria (Figure 3). These differences in the velocity confirm the ongoing north-south extension that occurs in the region, which is likely accommodated by the several roughly east-west–trending normal faults described in a previous section.

Finally, GNSS stations in Thassos and Samothraki islands and Turkish Thrace (THAS, SMTK, and RE01) bear a significant westward component of motion likely caused by the drag effect likely induced by the lithospheric dextral shear zone separating Europe from the Anatolian plate (Figure 3).

These observations suggested verifying the suitability and efficiency of an interpolation, driven by the knowledge of tectonic domains, and the local consistency of the velocity solutions in comparison to a standard continuum approach. Hence, we used an interpolation based on the Delaunay triangulation because of its fastness and low memory requirement for computation (e.g., Khatri Chhetri and Zafar, 2013) and the easiness of implementation in other contexts; additionally, it is robust, well tested, and of widespread use. A future improvement of the procedure will employ more rigorous methods, such as the weighted least-squares method, or kriging. For both approaches, we interpolated the data on 49 points, spaced 20 km (12.4 mi). This value is the maximum possible detail that could be reasonably obtained within the investigated area because of the data spacing between the available GPS stations, mostly in the range of 20–40 km (12.4–24.9 mi). For the blocky interpolation, we subdivided the GPS velocity data into clusters related to different tectonic blocks. By comparing Figures 5 and 7, the block subdivision method, compared to the continuum approach, highlights the strain velocity lateral variations bound to different tectonic regimes. We can also observe a better fit between the interpolated velocities with the original observations, like in Thrace, the neighboring of the Anatolia plate boundary, and the Albanian region (Figures 5B, 7B).

The location of the minimum and maximum displacement expected along the corridor slightly varies with the two approaches. In the first case, the minimum is recorded in Albania, Berat district, near the city of Poliçan (Figure 5), with values of 0.6 mm/yr (0.02 in./yr) for the velocity, implying a total displacement in 50 yr of only 30 mm (1.2 in.) (Table 2). With the second approach, the minimum velocity is obtained in Albania, but slightly to the east, close to the Ohrid Graben structure (Figure 7). For both approaches, the maximum expected ground displacement should be of 355 mm (14 in.) in 50 yr, between South Almopia (GRCS068) and Anthemountas (GRCS250) faults (Figures 5, 7), near Chalkidona (Greece).

However, for a pipeline integrity study, it is more crucial to recognize the sections where bending strains may occur. Research on the amount of permanent ground displacement causing pipeline breakage and the defensive measures is continuously improving based on new data and modeling. Based on the results reported by Porter et al., (1991), O’Rourke and Liu (1999), and Pour-Ghaz et al. (2018), a value of permanent ground displacement of 50 mm (1.97 in.) is the minimum value for the onset of damage and cracking. Recent tests on pipelines of 760 mm (30 in.) indicate a total displacement of 76–102 mm (3–4 in.) as the cause of the pipe wall tearing (O’Rourke and Liu, 2012). Hence, since the SGC–TAP will have a diameter of 122 cm (48 in.) for the onshore and 91 cm (36 in.) for the offshore section, we take as a threshold for the permanent ground displacement leading to pipe wall damage 100 mm (3.9 in.).

To identify possible bending strain points, we also calculated the difference between adjacent points along the SGC–TAP, from west to east (i.e., from Apulia to Turkey). Tables 4 and 5 report the results for the continuum and blocky approach, respectively. In the former case, the differential displacements are, necessarily, minimal as far as the interpolation over the whole route minimizes the relative differences. Nevertheless, the conservative threshold of 50 mm (1.97 in.) is exceeded between the points 17 and 18 (across the Ohrid Graben, Albania), 30 and 31 (between the Mygdonia [GRCS100] and Serres [GRCS145] composite sources), 36 and 37, and 37 and 38 (south of the Drama fault; GRCS140).

The differential displacement calculated with the blocky approach along the pipeline expected after 50 yr systematically exceeds the minimum value of 50 mm (1.9 in.) at all sections between blocks, but also within block 5 (points 30–31, as for the continuum approach). In particular, at four block boundaries (block 1–2, 2–3, 3–4, and 6–7 boundaries), the differential displacement exceeds the threshold 100 mm (3.9 in.) in 50 yr. Between the points 8–9 and 15–16 (block 2–3 and 3–4 boundaries, respectively), the long-term cumulative value is 181 (7.1 in.) and 184 mm (7.2 in.) (Table 5). From a tectonic point of view, these two maxima occur at both boundaries of block 3 that separate two distinct geodynamic realms: on the west, the accretionary wedge still characterized by a compressional regime and affected by several reverse faults; on the east, a region affected by crustal extension possibly induced by the local postorogenic collapse diachronically migrating westward (Caputo and Pavlides, 1993). Moreover, the calculated values of horizontal contractional (block 1–2 and 2–3 boundaries) and extensional (block 3–4 boundary) rates are in agreement with those estimated for the primary faults affecting the area (Hoffmann et al., 2010; Guzman et al., 2011).

As for the section between points 4–5 (blocks 1–2 boundary), the estimated cumulative displacement of 105 mm (4.1 in.) is related to the most external fault (basal detachment), separating the Albanian accretionary wedge from the Apulian foredeep (ALCS201; Figure 2). On the other side of the investigated area, between the points 45 and 46 (block 6–7 boundary), the 122 mm (4.8 in.) displacement is associated to the change in tectonic regime affecting the region between Thrace, which belongs to the broader extensional Aegean Region, and the lithospheric shear zone depicted by the North Anatolian fault zone–North Aegean trough. The orientation of the differential vector (Table 5) is indeed parallel to the latter zone, and the calculated modulus represents only a fraction of its present-day slip rate (among others, R. E. Reilinger et al., 1997; Kahle et al., 1998; McClusky et al., 2000; Kozaci et al., 2007; Gasperini et al., 2011). The uncertainty in the block boundaries can be estimated to be in the order of 10–20 km (6.2–12.4 mi), comparable to the distance between two subsequent measurement points. With the increasing number of observations, the spacing between the interpolated points can be reduced, allowing a more precise location of the maximum expected differential displacement.

CONCLUDING REMARKS

Starting from systematically revised and homogenized GPS measurements (Pérouse et al., 2012 and reference therein), we calculated the velocities (modulus and azimuth) as well as the relative uncertainties at numerous closely spaced points along the SGC–TAP route, from Apulia to the border between Greece and Turkey, following two different approaches. The first one considers the whole path as a continuum medium; the second approach assumes the presence of crustal blocks characterized by relatively homogeneous velocities, whose boundaries are caused and defined by the activity of the primary seismogenic sources and where the tectonic regime laterally changes across the area (GreDaSS; Caputo et al., 2012; Caputo and Pavlides, 2013).

The kinematic analyses of the points along the pipeline provide valuable information relative to the regions where differential movements and localized deformation could occur in the future, therefore damaging the infrastructure.

In the reference frame here adopted (Pérouse et al., 2012), our results show that the smallest velocity vectors are in central Albania, near Berat, close to the boundary between the blocks 2 and 3, whereas the highest values are recorded in the central part of the SGC–TAP route, in Greece within the block 5, close to Thessaloniki. These velocities imply that the maximum cumulative displacement forecasted in the 50-yr pipeline lifespan for the ending part of the SGC–TAP is 355 mm (14 in.) south of the South Almopia fault (GRCS068), calculated according to both approaches, with an uncertainty of 55 mm (0.4 in.).

However, for the structural integrity of the pipeline, it is of greater and utmost importance to recognize the sectors of this strategic infrastructure that could be affected in the long term by differential displacements, like shear deformation and bending. Indeed, these phenomena could induce tension or compression to cumulate between the pipeline’s adjacent segments, and this could lead to some damage to the infrastructure. For this purpose, we again considered a 50-yr time window.

The blocky approach to the data interpolation here proposed allows emphasizing the differential motion between neighboring blocks and, hence, between the various segments of the pipeline. Although the expected long-term differential movements are in the order of some centimeters (up to ∼18 cm [∼7.1 in.]), local stresses could rise if deformation will be concentrated along the pipeline. In such a case, even displacement as small as 50 mm (1.9 in.) can activate the breakage. To know whether the expected permanent ground displacement can overpass the indicative threshold of 100 mm (3.9 in.) in 50 yr can allow for taking the necessary counteractions at the right time. We recognized four sectors where the long-term relative movement between neighboring points could be as high as 181 ± 78 mm (7.1 ± 3 in.) and 184 ± 59 mm (7.2 ± 2.3 in.), near Berat (Albania) and in Greece, close to the border with Turkey, respectively (Figure 7; Table 5).