Optimal Routing of Gas Pipelines in Seismic Regions Using an Efficient Decision-Support Tool: A Case Study in Northern Greece

Featured Application

The proposed GIS-based decision-support tool can be reliably applied for the selection of the optimal routing of high-pressure gas pipelines, considering various design criteria and potential earthquake-related geohazards.

Abstract

High-pressure gas pipelines are significantly vulnerable to earthquake-related geohazards (tectonic faulting, slope instabilities, and/or soil liquefaction phenomena). Avoiding geohazardous areas is not always techno-economically feasible, as it would increase the length and cost of the infrastructure. Conversely, crossing these areas may adversely affect the structural performance of the pipeline, leading to unfeasible mitigation measures. Thus, selecting cost-effective, safe, and resilient routing is crucial. This study presents a GIS-based decision-support tool for optimal routing, taking into account, among other criteria, earthquake-related geohazards. The proposed tool considers not only the aforementioned but also more complex earthquake-related geohazards, such as secondary fault ruptures that are non-parallel or even perpendicular to the main fault, which might have been overlooked during the design of existing pipelines. To validate its effectiveness, the present tool is applied in a real case study in northern Greece, where the aforementioned earthquake-related geohazards coexist. Through a GIS-based multi-criteria decision method, various scenarios are examined by assigning different weights to the adopted criteria, and several cost-minimized routes are derived. This tool could be highly beneficial for the pipeline industry since it can assist operators and stakeholders in selecting the optimal pipeline route in geohazardous areas.

1. Introduction

High-pressure pipelines constitute expensive critical infrastructure (CI) that efficiently transport large volumes of natural gas onshore and/or offshore. The anticipated peak in global emissions from natural gas in the near future, combined with the volatile geopolitical situation in northeast Europe and the increasingly stringent environmental policies in response to environmental challenges, has created an urgent demand for the optimum design of environmentally friendly, resilient, and cost-effective CI. These aspects are in line with the sustainable design goals (SDGs) of the European Green Deal [1]. Route optimization is an essential and crucial process in the design phase of a CI since it can minimize the risk of failure and reduce the associated life-cycle costs, ensuring in parallel that the CI will operate safely and reliably throughout its lifetime. Nevertheless, route optimization is often affected by a variety of factors, including geopolitics, financial constraints, environmental considerations, and technical challenges.

Optimal route selection of CI is a complex and demanding task that is typically carried out “manually”, i.e., it heavily relies on the expertise and subjective opinion of engineers. Initially, a suitable routing corridor is identified based on the available and reliable data of the area under examination using existing analog maps. A preliminary survey is then conducted along this corridor, and several alternative CI routes are examined. After analyzing the survey data and ensuring that any potential hazards or obstacles have been avoided, the optimal route is chosen. Nevertheless, this manual approach can be time-consuming, error-prone, and limited in scope. It does not always guarantee the optimal path, as it does not take into account the various project-specific criteria, constraints, and design guidelines that usually complicate the selection of optimal routing.

To address this issue, researchers have developed mathematical models and algorithms to select the optimal route for CI [2,3]. Although such models are not affected by human bias, they encounter difficulties in handling the vast amount of geodata required for such complex spatial multi-criteria applications, as well as with the efficient visualization of the corresponding results. As a solution, the use of Geographic Information System (GIS) has replaced the aforementioned models. GIS offers a robust method for capturing, managing, evaluating, and visualizing information while enabling the efficient handling of geospatial data for multi-criteria analysis. GIS-based tools allow users to overlay different layers of geospatial data and simulate various scenarios, thereby ensuring that decisions are made with high precision and accuracy.

In applications such as the optimal route selection of CI, where, according to Allen [4], approximately 70% of a project’s financial administration and resource allocation are based on geospatial information, GIS offers a valuable solution that overcomes the aforementioned limitations of traditional methods and mathematical models. Furthermore, GIS provides the Least-Cost Path Analysis (LCPA) method, which utilizes a series of sophisticated geoprocessing tools to identify the most cost-effective route on a given cost surface between specified starting and ending points [5].

Relevant research has indicated that utilizing GIS for the optimal route selection of CI can lead to cost savings of up to 30% compared to the traditional process of CI routing [6,7]. Furthermore, the LCPA can be integrated with advanced Multi-Criteria Decision Methods (MCDMs) to improve the process of decision-making. The combination of GIS and MCDMs has been successfully implemented in various fields, including land-use analysis, urban planning and site suitability analysis. In the process of optimum CI routing, GIS and MCDMs have been effectively combined by many researchers (e.g., [5,7,8,9]).

Nevertheless, large-scale CI often cross seismic-prone areas for hundreds to thousands of kilometers, thus being considerably vulnerable to various earthquake-related geohazards, such as strong ground motion, slope instabilities, soil liquefaction phenomena, and seismic fault ruptures at the ground surface. Numerous significant failures of CI resulting from earthquake-related geohazards have been reported over the past few decades. For instance, remarkable structural damage has been identified after the 1971 San Fernando [10], 1999 Kocaeli, Turkey, and Chi-Chi, Taiwan [11], as well as the 2008 Wenchuan, 2009 Italy, and 2010 Chile earthquakes [12,13].

The recent 2023 M_w 7.8 earthquake in Turkey caused significant damage to the Kahramanmaras–Gaziantep natural gas transmission pipeline, leading to the occurrence of heavy explosions and extensive disruption of gas supplies in the cities of Gaziantep and Hatay [14]. In addition, it is worth noting that the impact of earthquake-related geohazards on CI extends beyond direct damage, with substantial indirect losses affecting industries worldwide. Indicatively, a significant part of the Trans-Ecuadorian Pipeline, approximately equal to 70 km, was seriously damaged by the 1987 Ecuador earthquake. As a consequence, 60% of Ecuador’s export revenue was lost, and combined with the 5 months of reconstruction, the economic losses reached 800$ million [15].

Considering the fact that the complete avoidance of seismically geohazardous areas is not always feasible due to various techno-economic constraints (e.g., uneven topography) and geopolitical conflicts, it is evident that earthquake-related geohazards constitute critical factors during the process of optimal route selection. Nevertheless, to the best of the authors’ knowledge, very few relevant publications can be found in the literature that fully take into account the presence of earthquake-related geohazards during the route selection of CI. In most cases, earthquake-related geohazards are either based on simplifications (i.e., named as soil and geology factors, geological and natural hazards, etc.) and rules of thumb (i.e., “vertically crossed obstacles”) without performing a more elaborate analysis, treated as absolute barriers (i.e., that cannot be crossed), or completely ignored [8,9,16,17,18].

Consequently, further research is required for the optimal route selection of CI in areas characterized by the presence of earthquake-related geohazards. Recently, the authors facilitated the process of optimal route selection of large-scale lifelines, such as offshore natural gas pipelines, fiber-optic and power transmission cables, subjected to the earthquake-related geohazards of active seismic faulting and slope instability. The advanced capabilities and functionality of GIS for geospatial processing, as well as the provided LCPA technique for optimal routing, have been successfully combined with MCDMs through the Python computer programming language. By allocating varying weight values to the design criteria, different scenarios were generated to accurately evaluate their effects on the CI and effectively analyze all possible CI routing options in seismic-prone areas [19,20,21].

This article is a revised and expanded version of a paper entitled: “GIS-based multi-criteria decision methodology for the optimal routing of large-scale lifelines against earthquake-related geohazards”, which was presented at the PTC-2024 Conference held in Berlin, 8–11 April 2024 [22]. The current paper presents an upgraded version of the aforementioned computational tool that aims to assist engineers in achieving the optimal route of onshore high-pressure gas pipelines, taking into account a wide range of design criteria. Emphasis is given to the presence of the aforementioned earthquake-related geohazards, including the geohazard of soil liquefaction, which has not been incorporated in previous versions. The main novelty of the present work regarding earthquake-related geohazards is that the proposed decision-support tool has also examined a more complex geohazard, i.e., the combined rupture of the main and secondary seismic faults being non-parallel or even perpendicular to the main fault. In general, this complex phenomenon has not been addressed in the process of CI optimal routing, and as it will be explained in the sequence, it can be crucial for pipelines.

Another valuable contribution of the current work is that the coexistence of all potential geohazards (standard and complex) has been successfully taken into account. Specifically, the effectiveness and reliability of the upgraded version of the computational tool have been validated through its application in a real case study in northern Greece. In this particular area, the aforementioned earthquake-related geohazards coexist within a few square kilometers, thus making the process of optimal routing demanding and challenging. The optimal route proposed by the tool has been compared with the route of an existing high-pressure gas pipeline, where the coexistence and presence of more complex earthquake-related geohazards (i.e., regarding secondary faults potential impact) might not have been adequately evaluated during its routing design phase.

Various scenarios corresponding to designer preferences have been considered by assigning different weight values to the adopted criteria. The proposed tool offers significant flexibility and efficiency to designers involved with energy lifelines since it can assist operators in the decision-making process when selecting the optimal route of CI in such critical areas. At this point, it is important to clarify that the term “optimal CI routing” usually refers to a length-minimized route, which in turn is directly related to a cost-minimized route. Therefore, the optimal CI routing will not only be the shortest in length but also the most cost-effective. Nonetheless, if the proposed routing has to cross geohazardous areas, it is crucial to assess the predicted strain levels on the pipeline and compare them against contemporary international standards and guidelines. In addition, a holistic life-cycle cost analysis should include other aspects, such as construction and maintenance issues in remote areas with steeper slopes, but these are beyond the scope of the current study.

Following this introductory section, Section 2 provides a detailed description of the examined problem and the proposed decision-support tool. Section 3 presents the application of the proposed tool, and in Section 4, the obtained results are discussed. Lastly, Section 5 concludes with the key findings of the study, as well as the advantages and limitations of the decision-support tool.

2. Materials and Methods

2.1. Critical Infrastructure and Earthquake-Related Geohazards

As it has been mentioned, CI is considerably vulnerable to earthquake events and earthquake-related geohazards due to the inertial distress caused by strong ground shaking (Transient Ground Displacements—TGDs) and/or kinematic distress caused by Permanent Ground Displacements (PGDs). Nevertheless, it should be stressed that none of the earthquake excitation components along the three directions are expected to significantly threaten an onshore buried gas pipeline, regardless of the earthquake type and resulting seismic waves. High-pressure natural gas pipelines are characterized by low mass; thus, their dynamic response is significantly different from that of typical above-ground structures, where high vertical forces (i.e., gravitational forces) in conjunction with high seismic accelerations lead to significant inertial forces and structural failures.

Hence, onshore buried gas pipelines are mainly distressed kinematically due to PGDs and External Loading (EL), especially when the fault rupture reaches the ground surface (i.e., in the case of shallow earthquakes). There is no doubt that the structural impact of TGDs can be detrimental for specific types of pipelines (e.g., water and oil) under certain local site conditions. Nonetheless, the distress of high-pressure gas pipelines due to PGDs and EL is more challenging and has a high probability of occurrence during their lifetime. Consequently, they play a key role in the design and optimal route selection [23].

The impact of earthquake-related geohazards, such as slope instabilities, soil liquefaction phenomena, and seismic fault ruptures at the ground surface, on CI is of paramount importance. This is evidenced by the fact that several international guidelines and standards have been developed for the seismic design of pipelines by combining existing analytical and numerical methodologies [24] (e.g., ASCE guidelines [25], American Lifelines Alliance [26], and Eurocode 8 [27]). The main goals of these standards are to (a) evaluate the performance of a pipeline subjected to earthquake-related geohazards, (b) establish pipe strain limits for pipe-soil interaction corresponding to different loading scenarios, and (c) propose several mitigation measures.

Slope instabilities are natural geological events that are often triggered by seismic or other dynamic loads (if any), fluctuations in the level of the groundwater table, soil erosion phenomena, etc., and can cause sudden or gradual mass movements on the ground surface [28]. Earthquake-related landslides are widely recognized as one of the most unpredictable and destructive geohazards that may directly or indirectly affect the structural integrity of CI located in or near mountainous or hilly areas.

Considering that it is not always technically and economically feasible for a CI to completely avoid potentially unstable slopes characterized by large inclinations, earthquake-related landslides may result in excessive EL (e.g., structural distress of the CI due to debris sliding), PGDs (e.g., exposure and/or spanning of the CI), or a combination of these actions across the alignment of the CI [29,30,31]. As a consequence, landslide-related incidents may lead to ruptures and leaks that may have severe environmental impacts, as well as long periods of service disruption [32,33]. In the mountainous regions of the Andes, nearly 50% of CI failures have been attributed to landslides [34], while the European Gas Pipeline Incident Data Group reported that 85% of gas pipeline accidents due to geological reasons in Europe from 2004 to 2013 occurred due to landslide-induced tension or buckling [35].

In parallel, soil liquefaction phenomena, which can occur as a result of repeated cyclic motions in saturated loose cohesionless soils, also play a key role in the structural response of CI. During a strong earthquake, where the shear strength of the soil degrades significantly, liquefaction-induced PGDs, such as ground oscillation, loss of bearing capacity, subsidence, lateral spreading, and buoyancy, may affect the structural performance of the crossing or adjacent CI. Buoyancy (which leads to the uplift of a buried CI approximately equal to its burial depth) and lateral spreading (which mobilizes the full passive pressure against the CI, thus resulting in a lateral movement of the CI) are the most critical effects on CI. Extensive distractions of CI caused by soil liquefaction phenomena have been documented in various earthquakes, including the 1964 Niigata Earthquake, 1983 Nihonkai-Chubu Earthquake, 1993 Kushiro-Oki Earthquake, 1994 Hokkaido-Toho-Oki Earthquake, and 2004 Niigata-ken Chuetsu Earthquake [36].

Finally, the problem of pipe-fault intersection has been adequately investigated with analytical [37,38,39,40], numerical [23,41,42], and experimental [43,44] methodologies. Nonetheless, the research effort, as well as the provisions of the aforementioned international standards, have focused on the intersection of the pipeline with single faults, treating them rather unrealistically as individual entities. Consequently, the presence of secondary faulting in the same region (which might also lead to unexpected failures of an adjacent or crossing CI) is usually ignored. Secondary faults, defined as subsidiary ruptures to a lesser extent in the proximity of the main fault, are mainly attributed to stress redistributions and/or heterogeneities in rock properties and can produce displacements comparable to those of the main fault [45,46]. Secondary faults can have multiple orientations (i.e., parallel or intersecting the main fault), while they can be activated either concurrently (i.e., during the same tectonic episode) or during different tectonic episodes [47,48].

The present study focused on a specific case where secondary faults intersect the tips of the main fault unilaterally, and their orientation is perpendicular to the main fault (Figure 1). It should be noted that both faults are assumed to rupture simultaneously, and that the displacements produced by both faults are of the same magnitude. The presence of such secondary ruptures has been recorded in several earthquakes, such as the 1992 M_w 7.3 Landers (USA), the 2010 M_w 7.2 El Mayor (Mexico), the 2016 M_w 6.5 Norcia (Italy), and the 2019 M_w 6.4 and 7.1 California (USA), causing damage to structures and infrastructure [49,50,51].

Figure 1a depicts a pipeline that perpendicularly crosses the main fault, parallel to and far away from the secondary fault. This constitutes a common case that has also been included in all seismic norms and guidelines used for the design of pipelines. On the other hand, Figure 1b shows a pipeline that is aligned parallel to the main fault but crosses a secondary fault. This constitutes a realistic scenario in engineering practice where a pipeline is aligned parallel to an active tectonic fault in order to prevent crossing it without properly assessing the presence and potential impact of secondary faults. Although this pipeline routing is hypothetically safer, recent studies by the authors’ group revealed that in such cases, the pipeline may even have the same risk as if it intersected the main fault [52,53]. Lastly, Figure 1c illustrates the potential scenario in which a pipeline crosses both the main and the secondary fault. It should be noted that in such cases, the design standards (e.g., Eurocode 8 [27]) suggest that a suitable intersection angle of the pipeline with the (main) fault should be chosen in order to minimize pipe distress. Nonetheless, the norms and engineering practice usually ignore the potential presence of secondary faulting and its consequent impact on the pipeline.

2.2. Description of the Proposed Decision-Support Tool

The proposed decision-support tool consists of three basic steps. Step 1: All required data and geospatial data (i.e., geodata) related to the topography, geology, seismicity, geotechnical conditions, and restricted parts (i.e., “no-go” zones) of the area under consideration are collected and organized into separate thematic layers within a GIS geodatabase. An overlay analysis is performed to generate a composite map that visually represents the spatial relationships between the collected geodata and provides stakeholders and experts with a rapid and detailed digital overview of the area. Having performed the overlay analysis, and after meticulous data processing, each thematic layer is then converted into raster format and transformed into a fundamental design criterion for the Multi-criteria Decision Analysis (MCDA) and subsequently for the LCPA technique.

Step 2: Once the basic design criteria are defined, the MCDA is carried out through the implementation of a proper MCDM that converts the designer preferences into weight factors. Herein, the Analytical Hierarchy Process (AHP) has been employed, in which an n × n comparison matrix (judgment matrix) is constructed, where n denotes the number of criteria. This matrix includes the ranking score that has been qualitatively assigned to each criterion by the designer based on an evaluation scale ranging from 1/9 to 9. Having defined the judgment matrix, the qualitative factors (i.e., ranking score) are then normalized and averaged in order to obtain an average weight for each criterion through pairwise comparisons. Additional numerical measures, such as the consistency ratio (CR), have been established with the aim of evaluating the consistency of the hierarchy. It should be noted that different scenarios can be created at this step based on the assigned weight factors.

Step 3: Having derived the weighting factors, the LCPA technique is employed within the GIS environment to generate multiple length- and consequently cost-minimized (i.e., optimal) CI routes corresponding to each scenario.

The flowchart shown in Figure 2 graphically illustrates the main steps of the proposed decision-support tool, with a comprehensive explanation of each step provided in the following sub-sections. It is worth noting that the proposed decision-support tool has been developed using the powerful capabilities of ArcMap 10.8, a basic application of ArcGIS Desktop [54]. The GIS-based LCPA technique has been automatically integrated with the AHP methodology through suitable scripts written in Python programming language.

2.2.1. GIS Geodatabase and Basic Criteria Definition

The integrated collection of various types of data and geodata that is expertly organized in separate thematic layers and managed into a GIS geodatabase as a single cohesive unit constitutes an essential step in such complex spatial multi-criteria applications. This step enables seamless data retrieval and facilitates efficient decision-making processes. Shapefiles and GeoTIFF have been a popular source of data for GIS, while in cases where digital data is unavailable, existing analog maps are scanned and digitized.

After collecting and storing the necessary geodata, an overlay analysis is conducted by combining the characteristics of the different geodata layers. Overlay analysis constitutes a spatial operation in which two or more different layers of geodata, which are registered to a common coordinate system, are superimposed to virtually represent the spatial relationships between features that occupy the same geographic space. This process results in a composite map, which is commonly used to find locations that are suitable for a particular use or are susceptible to some risk, and consequently, areas, patterns, and trends are effectively identified. It is noted that overlay analysis can be accomplished with any data type (raster or vector). In general, there are two methods for performing overlay analysis: feature overlay (overlaying points, lines, or polygons) and raster overlay.

In the context of optimum routing, a feature overlay has been initially carried out to provide stakeholders and experts with a rapid and detailed digital overview of the area under consideration and empower them to make well-informed decisions. However, after converting all vector geodata into raster format, raster overlay analysis is utilized to derive the desired raster cost (i.e., weighted) surface, which is essential for the LCPA (i.e., routing analysis) to handle areas with multiple overlapping geohazards and other hazards.

To ensure that the computational tool produces accurate and reliable results that can be used confidently in decision-making processes, the utilized geodata has been reclassified and resampled. Reclassification constitutes an attribute generalization technique provided by GIS software that can be applied to both vector and raster data. Specifically, during the reclassification process, the classes of a geodata layer change without altering its geometry, thereby allowing consistent comparisons and analysis across different datasets. In spatial multi-criteria applications, such as the one presented in the sequence, reclassification is crucial for the “criteria standardization” process, where each geodata layer is converted to a thematic layer. On the other hand, geodata resampling involves changing the spatial resolution of a raster dataset. Reprojecting an image to a different coordinate system is another important process that creates an image pixel grid on an alignment other than that of the original image. The value of each pixel in the new image must be computed by sampling or interpolating over a neighborhood of pixels from the corresponding position in the original image. Resampling and reprojecting geodata with different resolutions and mapping projections to a common resolution and projection constitute essential steps before combining and analyzing raster datasets. Herein, all the utilized geodata has been resampled and reprojected to have a uniform spatial resolution and coordinate system, thus eliminating any potential biases or errors [55].

Selecting the optimal route of a CI in real case studies constitutes a demanding decision-making procedure that typically takes into account various design criteria based on the project requirements and constraints, such as budget, accessibility, and environmental issues. Indicatively, “no-go” areas, such as military zones, archaeological sites, and urban and NATURA 2000 areas, are fundamental design criteria since CI crossing is strictly prohibited. Similarly, criteria related to areas characterized by steep slopes, high seismic activity, etc., are also critical and have played a key role in the current study.

2.2.2. Multi-Criteria Decision Method and Weight Factors

A complex decision-making system involves multiple sub-factors (i.e., criteria) that must be meticulously assessed based on their significance within the system. However, assigning weighting values to these factors is a crucial topic, as it often relies on the expertise and knowledge of the experts involved. In response to this challenge, unbiased MCDMs have been developed to effectively address various interrelated objectives. Practically, MCDMs “translate” the qualitative preferences of designers into quantitative weight factors by utilizing specific scientific procedures that systematically decompose complex decision-making systems into sub-factors.

AHP, developed by Saaty in the 1970s [56], provides a consistent framework for evaluating alternatives and selecting the best course of action by considering both qualitative and quantitative factors. Its flexibility and ability to handle complex decision-making systems have made AHP a well-known and highly effective MCDM method that has been extensively applied in numerous complex and multidisciplinary case studies [57]. Through this process, decision-makers are able to prioritize the criteria based on their importance, leading to a more systematic and structured approach to decision-making.

After identifying the criteria for the problem, the following n × n comparison matrix (i.e., judgment matrix), A, is constructed:

(1)

where C_i ($i=1,2,3\dots ,n$) denotes the considered criteria along the matrix rows and columns, while n represents the number of criteria. Each cell of the judgment matrix includes the ranking score, α_ij, which has been qualitatively assigned to each criterion based on an evaluation scale ranging from 1/9 to 9. It should be noted that the values of α_ij and its transpose α_ji are inversely proportional, i.e., $a_{ji}=1/a_{ij}$. A ranking score equal to 9 indicates that the criterion in the row of the matrix is much more important than the criterion in the column, whereas the reciprocal value (i.e., 1/9) implies the opposite. A value equal to 1 indicates that the two criteria are of the same importance and appear along the diagonal of the matrix.

Having defined the judgment matrix, A, the qualitative factors (i.e., ranking score) are then normalized and averaged to obtain an average weight for each criterion through pairwise comparisons. More specifically, the ranking scores in every column of A are summed up, and subsequently, each cell of A is divided by the sum of the specific columns, thus forming the normalized criteria comparison matrix Ā as follows:

(2)

where, ${\Sigma }_k=a_{1j}+a_{2j}+\dots +a_{nj}$ $k,j=1,2,3\dots ,n$ is the sum of each column of A. In the sequence, the priority vector W, which includes the weights for each criterion, is derived after calculating the average of each row according to:

$$\left\{W\right\}=\overline{W_i}={\displaystyle \sum }_{j=1}^n\overline{a_{ij}} i,j=1,2,3\dots ,n.$$

(3)

It should be noted that W determines which of the adopted criteria is the most dominant for the design process. Numerical measures, such as the consistency index (CI) and consistency ratio (CR), have been established to evaluate the consistency of the hierarchy. Initially, vector W_s is defined as the product of A and W, and subsequently, the consistency vector C_v is derived by dividing each cell value of the matrix multiplication result vector by the corresponding priority vector cell as follows:

$$\left\{W_s\right\}=\left[A\right]\left\{W\right\}$$

(4)

$$\left\{C_v\right\}=\left\{W_s\right\}\left\{{\displaystyle \frac{1}{W}}\right\}$$

(5)

The eigenvalue of the algebraic system λ is obtained from the average of the elements of the consistency vector C_v, and it is utilized to calculate the value of CI as follows:

$$CI={\displaystyle \frac{\lambda -n}{n-1}}$$

(6)

Moreover, the CR is calculated as:

$$CR={\displaystyle \frac{CI}{RI}}$$

(7)

where RI (i.e., random index) denotes the CR of a randomly generated matrix, and its value depends on the number of criteria and can be derived from related publications. Finally, when the resulting CR is less than 0.1, the pairwise comparison matrix is consistent; consequently, the AHP methodology is successfully applied without unreasonable pair combinations. More details regarding the theoretical and practical aspects of AHP were provided by Makrakis et al. [21].

The high functionality and advanced capabilities of GIS allow users to employ Python programming language to create customized scripts and tools that can automate and optimize the workflow in GIS applications, which have been utilized herein. Particularly, a Python script was programmed in order to successfully implement the matrices and relationships of AHP. This script has been integrated into the developed GIS model through the Python site package for analysis and data management, ArcPy.

2.2.3. GIS-Based Routing Analysis

Once the routing design criteria have been identified, the MCDA has been performed, and the corresponding weight factors have been derived, a multi-criteria routing analysis is conducted using the LCPA technique in GIS. LCPA constitutes an efficient methodology that can manage a variety of spatial criteria to enhance the effectiveness and reduce the costs of large-scale engineering projects, such as CI. LCPA is based on a variation of Dijkstra’s algorithm [58] and is a highly effective way to identify the shortest and most cost-efficient route on a cost surface, such as a raster surface, between the starting and destination points.

The algorithm first examines all the cells that are adjacent to the starting point in the horizontal, vertical, and diagonal directions, and the one with the lowest cost is identified. This cell is then considered the new starting point, and the next cell with the lowest cost is determined. This process is repeated iteratively until the user-defined origin and destination points are connected. Subsequently, a back-propagation process from the destination to the origin point is conducted using the previously identified cells with the lowest cost, which results in the least-cost path. It is worth noting that one of the most important advantages of Dijkstra’s algorithm and its variations is that they guarantee the solution of the shortest path problem in a computationally efficient manner compared to other route planning algorithms. This is attributed to the fact that they examine all possible cells in the cost surface, identify all suitable paths, and ultimately determine the most cost-effective one between two geographic points.

3. Results from the Application of the Decision-Support Tool

The effectiveness of the proposed computational tool has been confirmed through its application in selecting the optimal route for a high-pressure gas pipeline that crosses a specific area in northern Greece near the cities of Xanthi and Komotini. Figure 3a–c provide a general view of the examined area.

To realistically represent the topography of the examined area, a reliable Digital Elevation Model (DEM) has been utilized. More specifically, a DEM developed from the Copernicus earth observation component of the European Union’s Space Programme has been obtained as a GeoTIFF format file from OpenTopography.org [59]. The resolution of the utilized DEM has been carefully chosen to be equal to 30 m to ensure that the proposed decision-support tool can accurately handle extensive data sets and be successfully applied to large geographic areas while maintaining optimal computational performance. It is noted that global navigation satellite system (GNSS) advanced capabilities (e.g., precise point positioning (PPP)) can be used to increase the resolution of the geodata [60]. Through proper processing within the GIS environment, the map shown in Figure 4 has been generated. The administrative boundaries of Greece and neighboring countries (i.e., Turkey and Bulgaria) have been downloaded as shapefiles from the global database of political and administrative boundaries, namely geoboundaries [61]. The white points on the map denote the approximate locations of the major cities within the examined area.

Figure 5a–e depicts the relative geodata that have been collected and are essential for the application of the proposed decision-support tool. Figure 5a presents the hydrographic network (i.e., rivers and streams) of the examined area, as obtained from the open geospatial database for Greece, geodata.gov [62]. Figure 5b displays the seismic fault zones acquired as linear entities from the global database maintained by the Global Earthquake Model (GEM) Foundation [63]. In addition, by utilizing the adopted DEM and appropriate tools within the GIS software, the map in Figure 5c has been generated, showing the slope inclination of the examined area. Figure 5d illustrates the liquefaction susceptibility map. This map is part of a global liquefaction susceptibility map, which has been created based on the geospatial liquefaction prediction models of Zhu et al. [64]. Finally, Figure 5e presents a map of the epicenters of previous earthquakes over the last 100 years retrieved from the USGS database [65]. The epicenters of the earthquakes near the two secondary fault zones (see Figure 6a) are also highlighted.

Based on the aforementioned geodata, two additional maps have been created that roughly describe the more complex geohazard of secondary seismic faulting and the standard geohazard of soil liquefaction. Specifically, Figure 6a shows the two secondary seismic fault zones adopted in this study, as realistically possible, based on the geomorphology and seismotectonic activity of the examined region. These fault zones are situated on either side of the lake with a perpendicular orientation to the existing Kavala-Xanthi-Komotini seismic fault zone (1st and 2nd secondary faults, to the west and east of the lake, respectively). On the other hand, Figure 6b depicts the potentially liquefiable areas considered herein. The first area encompasses the wider region surrounding the northern part of Vistonida Lake, while the second one includes the coastal areas between the city of Kavala and the city of Alexandroupolis, including the southern part of Vistonida Lake.

Having defined the spatial distribution of both the standard (i.e., seismic faulting, slope instability, and soil liquefaction) and more complex (secondary seismic faulting) earthquake-related geohazards, four scenarios have been identified and examined for the CI routing analysis. It is noted that these scenarios do not include any other design criteria, as in the more holistic and realistic scenario that will be presented in the sequence. Accordingly, they consider three rather extreme cases and one intermediate case:

• Scenario 1: Both primary and secondary seismic fault zones are much more important than slope inclination and potentially liquefiable areas.
• Scenario 2: Potentially liquefiable areas are much more important than the slope inclination and seismic fault zones.
• Scenario 3: Slope inclination is much more important than potentially liquefiable areas and seismic fault zones.
• Scenario 4: An intermediate scenario in which the slope inclination, potentially liquefiable areas, and seismic fault zones are equally important.

Table 1. Weight influence of each of the adopted criteria and sub-criteria for earthquake-related geohazards.

Criteria	Weight Factors (%)				Sub-Criteria	Weight Factors (%)
	Scenarios					Sub-Scenarios
	1	2	3	4		1	2	3
Slope Inclination (°)	9.1	9.1	81.8	33.3	0°–10°	2.8
					10°–20°	24.3
					20°–30°	24.3
					30°–40°	24.3
					40°–50°	24.3
Potentially Liquefiable Areas	9.1	81.8	9.1	33.3	-	-
Seismic Fault Zones	81.8	9.1	9.1	33.4	Main Fault	81.8	9.1	9.1
					1st Secondary Fault	9.1	81.8	9.1
					2nd Secondary Fault	9.1	9.1	81.8

indicates the CI design criteria and sub-criteria that have been adopted, along with their respective weight factors (i.e., influence), as obtained from the implementation of AHP. Low values of influence indicate that a criterion or sub-criterion has a marginal impact on the CI routing analysis; thus, the resulting CI route could potentially cross this geohazardous zone. In contrast, a high influence value implies that a geohazardous zone must not be crossed when performing CI routing analysis. The criterion of slope inclination has been distinguished into five classes (sub-criteria) representing zones of slope inclination ranging from 0° to 50°, and the preference is given to pipeline passage via (almost) flat areas, i.e., all higher inclination zones are assigned a much higher weight factor. Additionally, three sub-scenarios have been examined regarding the geohazard of seismic faulting by assigning the specific weight factors to the main, the 1st secondary, and the 2nd secondary faults, aiming to their avoidance in the 1st, 2nd, and 3rd sub-scenarios, respectively.

After assigning weight factors to the selected design criteria and sub-criteria, the LCPA technique has been implemented to derive several cost- and length-minimized CI routing alternatives. Figure 7a–d presents the obtained CI routes corresponding to scenarios 1, 2, 3, and 4, respectively. The yellow points denote the starting and ending points of the proposed routings. In particular, in all scenarios, three optimal routes have been derived depending on which fault zone is considered more critical and has to be avoided. According to the weight factors in the right column of Table 1 in all subplots of Figure 7, Route I, Route II, and Route III refer to sub-scenario 1 (81.8% for the primary seismic fault zone), 2 (81.8% for the 1st secondary fault zone), and 3 (81.8% for the 2nd secondary fault zone), respectively.

A more holistic and realistic scenario (i.e., scenario 5) has also been defined and examined. In addition to the aforementioned earthquake-related geohazards, geospatial data pertaining to the existing road network, urban areas, and NATURA 2000 areas has been collected from the OpenStreetMap exports on Humanitarian Data eXchange (HDX) [66], Geodata.gov.gr [62], and the datahub of the European Environment Agency [67], respectively, and has been used as additional design criteria. The collected geodata is depicted in Figure 8, and

Table 2. Qualitative prioritization of the adopted criteria for scenario 5 with a consistency ratio equal to 0.08.

	Slope Inclination	Potentially Liq. Areas	Seismic Fault Zones	Road Network	Natura 2000 Areas	Urban Areas
Slope Inclination	1	2	2	3	3	3
Potentially Liq. Areas	1/2	1	2	4	4	3
Seismic Fault Zones	1/2	1/2	1	2	3	3
Road Network	1/3	1/4	1/2	1	1/3	1/4
Natura 2000 Areas	1/3	1/4	1/3	3	1	1/2
Urban Areas	1/3	1/3	1/3	4	2	1

provides further insight into the qualitative prioritization of the selected criteria. It is noted that similar prioritization tables have been used to obtain the weight factors in Table 1 according to the standard process of AHP. In the sequence,

Table 3. Weight influence of each of the adopted criteria and sub-criteria for scenario 5.

Criteria	Weight Factors (%)	Sub-Criteria	Rating Score (%)
			Sub-Scenarios
			1	2	3
Slope Inclination (°)	30.0	0°–10°	3.3
		10°–20°	11.8
		20°–30°	17.9
		30°–40°	28.9
		40°–50°	38.1
Potentially Liquefiable Areas	25.6	-
Seismic Fault Zones	17.8	-	-
Road Network	5.9	-	-
Natura 2000 Areas	8.7	-	-
Urban Areas	12.0	-	-

summarizes the influence of the considered geodata as obtained after the implementation of AHP. The results of the routing analysis for this scenario are shown in Figure 9. The proposed route has been compared with the existing route in this area of the high-pressure Trans Adriatic Pipeline (TAP) (derived from the corresponding pipeline routing report [68]), which transports natural gas from the Greek−Turkish border to southern Italy.

4. Discussion

Figure 3a–c, Figure 4, and Figure 5a–c clearly demonstrate the significant interest in the examined area owing to its distinctive topography and geomorphology. Specifically, the presence of Vistonida Lake in the central part of the examined area (Figure 3a,b) indicates significant seismotectonic activity, possibly resulting from lateral and vertical movements of the Earth’s crust, which has also been documented in previous studies [69]. The seismotectonic activity is further elucidated by examining the map in Figure 5b, which highlights the dense concentration of seismic fault zones. The approximate location of seismic fault zones has been based on the findings of Mountrakis et al. [70].

Herein, particular focus has been given to the Kavala-Xanthi-Komotini seismic fault zone, which runs adjacent to Vistonida Lake and extends for tens of kilometers in a west-east direction within the examined area [71]. In such areas, characterized by the presence of earthquake-related geohazards, a vast amount of geodata is needed to ensure reliable and accurate mapping, typically obtained via geophysical surveys and geological maps. In parallel, the examined area is characterized by a notable variation in elevation, with mountains of considerable height located just a few kilometers from sea level (Figure 3c). This fact is corroborated by the slope inclination map shown in Figure 5c, which reveals that the slope inclination ranges from 0° to 50°, with steep slopes to the north of the seismic fault zone and almost flat terrain to the south.

As shown in Figure 5a, the area is intersected by numerous wet valleys. It is worth noting that most rivers and streams in the examined area flow perpendicular (i.e., north–south direction) to the main seismic fault zone. Taking this into consideration, in conjunction with the considerable length of the Kavala-Xanthi-Komotini seismic fault zone and the potential stress redistributions and/or heterogeneities in rock properties, it is probable that this seismic fault may rupture in smaller parts, leading to the development of secondary seismic fault zones in an almost perpendicular direction. Therefore, two rivers located on either side of Vistonida Lake with a roughly perpendicular orientation to the main seismic fault zone have been crudely defined as secondary faults, designated as the 1st and 2nd secondary faults, respectively (Figure 6a). It is worth noting that, as depicted in Figure 5e, the earthquake epicenters located on both sides of Vistonida Lake provide evidence supporting the presence of secondary faults in these areas. Furthermore, the liquefaction susceptibility map presented in Figure 5d indicates that a significant part of the examined area located south of the main seismic fault zone (i.e., on flat terrain) is prone to earthquake-related geohazards of soil liquefaction. This is mainly attributed to the presence of a large number of rivers and streams, which are related to the transportation and deposition of loose and potentially liquefiable sediments. The current study focused on two zones of potentially liquefiable areas, as shown in Figure 6b, namely, potentially liquefiable areas 1 and 2. The former covers a wider area of the northern part of Vistonida Lake, while the latter includes a wider area of the southern part of Vistonida Lake, as well as a significant part of the regions near the coastline of the examined area.

At this point, it is important to mention that the reliability of the geodata sources, in conjunction with the complex nature of the examined problem (i.e., the large extent of the area under consideration, as well as the uncertainties related to earthquake-related geohazards), make their validation process challenging. Accordingly, there are some limitations regarding the employed geodata. For instance, seismic faults are represented within the GIS environment in a simplified manner as linear entities with starting and ending points and specific dimensions. An accurate mapping of critical factors, such as slip potential, path propagation, and secondary faults, constitutes a complex and demanding task. In parallel, the potential of soil liquefaction phenomena requires a comprehensive investigation of the mechanical and physical properties of the underlying soil layers. Hence, there is no doubt that apart from the vast amount of geodata required for such applications, conducting geophysical surveys and in-situ tests constitutes an essential factor in order to ensure the reliable and accurate assessment of all regional earthquake-related geohazards, which in turn leads to more precise hazard maps.

The first CI routing analysis has been focused solely on “geotechnical” criteria (scenarios 1 to 4), including the standard geohazards of primary seismic faulting, soil liquefaction, and slope instability, as well as the more complex geohazard of secondary seismic faulting. Slope instability constitutes a destructive earthquake-related geohazard. Undoubtedly, the probability of slope failure is directly associated with the mechanical and physical properties of the slope. Slope mechanical properties refer to the material properties, while physical properties involve the geometry of the slope (i.e., height and inclination). It is noted that a high slope inclination does not always indicate an increased risk of slope failure. For instance, a steep rock slope can be considered safe; however, rockfalls can occur in cases where the rock is not intact, which constitutes a major threat to the structural integrity of the crossing or nearby CI.

The current study focused on the geometrical characteristics of the slope, with a particular emphasis on slope inclination. Avoiding steep slopes is essential to prevent crossing via potential landslide zones, which may have an unfavorable impact on the structural integrity of CI [72,73]. Additionally, crossing slopes with high inclinations can significantly increase the cost of the CI, either by increasing its length or by requiring expensive and even impractical measures to mitigate the impact of potentially unstable slopes on the CI.

To address this issue in the current study, slope inclination has been classified into five classes (i.e., sub-criteria), ranging from 0° to 50°, as outlined in Table 1. Quite high weight factors (i.e., 24.3%) have been assigned to slopes with inclinations higher than 10°, while a relatively low weight factor (i.e., 2.8%) has been assigned to low slope inclinations (0°–10°). This variation in weight factors ensures that the proposed CI routes pass through almost flat areas and avoid all higher slope inclinations, similar to the strategy that was followed in the design of the existing TAP route in the examined region.

As far as the geohazard of primary seismic faulting is concerned, it has to be stressed that only the Kavala-Xanthi-Komotini seismic fault zone located near the Vistonida Lake has been considered as a design criterion. This is due to the fact that the other fault zones had no effect on the CI routing analysis. It is noteworthy that in engineering practice, reducing the fault-CI intersection length (i.e., vertical crossing) is considered preferable for minimizing the risk of CI failure [74]. Nonetheless, in reality, the intersection angle may be affected by uncertainties related to fault rupture characteristics. In addition, in such cases, the fault type plays a crucial role in the type and the magnitude of pipeline strains [75]. Regarding the earthquake-related geohazard of soil liquefaction, two potentially liquefiable areas have been identified. The main difference between these areas is that passing through the potentially liquefiable area 2 may be unfeasible for any engineering infrastructure applications due to potential soil erosion phenomena, as well as socio-economic and other constraints associated with the presence of CI in coastal areas.

The coexistence of all the aforementioned earthquake-related geohazards within a relatively small area intensifies the complexity of the routing analysis; thus, the process of optimal route selection is more challenging. Specifically, the primary seismic fault zone acts as a “physical boundary” between the flat terrain and steep slopes, while Vistonida Lake stretches from near the primary seismic fault zone to the coastal region. Thus, two narrow “corridors” are created between the seismic fault zone, the Vistonida Lake, and the seacoast. In addition, the western end of the main seismic fault zone lies in close proximity to the coast, while the secondary seismic fault zones are perpendicular to the primary fault zone and extend close to the coast. It is evident that avoiding all these fault zones, as well as potentially liquefiable areas, would not lead to an optimal alternative, as it would increase the length of the examined CI and the probability of failure (i.e., passing through steep slopes north of the primary seismic fault zone). As a consequence, alternative CI routings will inevitably cross one, two, or even three seismic fault zones and are expected to pass either within or at the boundaries of potentially liquefiable areas.

Results from Figure 7a–d clearly demonstrate that in all the examined “geotechnical” scenarios, the criterion of slope inclination played a significant role in routing analysis. Almost all of the proposed routings completely avoid steep slopes located north of the primary seismic fault zone. Similarly, crossing the potentially liquefiable area 1 has been successfully avoided in all routings. However, several of the proposed routings pass through either a part (small or large) or the entire potentially liquefiable area 2.

Routes I and II in Figure 7a correspond to scenario 1 (81.8% seismic fault zones) and particularly to sub-scenarios 1.1 and 1.2 (i.e., 81.8% primary and 81.8% 1st secondary seismic fault zone, respectively). These routes strategically cross either the primary or the 1st secondary seismic fault zones almost vertically, minimizing the length of the routing exposed to this geohazard. To achieve this, part of Route I passes through gentle slopes to the west of the examined area, running parallel and in close proximity to the primary fault before reaching the vertical CI-fault intersection point. The significant impact of the primary seismic fault zone in sub-scenario 1.1 is revealed by the fact that Route I crosses this seismic fault zone only once, in contrast to Route II. In sub-scenario 1.3 (81.8% 2nd secondary seismic fault zone), Route III, the longest route, successfully avoids the intersection with the 2nd secondary seismic fault zone. It is noted that Routes I and II intersect with the 2nd secondary seismic fault zone at different locations than Route III. Nonetheless, it can be seen that in contrast to Routes I and II, which successfully avoid both potentially liquefiable areas, Route III passes through potentially liquefiable area 2.

The 81.8% weighting factor assigned to the potentially liquefiable areas in scenario 2 (i.e., Figure 7b) has led to routes that pass through the narrow “corridor” between the primary seismic fault zone and the potentially liquefiable area 1. However, the proposed routings have not completely avoided potentially liquefiable area 2, thus crossing a part located to the west of the 1st secondary seismic fault zone. The steep slopes north of the primary fault zone are successfully avoided, while the low influence of the seismic fault zones has resulted in routes that cross both the primary and the two secondary seismic fault zones. Nonetheless, minor differences among the proposed routes are evident at the intersection of the two secondary seismic fault zones.

Figure 7c clearly illustrates the significant impact of slope inclination on the CI routing analysis. All the proposed routing alternatives pass through the flat coastal region to the south of the examined area, avoiding the steep slopes, the potentially liquefiable area 1, and the intersection with the 2nd secondary seismic fault zone. The only variation among the routes is the intersection with the 1st secondary seismic fault zone, with Route II crossing at a different point in this zone. This discrepancy can be attributed to the fact that Route II represents sub-scenario 3.2, where 81.8% of the weighting influence has been assigned to the 1st secondary seismic fault zone. However, it is stressed that the routing alternatives resulting from this scenario cross almost the entire potentially liquefiable area 2. Moreover, they are prohibited by the aforementioned restrictions in coastal areas.

The routings resulting from the intermediate scenario in Figure 7d are very similar since they all pass through the narrow “corridor” between the primary seismic fault zone and the potentially liquefiable area 1. Route II intersects almost vertically the 1st secondary seismic fault zone, while Route III completely avoids the 2nd secondary seismic fault zone. It is important to mention that scenario 1 resulted in the shortest proposed routings (i.e., Route II = 173 km and Route I = 177 km), while Route III from scenario 1, Routes I to III from scenario 3, and Route III from the intermediate scenario are equal to 187 km.

In order to conduct a more holistic and realistic CI routing analysis, additional geodata have been collected, processed, and utilized as design criteria for the process of optimal route selection. It has to be stressed that the rationale behind the qualitative prioritization and the resulting weighting factors presented in Table 2 and Table 3, respectively, was to follow as closely as possible the routing “design principles” of the existing TAP route, which crosses the area of interest, and all the aforementioned restrictions that reduce potential routing options that can be explored by the automated tool.

Accordingly, Table 2 presents the judgment matrix that summarizes the qualitative prioritization of the employed design criteria in terms of ranking scores ranging from 1/9 to 9. The CR of this matrix is lower than 0.1 (i.e., CR = 0.08), indicating the transitive property of the matrix, which is essential for the efficient implementation of AHP. The ranking score assigned to the adopted design routing criteria follows the same “pattern” as in TAP design, e.g., regarding the criterion of slope inclination, the hilly areas are considered unfavorable due to constructability and practical reasons and to avoid landslide-prone zones (i.e., soil slopes with higher inclinations). Figure 8 clearly illustrates that a significant part of the examined area belongs to NATURA 2000 areas. Specifically, the regions surrounding the 1st secondary seismic fault zone and the potentially liquefiable areas are NATURA 2000 areas. This makes optimal CI route selection even more complex.

In addition, the examined area is characterized by numerous urban areas with populations ranging from 2000 to over 10,000. These urban areas are evenly distributed on both sides of Vistonida Lake, with the cities of Kavala, Xanthi, Komotini, and Alexandroupolis having the largest populations. It is worth noting that many of these urban areas are located along the primary seismic fault zone, thus requiring further research beyond the scope of the present study. Finally, an extensive road network that connects the aforementioned urban areas can be seen in Figure 9.

Based on the qualitative prioritization shown in Table 2, the resulting weighting factors that correspond to all the considered design criteria are summarized in Table 3: slope inclination (30.0%), potentially liquefiable areas (25.6%), seismic fault zones (17.8%), as well as the additional geodata; i.e., road network (5.9%), NATURA 2000 areas (8.7%), and urban areas (12.0%). It is highlighted that, in accordance with the fundamental principle of AHP, the weight factor assigned to each criterion is divided among its sub-criteria (if any). As a result, the sub-scenarios ultimately receive a lower weight factor than the total weight factor of the specific scenario. Note that in holistic scenario 5, the design criterion of seismic fault zones has not been further divided into sub-scenarios, as in the geotechnical routing analysis. This decision was made to ensure that the impact of this geohazard is not reduced in the routing analysis.

Figure 9 shows the route derived from the realistic scenario 5. The proposed route has been compared with the TAP route in this area. Initially, both routes have similar lengths, with the TAP route being 5 km longer than the proposed route. It is evident that almost all NATURA 2000 areas are successfully avoided on both routes. Nevertheless, the proposed route almost vertically crosses the NATURA 2000 area near the origin point east of the examined area, while the TAP route extends further into this area. The two routes avoid both the potentially liquefiable areas by passing through the narrow “corridor” between the primary seismic fault zone and the potentially liquefiable area 1. Hence, the potentially problematic coastal region is completely avoided.

It can be noticed that the two CI routings are almost identical in the central part of the examined area, both crossing the 1st and the 2nd secondary seismic fault zones, with minor discrepancies being observed at the intersection point. However, it is noteworthy that the proposed CI routing passes through gentle slopes north of the primary seismic fault zone. This is attributed to the weighting factors assigned to the primary fault zone, the road network, as well as to the urban areas. The TAP route, on the other hand, passes closer to the urban areas in the eastern part of the examined area, crossing the primary seismic fault zone at a different location.

Nevertheless, it should be stressed that the comparison with this existing engineering project is only indicative due to the lack of full data derived from topographic surveys and geotechnical studies, information regarding the actual project’s routing, such as the adopted design criteria and the assigned weight values, and the resulting assumptions and simplifications. Furthermore, the inclusion of fragility curves in the tool, resulting from seismic vulnerability analysis (e.g., [76]), would enable designers to fully capture the dynamic distress and seismic resilience of the examined critical infrastructure.

5. Conclusions

High-pressure gas pipelines are highly vulnerable to earthquake-related geohazards (tectonic faulting, slope instabilities, and soil liquefaction phenomena). Since the complete avoidance of such geohazardous areas is not always feasible from a techno-economical perspective, crossing these areas may detrimentally affect the structural performance of pipelines, necessitating expensive and impractical mitigation measures. Thus, selecting cost-effective, safe, and resilient routing is of paramount importance. Optimal route selection, which is one of the main activities when designing critical infrastructure, is characterized by a multidisciplinary and complex nature due to the vast amount of data, as well as the complexities and uncertainties involved.

The current work presents a GIS-based decision-support tool that combines in a computationally efficient manner: (a) the advanced capabilities and functionality of GIS for geospatial processing, (b) a multi-criteria decision method that translates the qualitative prioritization of the adopted design criteria into weight factors, and (c) the sophisticated GIS-based LCPA technique that leads to the optimal path between two geographic points. The developed computational tool can be efficiently applied for the optimal route selection of critical infrastructure, taking into account, among other criteria, the presence of earthquake-related geohazards. This tool considers not only standard earthquake-related geohazards but also more complex geohazards, such as secondary fault ruptures being non-parallel or even perpendicular to the main fault, which might not have been fully considered during the design of existing pipelines.

The proposed decision-support tool is applied to a particular area in northern Greece. Four different “geotechnical” scenarios have been examined considering the presence of earthquake-related geohazards, while a holistic -and more realistic- scenario has been additionally identified, taking into account various design criteria. Several length- and cost-minimized routings are obtained. The application of the computational tool led to the following useful conclusions:

• The examined area in northern Greece, between the cities of Kavala and Alexandroupolis, is prone to earthquake-related geohazards. Both standard and more complex earthquake-related geohazards coexist within a few square kilometers. Thus, the optimal routing of critical infrastructure in such areas constitutes a very demanding and challenging process.
• The assigned weight value corresponding to the criterion of slope inclination considerably affected the orientation of the proposed routings. This is evidenced by the fact that almost all of the proposed alternatives resulting from the “geotechnical” scenarios pass through the flat area south of the Kavala-Xanthi-Komotini main fault.
• Crossing a potentially liquefiable area located near the coast may be additionally unfeasible in engineering applications due to various restrictions associated with the presence of any infrastructure in coastal areas.
• Assigning high importance to seismic fault zones, which results in pipeline routings that cross the fault zones almost vertically, has led to up to 8% shorter routes compared to the other examined “geotechnical” scenarios.
• The route resulting from the holistic scenario is qualitatively compared to the existing route of the Trans Adriatic Pipeline (TAP). Both routes have similar lengths; they completely avoid potentially liquefiable areas while crossing the two considered secondary seismic fault zones. Nonetheless, the proposed routing does not cross the western part of the primary seismic fault zone as it passes through the gentle slopes in the north of the study area.

Certainly, automated decision-support tools can considerably facilitate the process of optimal route selection by reducing the probability of failure and life-cycle costs of the examined infrastructure and minimizing the impact of human errors and bias. However, they should be used on a case-by-case basis, taking into account all relevant factors in conjunction with the expertise of engineering professionals.

The methodology of the proposed decision-support tool is versatile and highly adaptable since it can be effectively applied across a wide range of terrain types, regardless of geomorphology, topography, and other critical factors. The only variable that varies in regions with different geomorphologies is the considered criteria. For instance, in a completely flat region, the geohazard of slope inclination is not a serious concern, and consequently, the corresponding criterion could be “deactivated”. In such cases, other hazards and criteria, such as flooding phenomena, can be prioritized. Further improvements to the proposed computational tool are required to strengthen its present capabilities, focusing on the enrichment of the criteria adopted during the routing analysis, automatization of the whole process, and adopted optimization methodology.