Scenario-based seismic risk assessment for buried transmission gas pipelines at regional scale.

: Buried gas pipelines in seismic-prone regions may suffer leaks or breaks as a consequence of an earthquake, especially if the pipeline is subjected to large differential displacements due to geotechnical failures (e.g., landslide or liquefaction). This paper presents a methodology to assess the risk of a gas pipeline infrastructure at a regional level in the aftermath of a seismic event. Once earthquake characteristics, such as magnitude and epicenter, are known, seismic intensity measures (IMs), such as peak ground acceleration (PGA) and peak ground velocity (PGV), are estimated at the location of each pipe through a simulation-based procedure. The potential updating from real-time data coming from accelerometric stations is considered. These IMs are then used to study the cascading landslide and liquefaction hazards, providing a hybrid empirical-mechanical-based estimation of permanent ground displacements (PGD). With the aid of damage and fragility functions from the literature, loss figures and damage maps are derived as decision-support tools for network managers and stakeholders. Losses provide a preliminary estimation of repair costs, whereas damage maps support the prioritization of inspections in the aftermath of the event. The risk methodology is a novel combination of cutting-edge and consolidated approaches. Firstly, different cross-correlation models between PGA and PGV are included. Secondly, a new three-phase back-to-back geotechnical approach is provided for both landslide and liquefaction, representing (1) the susceptibility, (2) the triggering, and (3) the PGD estimation phases. The 1976 Friuli earthquake and the high-pressure gas network of northeast Italy are used as a test-bed scenario for the risk methodology aimed at emphasizing pros and cons of the different alternative options investigated. DOI: 10.1061/(ASCE)PS.1949-1204.0000330. This work is made available under the terms of the Creative Commons Attribution 4.0 International license, http://creativecommons.org/licenses/by/4.0/.


Introduction
Earthquakes can severely damage buried gas transmission networks (O'Rourke and Palmer 1996). As a result, reliable damage assessment in the aftermath of an earthquake and the development of enhanced fragility models are among the main challenges to be undertaken in seismic engineering of buried pipelines (Pineda-Porras and Najafi 2010). In recent decades, great efforts have been made to improve the methodologies for pipeline vulnerability assessment. For example, Lanzano et al. (2015) compiled one of the most extensive databases of observed damage to pipelines after earthquakes, providing an excellent base for the creation of empirical fragility models. It has been emphasized that buried gas transmission pipelines are particularly vulnerable to geotechnical failures such as landslides and liquefaction phenomena, because they generate large ground deformations, creating large stress concentrations at the joints (Wham and O'Rourke 2015). Edkins et al. (2016) gave a detailed description of observed damage in the aftermath of an earthquake to pipeline systems of different materials.
Damage to gas pipelines can lead to issues in energy production; if the connection with thermoelectric power stations is interrupted, it can affect everyday life, stopping energy flow to homes. Moreover, the damage can trigger major cascading accidents, such as fires and explosions, also known as natural hazard triggering technological disasters. Thus, there is a pressing need for the development of a reliable risk-based decision support system (DSS) that can generate loss figures and damage maps of infrastructure immediately after an earthquake event. Such information can then be used to prioritize inspections and interventions to rehabilitate the gas transmission service as fast as possible (Elsawah et al. 2016).
Several relevant studies are available in the literature on the seismic risk assessment of buried pipelines. Initially, only scenariobased ground shaking was considered (Hwang et al. 2004;Toprak and Taskin 2007). Subsequently, the importance of geotechnical failures was emphasized and geotechnical models have been incorporated in the risk methodologies (e.g., Mousavi et al. 2014). O'Rourke et al. (2014) were among the first to emphasize the importance of having reliable underpinning data to produce robust and comprehensive risk figures. Finally, Esposito et al. (2015) emphasized the importance of correlation models among intensity measures and envisaged the possibility of adopting simulationbased hazard-assessment procedures.
Building upon the seminal research mentioned previously, this paper presents a Monte Carlo simulation-based methodology that can generate preliminary information about the status of a buried gas pipeline infrastructure after a seismic event. This methodology introduces several novelties with respect to previous works. First, the intensity measures at several locations are investigated considering different forms of geographical correlations and crosscorrelations between intensity measures, i.e., peak ground acceleration (PGA) and peak ground velocity (PGV). Second, a procedure to anchor the ground shaking simulations to observed values is presented. Third, more-refined large-scale geotechnical procedures are developed through a three-phase study of (1) susceptibility, (2) triggering, and (3) ground displacements caused by both landslide and liquefaction. Previous works focused mainly on landslides and adopted the simplified procedure recommended in HAZUS (FEMA 2004). Finally, existing empirical vulnerability and fragility functions are used to compute loss curves and damage maps, respectively. This paper represents for the first time the earthquake consequences on gas lifelines as expected damage maps, which is considered to be a key deliverable for informed decision making.
The following sections describe the methodology and apply it to a virtual case study; a portion of the high-pressure gas transmission system in the northeast region of Italy is analyzed under the effect of the 1976 Italian Friuli event. This case study is virtual because the pipeline network considered herein (based on 2013 information) was not completely built at the time of the earthquake event (i.e., 1976). The combination, however, of the model of a real gas network with an actual seismic scenario is deemed to be a useful platform to demonstrate the efficiency and limitations of the methodology.

Methodology
Risk assessment for spatially-distributed infrastructures, such as road networks and buried pipeline systems, requires (1) the identification of the main properties of critical components of the infrastructure, their failure mechanisms, and the associated vulnerability models; (2) the acquisition of the underlying data complementary to the infrastructure (e.g., the soil geotechnical properties); (3) a hazard assessment procedure capable of providing seismic intensity measures at regional scale; and, finally, (4) a risk quantification procedure consisting of the convolution of hazard and vulnerability models. Figs. 1 and 2 provide a graphical representation of the framework developed in this study involving the steps to assess the network hazard ( Fig. 1) and risk (Fig. 2).
For a buried gas pipeline system, specific information is needed [ Fig. 1(a)]: the geometry of the system; the buried depth; the service pressure; and the diameters, materials, and joint typologies of the pipes. Such characteristics are necessary for choosing a suitable vulnerability model among those available in the literature (e.g., ALA 2001; Piccinelli and Krausmann 2013;Lanzano et al. 2015).
Once the system geometry is available, a discretization of the pipeline is needed; specifically, each element of the system is divided in segments with a length (L) which is small enough to neglect behavior/property variability along the segments. Each pipeline segment is represented by its midpoint. Complementary data to the infrastructure is all the information necessary to define the topographic, geologic, and geotechnical problems of the buried system [ Fig. 1(b)]: the digital elevation model (DEM) map; the slope map; the shear wave velocity map of the upper 30 m of soil (V s;30 ); the geolithology map, representing the local lithology from which, in the absence of detailed surveys, it is possible to obtain the geotechnical properties (friction angle ϕ, cohesion c, and selfweight γ); and the groundwater table map.
For a postevent seismic risk assessment, information about the earthquake also is necessary [ Fig. 1(c)]. Essential data are the event magnitude and the epicenter location, or the fault plane, which are normally publicly available soon after the event. Additional data, allowing the refinement of the hazard model, are the strong motion records from seismic stations. Once earthquake data are available, the hazard model can be developed based on predefined ground motion prediction equations (GMPEs) [ Fig. 1(d)]. The hazard model is then used to generate the seismic intensity measures (IMs): the peak ground acceleration and the peak ground velocity. For buried infrastructure systems, the permanent ground displacement (PGD) also is required; this is the IM typically adopted for the interpretation of geotechnical failures and should not be confused with peak ground displacement. This study presents a new back-toback three-phase geotechnical model allowing the quantification of the PGD due to landslide and liquefaction phenomena.

Seismic Hazard Model
To account for the uncertainty in seismic actions due to a specific earthquake, this study adopts a stochastic simulation-based procedure to estimate seismic IMs (in terms of PGA and PGV) at the midpoint of each pipeline segment ðIMÞ. Specifically, IMs are assumed to follow joint lognormal distribution, with central values ðIMÞ computed by a GMPE and covariance calculated with a correlation model. The following formulation is adopted: where σ I = intraevent standard error of the selected GMPE; C = correlation matrix; and σ 2 I · C = covariance matrix (Σ) of the lognormal distribution. Such an approach is largely adopted to predict IMs probabilistically due to a specific earthquake scenario (e.g., Wald et al. 2006).
This study investigated three main correlation models [ Fig. 1(e)]: (C1) uncorrelated IMs, (C2) IMs spatially correlated according to Goda and Hong (2008), and (C3) IMs correlated and crosscorrelated according to Weatherill et al. (2014). The spatial correlation provides the correlation between the same IM (i.e., PGA or PGV) at different geographical locations. The cross-correlation provides the correlation between values of PGA and PGV. For the first case, C1, the correlation matrix is equal to the identity matrix, therefore IMs at different locations will be completely uncorrelated. For C2, terms of the correlation matrix (ρ i;j ) are obtained as a function of the site-to-site distances (R i;j ) between the midpoints of segment i and segment j Several equations for calculating the correlation coefficients are available in the literature (e.g., Goda and Atkinson 2010;Iervolino 2011, 2012). The functions to compute ρ i;j are different for PGA and PGV and they differ for different geographical region and seismogenetic contexts. Correlation terms are equal to 1 if the interdistance is equal to 0, and they decrease to 0 for very large values of R i;j .
Finally, for C3, the correlation matrix contains both spatial and cross correlation terms where C PGA and C PGV are obtained according to Eq. (2), and the generic term of the matrix C PGA;PGV is equal to where ρ PGA;PGV = correlation term between PGA and PGV, and can be estimated by bespoke statistical analyses based on recorded data in a given geographical region. If information from seismic stations is available, the hazard model can be further refined by updating the probability model with the observed IM values. To demonstrate this possibility, correlation models C2 and C3 were modified according to the procedure described by Miano et al. (2016) in order to consider available information from seismic stations. These two additional models are named C4 and C5, corresponding to C2 and C3, respectively. In both C4 and C5, the central values of the distribution and the covariance matrix are defined according to In Eq. (5), IM 1 = intensity measure calculated with the GMPE for the points of interest (e.g., the midpoint of the pipes in this study); IM 2 = intensity measure calculated for the seismic stations; Σ IM1 and Σ IM2 = covariance matrixes for the points of interest and for the seismic stations, respectively; and Σ IM1;IM2 = crosscovariance matrix between the points of interest and the seismic stations' locations. Let X indicate the intensity measure acquired at the seismic stations; the statistics of the joint distribution of the IMs at the sites of interest, updated with the acquired data X, can be computed according to Eqs. (6) and (7) Therefore IM 1j2 and Σ IM1jIM2 can be used in Eq.
(1) to simulate IM fields that are anchored to real observations.

Geotechnical Model
For each of the five correlation models described previously [ Fig. 1(e)], it is possible to obtain simulated intensity measure fields of PGA and PGV [ Fig. 1(f)], however, such IMs are not enough for the case of buried infrastructural systems. This is because, for such kinds of infrastructure, PGD is also required because several vulnerability/fragility models are expressed in these terms. PGD cannot be calculated by GMPEs but it is inferred from PGA and PGV because none of the GMPEs predicts PGD. This is mainly due to a lack of confidence in the low-frequency content of recorded ground motions, leading to lack of confidence in peak ground displacements and consequently PGD, with which the latter is strongly correlated (Malhotra 2015). PGD is the combination of three main contributions: (1) coseismic, (2) landslide, and (3) liquefaction displacements. This study neglects coseismic effects because, in general, they do not present localized differential displacements that can influence pipelines; obviously, additional  studies are necessary to better understand if specific conditions affect buried pipelines.
This study develops two new geotechnical models for landslide and liquefaction phenomena. These two models have a similar algorithm consisting of three assessment phases: (1) the analysis of the susceptibility to landslides and liquefaction (2) the probabilistic assessment of the triggering, and (3) the estimation of the resulting PGD. The susceptibility analysis is based on empirical failure domains; such domains are functions of the magnitude of the earthquake and site-to-source/epicentral distance and they discriminate whether a geotechnical failure can or cannot occur according to historical observations. The triggering assessment is conducted for those locations indicated as susceptible, and it is performed adopting physical-based or semiempirical large-scale methodologies. Finally, PGD is calculated adopting empirical formulations for those locations in which the geotechnical failure is triggered. For each point of interest, the final PGD is the maximum of the PGD values due to landslide and liquefaction. The two new geotechnical models are described in more detail subsequently.

Landslide Model
The landslide susceptibility analysis is based on the empirical failure domains proposed by Keefer (1984), who defined three empirical domains corresponding to three different landslide typologies [ Fig. 1(g)]. If information on the landslide typologies is not available, the more conservative domain can be adopted [i.e., the Type I line in Fig. 1(g)]. If the coupled magnitude-epicentral distance falls into the domain of potential failure, then the triggering phenomenon is analyzed; conversely, outside the failure domain, the landslide is deemed negligible and the PGD due to landslide is equal to zero.
This study analyzed the triggering only with reference to translational landslides, i.e., landslides displacing along a planar or undulating surface of rupture, sliding out over the original ground surface (Varnes 1978); rotational landslides were neglected because their safety factor (SF) is generally higher than the SF for planar landslides (Ferentinou et al. 2006). The triggering is analyzed with Newmark's method (Newmark 1965). This model treats the landslide as a rigid block that slides on an inclined plane [ Fig. 1(h)] when the critical acceleration a c (i.e., the threshold base acceleration required to overcome the resistance and initiate sliding) is exceeded. For analyses at the regional level, a simple limitequilibrium model of an infinite slope in material having both friction and cohesive strength can be applied; under these hypotheses, the factor of safety is given by where ϕ 0 , c 0 , and γ = effective friction angle, effective cohesion, and soil unit weight, respectively; u = saturation ratio; γ w = unit weight of water; and α and t = slope angle and slope-normal thickness of the failure slab. Parameters t and u are the most difficult to estimate at large scale; this paper studies u as a parameter; in the absence of detailed information, t can be estimated using models that correlate topographic characteristics with slab depth (Saulnier et al. 1997;Tesfa et al. 2009;Catani et al. 2010;Shafique et al. 2011).
Once the SF is computed, the critical acceleration a c can be calculated according to where g = acceleration of gravity. For a given site of interest, if the simulated acceleration exceeds a c , then lock sliding initiates, the landslide is triggered, and the PGD can be calculated be means of a proper empirical formulation [ Fig. 1(i)]. This study adopted the formulation proposed by Saygili and Rathje (2008), presenting the lowest predictive error [Eq. (10)]. The regression coefficients adopted are a 1 ¼ −1.56, a 2 ¼ −4.58, a 3 ¼ −20.84, a 4 ¼ 44.75, a 5 ¼ −30.5, a 6 ¼ −0.64, and a 7 ¼ 1.55 Liquefaction Model Several liquefaction susceptibility domains are available in the literature for different geographical regions (Papadopoulos and Lefkopoulos 1993;Papathanassiou et al. 2005;Martino et al. 2014). This study conducted the liquefaction susceptibility analysis according to the three failure domains proposed by Galli (2000), based on Italian data and suitable for the subsequent case study. The envelope of the three domains in Galli (2000) was adopted [i.e., envelope line in Fig. 1(j)]. For a given seismic event and for a specific site, the triggering analysis can be carried out only if the magnitude and the epicentral distance fall in the failure domain.
The triggering analysis [ Fig. 1(k)] was performed through the well-consolidated semiempirical method proposed initially by Seed and Idriss (1971) and subsequently improved in several studies in order to create a shear-wave velocity-based liquefaction triggering assessment (Goda et al. 2011;Kayen et al. 2013). Goda et al. (2011) presents an example of application of such a methodology to the large scale. This paper adopted the probabilistic formulation proposed by Kayen et al. (2013). Specifically, the safety factor for liquefaction can be calculated as where CSR= cyclic stress ratio; and CRR = cyclic resistance ratio. According to Seed and Idriss (1971), the CSR at a given depth z is given by where PGA is the result of the simulations discussed in section "Seismic Hazard Model" [ Fig. 1(f)]; σ 0 v ðzÞ and σ v ðzÞ = vertical effective and total stress, respectively, at depth z; MSFðMÞ = magnitude scaling factor, defined by Youd et al. (2001) as ðM=7.5Þ −2.56 ; and r d ðzÞ = shear-stress reduction factor, given by Kayen et al. (2013) according to the following expression: where V s;12 = average shear wave velocity in the upper 12 m of the soil column, obtained from the available V s;30 through the methodology proposed by Boore (2004). According to Kayen et al. (2013) the CRR is given as where Vs 1 = stress-corrected shear wave velocity defined by where P a = reference stress equal to 100 kPa. In Eq. (14), FC is the percentage of fines content, Φ −1 ð·Þ indicates the operator inverse of a cumulative normal distribution, and P L is the liquefaction probability term. Several formulations of P L are available in the literature (Juang et al. 2002;Kayen et al. 2013); this study adopted the PGA-based formulation proposed by Santucci de Magistris et al. (2013) because it is based on Italian data. It is then possible to conclude the triggering analysis calculating the liquefaction potential index (LPI), which is an indicator that succinctly captures the potential liquefaction damage to subsurface structures. LPI was originally proposed by Iwasaki et al. (1978); it is calculated according to Eq. (16), where wðzÞ ¼ 10 − 0.5z and SF is set equal to 0 if SF is greater than 1 Once the liquefaction potential is computed, then the liquefaction-induced PGD required for the risk and loss analysis can be calculated [ Fig. 1(l)]. Toprak and Holzer (2003) identified specific values of LPI corresponding to specific damage scenario. If LPI is less than 5, then the liquefaction is not triggered and the PGD due to liquefaction is equal to 0. If LPI is between 5 and 12-15 (this study assumed 13.5 as a boundary limit), sand boiling and consequently vertical settlements are expected, and therefore PGD is equal to the vertical settlement. Finally, if LPI is greater than 12-15, lateral spreading also is expected to occur, and the PGD is calculated as the square root of sum of square (SRSS) of the vertical settlements and the lateral spread displacements.
The vertical settlement due to sand boiling can be computed according to Tanabe and Takada (1988) from Eq. (17), where H is the thickness of the liquefiable layer and N is the number from the standard penetration test. In the absence of detailed information, N can be empirically correlated with the friction angle (Schmertmann 1978) The lateral spread is calculated as proposed by Bardet et al. (1999) using Eq. (18), where the regression coefficients are

Three-Phase Geotechnical Model
The three-phase geotechnical approach for landslides and liquefaction assessment brings together well-established methodologies with recently developed probabilistic formulations specialized for the Italian geographical context. Because of the modularity of the approach, if new models are available for each phase or for different geographical contexts, the procedure can be further improved and adapted to other areas. Such versatility of the procedure is one of the main novelties of this work, and it fits perfectly within the performance-based engineering framework and its general principle of deconstructing the risk problem.

Risk Assessment
For each of the five correlation cases (i.e., C1-C5) presented in section "Seismic Hazard Model," PGA and PGV were simulated for all the pipeline segments [ Fig. 2(a)]. The geotechnical approaches for landslide and liquefaction allowed the calculation of the PGD (i.e., the maximum among the two geotechnical phenomena). The three IMs were then convoluted with vulnerability and fragility models. Specifically, vulnerability models related directly the IM with the losses (e.g., number of breaks and leaks), whereas fragility functions related the IM with the probability of exceeding a specific level of damage. This study adopted empirical vulnerability and fragility models from the literature. Vulnerability functions [ Fig. 2(b)] were convoluted with the hazard to obtain loss curves [ Fig. 2(c)], i.e., curves representing the probability of reaching or exceeding a specified level of loss; instead, fragility functions [ Fig. 2(d)] were adopted in a nontraditional manner to build damage maps [ Fig. 2(e)] aimed at supporting the infrastructure manager in deciding where to deploy inspection teams to verify potential damage to the system in the aftermath of an event.

Vulnerability Functions and Loss Curves for Natural Gas Pipelines
The HAZUS manual (FEMA 2004) identifies two damage states (DS) for pipelines: (1) leak and (2) break. Such damage states can be induced both by transient actions, also known as strong ground shaking (SGS), and by ground failure (GF) phenomena, such as earthquake-induced landslides and liquefactions. According to HAZUS, if the damage is induced by GF, the percentage of leaks and breaks is estimated as 20 and 80%, respectively. Conversely, if the pipeline is damaged by SGS, the percentage of leaks and breaks is reversed, to 80 and 20%, respectively. The American Lifeline Alliance (ALA 2001) provides damage functions for buried pipes that take into consideration different damage sources (i.e., SGS and GF), materials, diameters, and joint typologies. Such functions permit the calculation of the repair rate (R R ) expressed in terms of normalized length (1=km) as follows: where, for example, K 1 and K 2 for iron continuous pipelines are both equal to 0.15. Eqs. (19) and (20) correspond to SGS and GF, respectively; therefore, according to HAZUS indications, the expected number of leaks and breaks along the infrastructure can be calculated according to Eqs. (21) and (22), where L i is the length of the ith segment of the infrastructure i ð0.8 · R R;SGSi þ 0.2 · R R;SGSi Þ · L i ð21Þ The empirical complementary cumulative distribution function of the numbers of leak and breaks [ Fig. 2(c)] obtained from the simulations of IMs will be the first result available to the infrastructure managers. This is a novel result available to infrastructure managers who typically make use only of the PGA shakemaps provided by national geological survey services (e.g., United States Geological Survey or Italian National Institute of Geophysics and Volcanology). Each simulation produces the values of PGA, PGV and consequently PGD (through the three-step geotechnical module) for each pipeline segment (Fig. 2). These values are used in Eqs. (19) and (20) to obtain the repair rate for SGS and GF. Repair rates are transformed into breaks and leaks according to Eqs. (21) and (22) and then summed for the entire network considered. The simulation results are then represented as empirical complementary cumulative distribution function of breaks and leaks over the whole portion of the network considered, representing the probability of exceeding a specific number of breaks or leaks, respectively.
If the repair cost for each damage typology is available, loss curves can be derived by multiplying Eqs. (21) and (22) by the respective repair costs. As an example, the HAZUS manual (FEMA 2004) suggests a repair cost for a break of US$1,000 for a welded steel pipe with gas-welded joints (category NGP1). The literature provides more recent figures (e.g., INGAA 2016). For each simulation, the total number of breaks and leaks can be transformed into costs and the empirical complementary cumulative distribution function for losses can be provided.

Fragility Functions and Damage Maps for Natural Gas Pipelines
Several fragility curves are available in the literature (Piccinelli and Krausmann 2013); this study adopted models proposed by Lanzano et al. (2014) because of their large and up-to-date underpinning database. Lanzano et al. (2014) presented lognormal and logistic fragility functions for different joint typologies and for SGS and GF corresponding to three damage states: DS0, corresponding to no damage; DS1, corresponding to longitudinal and circumferential cracks and potential compression joint breaks; and DS2, for tension cracks along continuous pipelines and joint loosening in segmented pipelines. Fragilities for SGS and GF are functions of PGV and PGA, respectively. In general, fragility functions are convoluted with the hazard to obtain the risk of the system; this study adopted fragilities to draw damage maps. In particular, for each scenario simulation, i.e., for each PGA and PGV simulated in the hazard module for the specific pipeline segment, it is possible to calculate the probability of reaching DS0, DS1, and DS2 for both SGS and GF according to Lanzano et al. (2014) [Fig. 2(d)]. The specific segment of the infrastructure is characterized by the damage state having the highest probability of occurrence. The simulations provide the frequency of damage for each segment of pipeline, and the modal value among all the simulations identifies the expected DS for the specific segment, resulting in a map of the expected damage states on the infrastructure [ Fig. 2(e)] that can be used as tool to prioritize the inspections in the aftermath of an earthquake.

Critical Notes on Risk Assessment
The preceding risk assessment represents a cutting-edge enhancement of alternative literature methodologies for pipeline networks (e.g., Esposito et al. 2015) because it includes, for the first time, different correlation and cross-correlation models between IMs, and it adopts two three-phase mirrored geotechnical approaches for landslides and liquefaction. Moreover, the adoption of fragility functions to generate damage maps is an innovative adaptation of the fragility tool to the DSS context.
In its current development stage, the methodology includes the treatment of uncertainties related to the hazard and the vulnerability through the use of Monte Carlo simulation based on the GMPE and the fragility curves, respectively. The three-step geotechnical model and the costs do not include the implementation of uncertainties at this stage.
To demonstrate the efficiency of the methodology and quantify the significance of the novel aspects introduced, the following sections present a detailed case study.

Italian Gas Pipeline Network
The Italian gas pipeline network is formed by the national transmission network (spanning approximately 8,800 km) and by the regional distribution networks, which extend more than 22,600 km. Fig. 3(a) shows the Italian transmission network and the regional distribution network. The first plays a vital role because it is the backbone of the gas flow imported from abroad; in particular, gas is injected in the national network systems through eight terminals, five import points [Gela, Gorizia, Mazzara del Vallo, Passo Gires, and Tarvisio; Fig. 3(a)] and three regasification stations [Cavarzere, Livorno, and Panigaglia; Fig. 3(a)]. Pipelines belonging to the national network are characterized by large diameters, buried depth greater than 90 cm (on average 1.5 m), and nominal and operating pressures of 70 and 50 kPa, respectively (Vianello and Maschio 2014).
This study focused on the national network crossing the Friuli Venezia Giulia region [shaded area in Fig. 3(a)]; this region is simultaneously strategic because of the presence of the Tarvisio and Gorizia import points, and also an area of high seismic activity. Fig. 3(b) shows the portion of national gas pipelines network analyzed in this study. The total pipeline length is about 510 km, and diameters range from 26 to 48 in. (about 660-1,220 mm).
The geometry of the system is realistic and is provided publicly by the company managing the network, Snam Rete Gas. For the purposes of sampling, the pipelines were discretized into 10,975 segments having an average length of about 50 m. Fig. 3(b) shows also the topography of the region. The connection from Tarvisio takes advantage of a valley to cross the Alps; from the Tarvisio input point, the pipeline network has three parallel branches, and from Gorizia it has only one branch. The higher number of parallel branches increases the redundancy and therefore the resilience of the system (Golara and Esmaeily 2016).

Friuli Earthquake
To assess the risk of the case-study network, the Mw 6.4 earthquake that occurred in Friuli on May 6, 1976, was used as the base scenario. This event was used because it was the strongest instrumentally recorded event in the northern Italy region (Moratto et al. 2011). Ten seismic stations recorded the event across Italy, and measurements are freely accessible ;   Fig. 3(b) represents four of the ten stations as triangles. Fig. 3(b) also depicts the fault plane as defined by the Italian database of individual seismogenic sources (DISS Working Group 2018). For the Monte Carlo simulation, this paper adopted the GMPE proposed by Bindi et al. (2011) [Figs. 1(d-f)] as the most suitable and recent GMPE based on the Italian event database. Bindi et al. calibrated an intraevent standard deviation for PGA and PGV of 0.29 and 0.27, respectively.

Available Data
Several maps were collected for the case study and a geodatabase was assembled in a geographical information system (GIS) framework: first, a 30-m resolution digital elevation model for the analyzed region [ Fig. 3(b)] was obtained from the shuttle radar topography mission (SRTM) project (Farr et al. 2007). Then, based on DEM data, the map of slopes [ Fig. 3(c)] was calculated. For V s;30 [ Fig. 3(d)], the USGS Global Vs30 Map Server was used (Wald and Allen 2007) in combination with the soil map adopted by the Italian Institute of Geophysics and Volcanology (INGV) (Michelini et al. 2008). Finally, maps of soil friction angle synthesized on the basis of the 1:500,000 Italian geological map (MEPLS 2018) and the 1:150,000 regional geolithological map (IRDAT-FVG 2018). This information needs to be gathered in advance for the area of interest in order for the postearthquake risk assessment tool to be ready when an event occurs. A fair amount of information is often available to the network manager, and it can be used for the set-up of the risk assessment tool.

Results
To obtain stable low/high percentiles of final losses, 10,000 simulations are needed. The entire simulation procedure and the generation of the loss curves and the damage maps can be performed in less than twenty minutes, and therefore it is very suitable for a postevent DSS. Focusing on the cost estimates of Fig. 4, from a geotechnical point of view, an increase in soil saturation (i.e., from 50 to 100%) leads to an increase of the expected losses because of its significant influence on both landslide and liquefaction. Therefore, it is very important to have detailed information about the soil saturation along the infrastructure route (e.g., Lillesand et al. 2014).

Loss Curves
From a statistical point of view, neglecting the correlation (Model C1, Fig. 4) leads to large underestimation of the final losses, as also observed by other studies for other typologies of infrastructures (e.g., Miano et al. 2016). Moreover, when the spatial correlation and the cross-correlation between IMs are taken into account, the dispersion of losses tends to decrease with respect to C1, leading to a more accurate estimation. In this case study, the results obtained considering the spatial correlation and neglecting the updating (i.e., Case C2) tended to slightly underestimate the losses with respect to the case which included updating (C4). Differences between Cases C3 and C5 were negligible; when both spatial correlation and cross-correlation were considered, the further updating with real data did not improve the final assessment. Thus, neglecting the cross-correlation of the intensity measures may lead to an underestimation of postearthquake losses. Fig. 5 shows the maps of potential damage induced by SGS and GF. The predominant damage state is DS1 and the localization of the damage is mainly in the projection of the fault plane (the rectangle in Fig. 5).

Damage Maps
In this study, it was extremely rare to attain DS2 for such types of infrastructure under earthquake actions, as also emphasized by the European Gas Pipeline Incident Data Group in their annual report on gas pipeline incidents (EGIG 2015). Nonetheless, regional networks can be more susceptible to earthquakes and potentially can incur major damage (DS2). The same map also plots historical observations of landslides and liquefaction, according to the Italian Catalogo italiano degli Effetti Deformativi del suolo Indotti dai forti Terremoti (CEDIT) database (Martino et al. 2014). All the geotechnical failures historically recorded were very close to pipelines locations characterized as damaged (DS1) according to the methodology. In the future, the proposed damage maps can be integrated with satellite surveys of geotechnical failures [which are increasingly available nowadays and are often adopted in postearthquake recognition (Sextos et al. 2018)] in order to increase the level of confidence in the localization of potential damage.

Conclusions
This paper presented a methodology to quantify losses in terms of repair costs, leaks, and breaks and to localize the damage of a gas transmission infrastructure in the first critical period after a major earthquake. The methodology is composed of a Monte Carlo simulation-based procedure for the generation of plausible intensity measure fields that are spatially correlated and cross-correlated; two new back-to-back three-phase geotechnical approaches accounting for (1) susceptibility, (2) triggering, and (3) ground dislocation due to landslide and liquefaction; and a risk approach based on two different literature vulnerability and fragility models. The application to a real case study verified the coherence and computational efficiency of the methodology, and led to the following conclusions: 1. Considering spatial correlation and cross-correlation of the intensity measures required to assess the seismic risk of natural gas pipeline networks leads to a more reliable estimate of losses, avoiding the potential underestimation that was observed when the intensity measures were taken as fully uncorrelated or simply spatially correlated; 2. The updating of the intensity measure fields from accelerometric stations in real time, if available after the earthquake, is beneficial because it allows anchoring the calculations to real observations; 3. The methodology presented herein provides the infrastructure stakeholder with a breakdown of costs and type of interventions required in case of damage in the aftermath of the earthquake event; and 4. This methodology also produces maps localizing the expected damage, thus facilitating inspection missions and reducing the recovery time in the aftermath of a seismic event. Further study is needed for the implementation of coseismic effects considering stochastic simulation of the fault slip, as recently suggested by Goda (2017). New vulnerability models are also required to better reflect the damage typology, potentially using machine learning algorithms. This work can be extended combining the damage maps with (nearly) real-time data on landslides that can be retrieved through satellite imagery and remote sensing capabilities.