Communication Network Reliability Under Geographically Correlated Failures Using Probabilistic Seismic Hazard Analysis

The research community’s attention has been attracted to the reliability of networks exposed to large-scale disasters and this has become a critical concern in network studies during the last decade. Earthquakes are high on the list of those showing the most significant impacts on communication networks, and simultaneously, they are the least predictable events. This study uses the Probabilistic Seismic Hazard Analysis method to estimate the network element state after an earthquake. The approach considers a seismic source model and ground prediction equations to assess the intensity measure for each element according to its location. In the simulation, nodes fail according to the building’s fragility curves. Similarly, links fail according to a failure rate depending on the intensity measure and the cable’s characteristics. We use the source-terminal, and the diameter constrained reliability metrics. The approach goes beyond the graph representation of the network and incorporates the terrain characteristics and the component’s robustness into the network performance analysis at an affordable computational cost. We study the method on a network in a seismic region with almost 9000 km of optical fiber. We observed that for source-terminal that are less than 500 km apart the improvements are marginals while for those that are more than 1000 km apart, reliability improves near a 30% in the enhanced designs. We also showed that these results depend heavily on the robustness/fragility of the infrastructure, showing that performance measures based only the network topology are not enough to evaluate new designs.


I. INTRODUCTION
Telecommunication networks have become increasingly ubiquitous during the last two decades. By 2013 the importance (for regulators, consumers, providers, and researchers) of studying the risk in these networks was pointed out in [1]; lately, this discussion has become even more relevant. During the 2020 pandemic, we relied more on the telecommunication The associate editor coordinating the review of this manuscript and approving it for publication was Jesus Felez . network to carry out our activities. After the initial infection waves and particularly after the wide vaccination effort many activities have returned to being in person, while many others remain online. The expectation is that the tendency to increase online activities (which was present before the pandemic) will continue into the future. This tendency has been boosted by the pandemic. For example, in e-business, the changes in consumers' behaviors are expected to last [2]. This dependence means that more of our daily tasks rely on telecommunications networks, and the lack of their VOLUME 11, 2023 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ availability impacts harder on every activity in our society [3]. Consequently, the study of reliability in networks continues to be of utmost importance and additionally faces new challenges. Network reliability had its beginnings in the second half of the twentieth century when the problem of building reliable circuits was addressed. Then came the need to analyze the reliability of networks deployed in large geographical areas such as telegraph, telephone, and data networks see [4] and [5] for a historical overview. In 1983, Ball and Provan [6] established that network reliability computation was an NP-hard problem, even in a setting with several simplifications. Nevertheless, the network reliability analysis is of practical necessity; therefore, researchers developed methods to overcome computational complexity. For the last 60 years, much research has been accumulated on this topic, with techniques from different areas such as graph theory, Monte Carlo simulation, and complexity, to mention a few.
As telecommunication networks become more widespread, higher service levels are needed. Coupled with these new service standards, interest in how large-scale disasters affect network reliability is becoming more critical in the study of reliability, see [7], [8], [9], and [10]. As mentioned in [11] the risk assessment of networks must go beyond the graph representation and account for the geographical region where they are embedded. Network failures that arise from large-scale hazards such as earthquakes, wildfires, or extreme weather events show geographical failure correlations-outages due to this type of hazard present different behavior than a typical outage. Natural disasters provoke multiple simultaneous faults, failed components have spatio-temporal correlations, and are harder to repair (see [8], [12]). The strong correlations these faults show make them unsuited for analysis under the independent failures assumption. If independence is assumed, it tends to overestimate the reliability [13] and could lead to less reliable designs [14]. As communication networks are undoubtedly an essential part of the critical infrastructure, their proper functioning in post-disaster periods is crucial. The network vulnerability study is combined with disaster-resilient network design, restoration methods, and emergency network deployment to attain communication during the recovery period [10], [15].
Among geographically correlated failures, we focus on earthquakes for two reasons: the damage they produce is an economic burden even in countries prepared for them, and there are well developed methods to assess the infrastructure risks (see [16]). In [8] the worst catastrophe in Japan's history is discussed. This was in March 2011, consisting of the Great East Japan Earthquake of magnitude M w = 9, together with the subsequent tsunami and nuclear accident, it damaged around 2,200 telecommunication buildings. Similarly, as a consequence of the magnitude M w = 8.8 earthquake in southern Chile in 2010, the Chilean telecommunications network was severely affected by a one-day outage, as reported in [17]. Regarding Internet service, Chile was practically cut off from the outside world. A later study of network performance showed that the Chilean Internet downtime on the day of the earthquake was close to 68%, meaning that only 32% of Chilean networks were accessible from the rest of the world.
Earth science, structural engineering, seismology, and many other areas have studied seismic hazard and risk analysis. Several research centers worldwide are dedicated to developing new methods to better prepare us for a seismic event. A natural question therefore is: Can these models be adapted for the evaluation and design of telecommunications networks?
The telecommunication network was modeled as a graph where nodes represent data exchange facilities and edges correspond to optical fiber lines. To analyze the network's performance, the network's probability of remaining connected after the seismic event was estimated. In our model, failures can occur in both edges and nodes. We consider that all the components can assume only one of the two states: operational or not. For nodes, following the recommendations of the HAZUS technical manual for earthquakes [18], given the seismic characteristic, we compute the intensity measure at the node location and, given according to the structure fragility, we assign a failure probability. In the case of edges, we assume the optic cables are buried underground. We compute a failure rate for each link based on the seismic intensity at the edge geographical location, the type of terrain, and the edge length according to the methodology proposed in [19]. We use the source-terminal and the diameter constrained reliability metrics. Whether the post contingency network is declared operational or not is based on the connectivity (see [20] for details on this and other alternative criteria for defining operational states).
We compute a failure rate for each link based on the seismic intensity at the edge geographical location, the type of terrain, and the edge length according to the methodology proposed in [19]. We use the source-terminal and the diameter constrained reliability metrics. Whether the post contingency network is declared operational or not is based on the connectivity (see [20] for details on this and other alternative criteria for defining operational states). To be precise, this work differs from [19] in various aspects, including the following ones: our analysis focuses on a more extensive network. Because of this, we used connectivity-based reliability as a performance measure instead of a capacity analysis. Finally, we use an alternative Ground Motion Prediction Equation adequate to the studied network's geographical location.
As a case study, we use the Chilean National University Network (REUNA). REUNA's high-speed optical network extends 3000 km from the north to the south of Chile (Arica to Puerto Montt). It privately connects all affiliated institutions to facilitate the exchange of data (traffic) between them. Using the methodology proposed in the current work, we evaluate the base REUNA network topology, as well as alternative network designs proposed in [21] to improve the network fault tolerance. We estimate different reliability metrics and we discuss the characteristics and relative advantages of the proposed designs.
The contribution of this work is two-fold. First, we integrate methodologies, which although standard, have not been used together before in the simultaneous assessment of node and edge failures according to the probabilistic seismic hazard analysis (PSHA) method and the ideas proposed in [19], which are the state-of-the-art methods respectively for infrastructure failure and underground cabling failure when affected by an earthquake. The PSHA allow us to incorporate the networks' geographical characteristics into the network performance analysis, complementing other methods that emphasize different aspects of the network. Second, the reliability metrics, based on the connectivity, allow us to quantify the contribution of additional components to the robustness of the network, and to identify those pairs of nodes that benefit the most from the additional elements.
The rest of the paper is organized as follows. In Section II we show previous works that have served as a foundation for this study. In Section III we describe the model used for the fiber optic network studied, in Section IV we present the numerical results obtained, and finally in Section V the conclusions of the research.

II. RELATED WORKS
The reliability of networks exposed to large disasters has captivated the research community's attention and has become a critical concern in network studies during the last decade (see surveys [8] and [22]). Looking at the damage caused by recent disasters, earthquakes top the list of those showing the most significant impacts on communication networks, and simultaneously, they are the least predictable events [12]. Telecommunication network reliability specifically under seismic failures has recently been addressed in the literature. We discuss those related to the one presented in this article.
In [21] the authors propose a method to attain a robust network design with seismic events in mind. Their approach is based on graph theory concepts such as bridges, cut-points, and edge betweenness centrality, among other metrics. The authors proposed an algorithm to solve the problem of adding a restricted number of links to a given topology, in such a way that it becomes biconnected and optimizes each metric of interest. An overview of algorithms available for the enumerations of shared risk link groups (SRLG), can be found in [23]. This method captures the dependency between failures by grouping network components that are prone to fail together. This approach has been used to study the network resilience against earthquakes. In [24] the authors used SRLG to study the source-terminal reliability of two Japanese networks using historical earthquake data. In this setting any component that is near enough to the earthquake epicenter fails. In [25] the SRLG were used to study the vulnerability and availability of the Italian portion of the Interoute optical fiber network. Using seismological research and historical data, they employ a seismic source model that allows for estimation of the earthquake activity rates for a given epicenter and magnitude. Again, any components in a disk around the epicenter fail but the radius of the disk depends on the earthquake's magnitude. A set of SRLGs can be very useful for network recovery/protection mechanisms but for network design it cannot account directly for the effect of improving the seismic infrastructure for the components. As mentioned in [16], a failure model should consider that optical fiber cables fail with a probability which depends on the length and on its distance to the epicenter, but it should also depend on the component characteristics.
In the present work we use the probabilistic seismic hazard analysis (PSHA) method to study the reliability of an optical network deployed in a seismic region. In this method the post contingency state of each component is computed through a three-stage process, we follow the description in [26]. The first stage is to obtain the earthquake characteristics from a seismic source model (epicenter, depth, and magnitude). The second stage corresponds to obtaining the intensity measure of each component of the network, we will compute the spectral acceleration (SA) for the nodes and the peak ground acceleration (PGA) for each portion of a link. Finally a component will fail randomly with a probability that depends not only on the intensity measure but also on the characteristics of the components. The advantage of using PSHA method is that the computation of the intensity measure of each component allows for consideration of site conditions where it is located, it also incorporates the characteristics of the component, such as the height of the building in the nodes case and the material of the optical fiber in the link case. By giving the possibility of considering how these component characteristics will impact the reliability of the network this approach differs to those previously mentioned. Specifically, this approach allows for the consideration that more robust components will tend to fail less.
Although PSHA is a well established method, with a good availability of the needed parameters, it has rarely been employed for network reliability evaluation. One exception in a related area is for power networks [27], the focus in this work is to show the impact of geographical correlations on a power network design. Another exception is for a complex transportation network in [28]. Here, the authors consider the San Francisco Bay Area transportation system; the network performance metric is the accessibility to the users, which is only estimated for a subset of scenarios due to the computational cost. These works do not focus on communication network reliability, usually expressed by a connectivity measure.
One reason that can explain why this seismic risk assessment hasn't been applied earlier to telecommunication networks is because, although there are numerous fragility functions for predicting damage for buildings, fewer fragility functions exist for infrastructure systems, and in particular, almost none for buried cables. In [29], the authors produced some of the first empirical failure rate functions for buried (electrical) cables with respect to ground shaking and liquefaction-induced ground deformation. In 2021 this work was adapted for optical fiber cables (see [19] with a preliminary version in [30]). In that work, an estimation of the failure rate (FR) per km is proposed, where the FR depends on the earthquake intensity measure (PGA) and the site characteristics (landslide and the liquefaction potential), and optical fiber vulnerability. For nodes (data exchange facilities) their failure/working state is computed using the fragility curves given by the characteristics of the building. The authors focused on the traffic loss after the seismic event, to analyze the RIMAC network (a small 4 nodes optical network). The contribution of this work is remarkable, being the first to extend the PSHA method to buried optical fibers. We will follow this methodology to assess the risk for the optical fiber link in our setting.
The present work is the first to evaluate the reliability, given by the connectivity where failure states are simulated with the PSHA method for nodes and links. This performance measure is better suited for bigger networks. In this study the reliability performance measure is applied to a real network with 18 nodes and 19 links with almost 9000 km of optical fiber located in a seismic region connecting, Chile, one of the longest countries in the world (see Fig. 1).
First, we discuss that the methodology shows coherent results with other methodologies employed for the same network. Then, we show that the reliability metric allows for obtaining a more insightful impact for the proposed designs, by distinguishing the source-terminal reliability for nodes that are farther apart. Finally, we show that the method allows for quantifying the impact of improving the seismic preparedness of a component on the whole network reliability.

III. MODEL AND IMPLEMENTATION
We consider a telecommunications network subject to failures due to seismic events. A graph with nodes and edges represents the telecommunication network. The nodes represent the buildings that host data exchange facilities for the correct operation of the network (without considering intermediate components), while the edges represent the optical fiber cabling that we assume is buried. Only the faults that are explained in the following subsections are considered. The optical fiber cables were discretized in portions of 10 km, to estimate a failure rate for each portion separately.

A. THE PROBABILISTIC SEISMIC HAZARD ANALYSIS METHOD
The PSHA is a framework to assess seismic hazards and risks [26]. The word ''Probability'' in the name emphasizes that the approach incorporates uncertainties in the analysis. The elements of the method can vary to consider different users' needs or to be region specific. We first give a general description of the method, and later we review the details of the version used in this study.
The PSHA method allows for the simulation of each component's state in a three-stage process. In the first stage, a seismic source model estimates the magnitude, location, and other seismic characteristics; the model uses historical seismic data and/or geological models to calculate the occurrence rate. The second stage uses a ground-motion prediction equation (GMPE) to simulate the resulting intensity measure at each location based on the characteristic of the chosen seismic event (in our case spectral acceleration). GMPE accounts for the component's height and the site conditions. In the third stage, the probability that the infrastructure suffers a specific damage level is computed through the fragility. More precisely, the infrastructure damage functions take the form of lognormal fragility curves that relate the probability of exceeding the component's damage level for a given intensity measure, in general, the Spectral Acceleration (SA) although related inputs could be also considered such as the seismic duration or liquefaction (see the discussion in Ch. 9 from [26]).
The integration of the three-stage process gives the probability that an infrastructure suffers certain damage considering the annual seismic rates from different sources, the different intensities, the site characteristics, and the building type.
In this work, we show the potential of the PSHA method combined with the recently developed risk assessment for buried cables for the study of telecommunications network reliability. In the second and third stages of the method it is necessary to take into account the specific vulnerability levels of the components. So we focus on these two stages and consider a straightforward seismic source model.
The GMPEs have significantly evolved in the last decades. Countries exposed to seismic events have installed seismometers that form a network to collect detailed data on the seismic impact. In particular, in Chile, the Integrated Plate Boundary Observatory Chile (IPOC [31]) has more than 100 stations, generating a strong ground motion dataset of more than 400 seismic events since 1985. Moreover, it includes three large-magnitude interface earthquakes that struck Chile in recent years (2010 M w 8.8 Maule, 2014 M w 8.1 Iquique, and 2015 M w 8.3 Illapel). The study and analysis of these well-characterized seismic events lead to improved GMPE, attaining smaller errors and increasing the validity of the prediction to higher magnitudes.
As the REUNA network is located in Chile, we used one of the most recent calibrated GMPE, proposed in [32]. The applicability of this GMPE is reasonable for a range from magnitudes M w 5.0 to 9.0 and distances from the epicenter up to 1000 km.
The Spectral Acceleration, SA is a magnitude which depends on the period T , and such that ln(SA(T )) follows a normal distribution with the median given by Equation 1 and the variance given by Equation 2. Following the notation in [32] for a given period T (in seconds) the median of the SA logarithm is: where M w is the moment magnitude of the earthquake; Z h is the hypocentral depth in kilometers; R is the source-to-site distance, V S30 is the average shear-wave velocity between 0 and 30-meters depth, ⃗ θ are the regression coefficients that in this model depend on T .
The total variance of the model σ 2 (T ) is given by the sum of the between-earthquake residual variance τ 2 , singlestation residual component variance φ 2 SS , and into a site-tosite variance φ 2 S2S : For a detailed discussion on the model and the parameters, we refer the reader to [32]. The precise formula and parameters used are reproduced in the Appendix, Section A. These equations allow us to sample the intensity measure for each component of interest, that is, the Peak Ground Acceleration (PGA) for links and the Spectral Acceleration (SA) for the nodes.
Except for V S30 , all the other parameters depend on T , see Table 1 in [32]. The V S30 allows for characterizing the seismic amplification of the ground. In the main cities, this parameter is available from detailed studies. Table 1 shows the references to obtain the values for each node. For the rural areas of Chile, in [33] estimations for V S30 are given, which were used for the links and for those nodes where a detailed study was unavailable.

B. NODES RISK ASSESSMENTS
For the third stage on the nodes, the failed/working state is simulated based on the building characteristics of the components. In the communication substation case, we used fragility curves that account for the robustness of the building characteristics. Given an intensity measure, the fragility curve allows for computing the probability that the building has suffered a certain damage level. The fragility curves give the vulnerability of infrastructure according to the structural typology. In our case, we assume two types of buildings: one built with concrete shear walls with a HAZUS C2L classification, and a second one with a HAZUS W2 classification, corresponding to buildings more susceptible to damages. Although we are working with this specific structure, this approach is perfectly adaptable to any other type of infrastructure classified by the HAZUS manual. The parameters used to generate the fragility curve are summarized in Table 2 obtained from the HAZUS manual. In C2L type buildings a complete structural damage is observed when the structure has collapsed or is in imminent danger of collapse due to failure of most of the shear walls and failure of some critical beams or columns. In W2 type buildings a complete structural damage is observed in some of these situations: large permanent lateral displacement, collapse or imminent danger of collapse due to failed shear walls, broken brace rods or failed framing connections; fall of foundations; large cracks in the foundations. Fig. 2 shows the fragility curves for both types of infrastructures; it is easy to see that the complete failure probability is higher for W2 than for C2L, showing that the W2 structure type is more vulnerable than C2L.

C. LINKS RISK ASSESSMENT
After an earthquake severe ground deformation occurs, which will produce damage by bending cables or shearing due to the ground displacement caused by the earthquake as shown in [41]. The failure probability of a link increases if links are VOLUME 11, 2023 longer, and also if any link portion is closer to the seismic source. This behavior is captured in the approach developed in [19] that will be used. For completeness we reproduce the details here.
We start by computing a failure rate FR l for a small portion of buried cable l that will depend on the local characteristics and then we integrate all the small portions that compose the link. This rate depends on the cable robustness and the peak ground displacement (PGD) 1 which is an earthquake intensity measure, see Equation 6 below. Depending on the model PGD depends on PGA and other terrain characteristics. Models to estimate PGD are calibrated for specific geographical areas. For our study we picked the one proposed by Saygili and Rathje in [42] given in Equation 5, because it was tested with other PGD estimations for the Chilean seismic records in [43], and showed a good fit. Then, aside from those terrain characteristics considered for the SA, we also need the critical acceleration k c . The k c parameter describes how much horizontal acceleration is needed for a sliding displacement to occur with a certain probability. The critical acceleration k c is a building design criteria see [44] and references there in.
We divide each link in portions of 10 km (or less in the extremes). The (random) number of failures that occurs in a link L is given by a Poisson random variable with parameter computed as sum of the FR l over all of the portions l that compose the link, that is where |l| is the length of the portion l.
If there occurs at least one failure in the link, the whole link will be in a failed state; that is, the link failure probability p L conditional on the seismic parameters is given by We compute the PGA from Equation 1 with T = 0 and together with k c we plug them into Equation 5 to obtain the 1 Note that failure rates are estimated based on observed data arising from repairs, see the discussion in [29].
The resulting PGD is used in Equation 6 together with α, the cable typology robustness parameter, to obtain the failure rate for each portion: Two values of α are adopted to discuss its influence on the network performance, α = 0.26 and α = 1.07, associated to the least and the most vulnerable cable typology, respectively. In our case study in Chile, there isn't a detailed liquefaction study available for all the non-urban areas where the optical fibers are buried, so we assume a single parameter for the whole affected area under consideration.

D. RELIABILITY PERFORMANCE MEASURE
Defining reliability metrics will depend on the context where we want to use them. For example, flow capacity metrics may be of interest for electrical or transport networks, to name a few, see [45] and references therein. In our case, we use the reliability metrics associated with network connectivity for telecommunications networks. Let G = (V , E) be a graph with vertex set V and edge set E, given vertices s and t in G, the st reliability is the probability that a continuous path exists between the source node s and a destination node t. Since routing decisions are a path from a source node to a destination node, the notion of st reliability has been considered, see [46].
An alternative model of topological reliability is the diameter constrained reliability of a network G; this measure considers not only the underlying topology, but also imposes a bound on the diameter, which is the maximum distance between nodes in the post-contingency network. The (st − DC) diameter constrained reliability with diameter d is the probability that the nodes s and t remain connected by paths composed of d or fewer working edges, see [20]. This reliability metric becomes more relevant for services that require synchronous communication such as telehealth or online education.
We will focus on these two reliability metrics but the methodology could also be used with other reliability measures based on connectivity criteria.

E. SIMULATIONS
The historical information of Chilean earthquakes is used as a seismic source model to generate seismic events' characteristics randomly. In particular, we used a list of 382 earthquakes with M w ≥ 6 which occurred since 1960 that are included in the catalog U.S. Geological Survey (USGS), see [47]. We integrate these models within a Monte Carlo simulation scheme that allows us to estimate the source-terminal connectivity on the post contingency network for each pair of network nodes for each sampled seismic event. We compute the st and st diameter constrained reliabilities by average across the sampled scenarios. We describe below the model in detail.
We propose a four-stage simulation model. The first three stages are performed according to the PSHA method and allow us to obtain the components' working/failure states; the fourth stage is added to compute the source-terminal connectivity of the post contingency network, see Fig 3. In the first stage, as a seismic source model, we randomly choose the characteristics of the earthquake from a list of historical seismic events. Then in the second stage we sample the intensity measures, SA(T ) and PGA, from a lognormal distribution with parameters obtained from Equation 1 and Equation 2. In this step, Equation 1 accounts for the in site characteristics (V S30 ), geographical position (R), and the height of the building directly related to the period (T ).
In the third stage we obtain the failure/working state for nodes and links. For the nodes, following the left-hand side of the flowchart in Figure 3, given a damage level and the SA(T ), a conditional failure probability for each node is computed from the fragility curve. We obtain the state of nodes simulating their failures according to this probability. Similarly,following the right-hand side of the flowchart in Figure 3, for each link portion we simulate the PGA from a lognormal distribution with parameters given by Equations 1 and 2 with T = 0 and R given by the location. The PGD is computed using Equation 5 and the failure rate FR l is obtain from Equation 6. The failure rate for the link is obtained by adding failure rate for the portions on the link (see Equation 3). Finally, we obtain the link state by simulating their failures with probability p L given in Equation 4. Once the nodes and links state is obtained from the simulation, in the fourth stage, the connectivity of each post-contingency graph is analyzed. For each pair of nodes, we perform a shortest path algorithm over the survival network. As the shortest path algorithm assumes nodes remain working, a failed node is incorporated by declaring nonworking any link that is adjacent to it. This is a standard approach.
The simulation process is repeated several times to obtain a number of statistically representative scenarios. A record of the connectivity of each source terminal pair of nodes for each scenario is kept. Then, the method finishes by averaging across all the sampled scenarios. We obtain the st reliability and the st − DC reliability estimations for each pair of nodes. We discuss in detail the obtained results in the following section.

A. EVALUATION OF THE REUNA NETWORK
We apply the developed Monte Carlo methodology to the REUNA network shown in Fig. 4. The Monte Carlo simulation employs a sample size of 10 5 . In each simulation, we generate an earthquake epicenter and depth using the historical Chilean seismic data. As we are interested in events that cause significant damage, we generate seismic events of magnitude eight (i.e, we will be looking at the reliability of the network conditionally in the occurrence of a magnitude eight earthquake). Which, as we mentioned, is an event occurring in Chile approximately once every ten years or more. VOLUME 11, 2023 As this study is a proof of concept, some simplifications were made. Without loss of generality: we assumed that all nodes correspond to infrastructures of the same type (buildings with a HAZUS C2L classification), so that they have the same fragility curves; and we assumed that the length of the links was equal to the distance (Euclidean) between the nodes. For the α parameter associated to the link failure rates (Equation 6) we assume a value associated with more robust optical fiber (α = 0.26). We also assumed a k c value equal to 0.1. Recall that this parameter, which depends on the slope of the construction site, is a building decision, which depends on the construction decisions and investments; we use here a typical value for all the sites (we discuss the sensibility with respect to other possible values mainly in the Appendix B). These assumptions can be easily replaced with real life data (as available) without changes in the computational methodology. The model was implemented using the R language, and the experiments were run on a DELL computer, IntelC,B) Core i5-8250U CPU @ 1.60GHz CfB-4 with 8GB RAM. The total time for running an experiment with sample size 10 5 replications was 210 minutes.
We used two measures of interest; the st and the st diameter constrained reliability with diameter d = 5. In Table 3 we can see some of the results. We considered all pairs of nodes in the network; the table presents the average estimated reliabilities (average over all pairs of nodes), as well as the smallest and largest reliability values observed. We see that the methodology can be applied to obtain, in a reasonable computational time, estimates for the reliability of the network. The results show the importance of choosing the adequate reliability measure, as there are clear differences in the estimated values of the two measures considered. Additionally, they show that it is important to take into account, not just a global network reliability measure, but also the particular reliabilities corresponding to different source-terminal pairs, as some node-pairs have quite different connectivity probabilities from the mean reliability values.
From the experiments, we were able to estimate the marginal failure probabilities for each component and the correlations by pairs of links or pairs of nodes estimated for all the network's components. Table 4 summarizes these results, which show a significant failure correlation. These findings are consistent with the discussed studies.

B. COMPARISON OF DIFFERENT TOPOLOGIES
In Table 5 the following six alternative topologies result from adding three links to the original REUNA network in order to obtain a biconnected topology:  Moreover, we considered four additional topologies R−3A+1, R-3B+1, R-3C+1, R-3F+1, resulting from adding the link between nodes (Iquique -Temuco) to the topologies R-3A, R-3B, R-3C and R-3F. Therefore, we analyze ten new topologies: six with 3 additional links and three with 4 additional links. We employed the methodology proposed in this work to evaluate the reliability of these alternative networks and compare them to the original REUNA network.
The mentioned designs were proposed in [21], where the robustness of the REUNA network was compared against the enhanced designs, in order to find the most robust of them. There the network robustness was measured with the edge betweenness centrality, relative number of cutsets, and node Wiener impact. We will compare these results with the ones obtain from the methodology proposed in our work. The reliabilities discussed in this work were estimated as an average across the simulated scenarios. Given that all variances are bounded by 0.003, with a 95% confidence level, all the estimations have a precision of 0.01 or better.
Tables 6 and 7 display the st and st − DC reliability of the original network as well as the various proposed designs. For the st reliability all three-edges designs have reliability between 0.654 and 0.745 and improve the results of the original network. Analogous for the st − DC reliability in the range of 0.563 to 0.606. This shows that our results are coherent with those obtained in [21]. The proposed designs with three extra edges significantly improve the reliability of the original network. In our analysis among the threeedge designs, model R-3B has one of the best performances. Notice that R-3F has a very good performance while in [21] it had one of the worst. A similar conclusion is obtained for the st − DC reliabilities. As we will detail in the following subsection, an analysis that considers the distance between the source-terminal nodes allows a deeper understanding of the contributions to the robustness of the additional components.

C. QUANTITATIVE IMPROVEMENT OF THE RELIABILITY
The R-3B and R-3C new designs show a 20% improvement when we consider the average st reliability and 22% st − DC reliability which is significant. For the four edge design R-3C+1 shows an marginal improvement with respect to the R-3C for the st (less than 1%). The same results and analysis are obtained for R-3B and R-3B+1. An other interesting observation is that R-3B outperforms all four link designs except for R-3B+1. For the st −DC reliability a 7% improvement is attained with respect to the R-3B and R-3C designs. When we disaggregate with the distance between the pairs of source-terminal nodes, the impact of the three and four edge proposed topologies can be better quantified see Fig. 5 and Fig. 6. For source terminal nodes separated by 1000 km or more the st reliability improves near 30% for the R-3F and near 40% for the R-3B. For the same set of source-terminal nodes improvements attain near the 80% for the st − DC reliability when the R-3B+1 is analyzed. An other interesting fact is that for pairs of nodes that are less than 2000 km apart the best design is R-3B+1 (very similar to R-3B) but for those pairs that are more than 2000 km apart the R-3F+1 significantly outperforms R-3B+1. We plotted the reliabilities of all source-terminal pairs 1001 km apart or more to unveil which pairs benefit more from each design; 48 pairs are between 1000-2000km, and 24 pairs are more than 2000 km apart. Figure 7 confirms that for the st reliability, the R-3B+1 design outperforms R-3F+1. On the other hand, Figure 8 shows that there is an important level of dispersion for st−DC reliability, which also means that it is not clear which design is better, in particular for those pairs of nodes in the category 1001-2000 km; this issue merits further analysis.
This detailed information is valuable for the public policy designers showing how the investment in the telecommunication network will impact different regions and contributes to decentralization.  Tables 8 and 9 summarize the reliability of the networks, showing both the classical st reliability and the st − DC reliability with diameter d = 5. We show the results for the original network, for the R-3F and the R-3B three links designs and the R-3F+1 and R-3B+1 designs. The added links improve the reliability in those nodes that are further apart, giving the perspective of the final user. Our analysis of the networks reliability improvement complements those performed in [21] where the robust measures edge betweenness centrality and node Wiener impact focus in the relevance of critical components.

D. IMPACT OF BUILDING DECISION CRITERIA ON THE RELIABILITY
Robust metrics based on the network topology cannot account for the reliability improvement given by better prepared infrastructure. The approach in this work can directly measure the impact on building decisions. Table 10     It shows how the same topology can attain very different performance with more prepared infrastructure. Moreover, when we analyze the reliability with the more robust infrastructure (k c = 0.4) the best design is R-3F+1 and not R-3B+1. While when we change the nodes' infrastructure to the W2 infrastructure type R-3B+1 remains the most reliable design (see Table 11). The method allows us to perform a similar analysis with other site characteristics. For example for the nodes decision considering multiple locations in the same city would impact the V S30 and different building heights change the period T . Similarly, for the links, changes in the vulnerability of the optical fiber α where the cables are lying will affect the reliability performance.

V. CONCLUSION
We study the reliability given by the connectivity of the telecommunications network when affected by seismic events. We used a network failure model based on stateof-the-art methods for underground cables and infrastructure failures. Although each of them have been previously been studied in the literature, this is the first time they are implemented together to study network reliability.
The modeling framework enables the evaluation of enhancing network reliability by considering the optical cable type, the infrastructure's characteristics, and geography particularities. While at the same time, it is simple and abstract enough to have an amenable computational cost. The method addresses the need of the decision-makers to allocate funds for infrastructure improvements to upgrade the telecommunication networks' performance. It also enables the identification of the users who benefit the most and quantifying the improvement.
We implemented the method through a Monte Carlo simulation and we studied the REUNA network's behavior in Chile. The original network and a number of alternative topologies were analyzed by estimating their reliability. Our results show that some of the proposed alternative network topologies improve the source-terminal reliability by almost double, for pairs of nodes that are farther apart, while there are only marginal improvements for nearby pairs. We compare our results to those obtained in [21], where they study the same network and modification with a different fault-tolerance method. Our results bring complementary insight with the perspective of the end-user to those in the previous work that focused on the component relevance.
As discussed in Subsection IV-D, the proposed method also allows us to understand how much the components prepared to withstand earthquakes enhance the entire system's performance. A prepared link can be obtained from using more robust optical fiber cables or by improving the site conditions reflected in the critical acceleration parameter. The node preparedness is represented by its fragility curve; specifically, less vulnerable nodes will decrease the (conditional) failure probability.
Our study shows that we can incorporate the geographical characteristics and the component's robustness into the network performance analysis at an affordable computational cost without any simulation enhacement. We discussed the importance of considering the vulnerability of the infrastructure as a whole, by showing that the effect of adding an additional link in the topology is not always significative. For experiments considering weaker infrastructure, adding a link generates neglectable reliability improvements; while for stronger infrastructure, adding the same link results in a significant reliability improvement.
Although, several extensions and improvements can be considered for the proposed method, we remark that data availability is the main challenge for implementing this technique in real-life networks. In this analysis, we used information gathered and validated by studies in the area where the REUNA network components lie. We found several detailed studies for urban areas. However, the lack of information in rural areas forces us to make some assumptions that do not change the complexity of the study. For instance, at present, there is no detailed information on some parameters; as incredible as it might be, there are areas where not even the exact location of the buried optical fiber is known. Governments have realized the strategic value of having it, and nowadays, efforts are being made to have a more accurate map of the fiber optic network (see the 2018 study in the Chilean case [48]). Our work contributes to making visible the advantages of having this information. Also, a more detailed seismic model and the SA between site correlations should be incorporated. Additional aspects of inter-infrastructure dependence (for instance between the data network and the electrical network) should be addressed since for seismic hazards these dependencies have a significant impact on reliability.

APPENDIX A GMPE EQUATIONS
For completeness' sake, we reproduce the GMPE equations used in our study. Following the diagram described in Figure 3 the SA(T ) is calculated using the equations proposed in [32], these are adequate for this research as they are calibrated for the Chilean geographic characteristics. The SA(T ) is calculated according to the following equation, (the SA is equivalent to the PGA at T = 0): ln(SA(T )) = θ 1 + f source + f path + f event/depth where M w represents the seismological scale of moment magnitude, Z h is the hypocentral depth measured in kilometers, R is the distance between the hypocenter and the considered component, PGA M for locations with V S30 equal to 1000 m/s whose parameters have been mentioned in Table 1, values of coefficients C 1 , θ 9 , C 4 , V lin , b, c and n are values adopted directly from the model in [32]. T corresponds to the spectral periods, and takes values between 0.02 and 10 s; the values of the θ parameters depend on T . We use T = 0 for the edge acceleration since we consider that they represent buried fiber optic cables; for the nodes we use T = 0.25 associated to the type of infrastructure used.

APPENDIX B RELIABILITY FOR OTHER CRITICAL ACCELERATION VALUES
The critical acceleration describes how much horizontal acceleration is necessary for landslide displacement to occur with a given probability. As we have indicated no detailed study of liquefaction in Chile is available, so we assume that the entire area under consideration is affected. In this appendix we show how the reliability changes with a critical acceleration of k c = 0.4. Table 13 summarizes the marginal failure probabilities for each component and the correlations by pairs of links or pairs of nodes obtained for all the network's components. The results show a significant failure correlation.     obtained. Lower critical acceleration values translate into more vulnerable seismic zones. Table 16 shows no statistically significant improvement between the R-3F and R-3F+1 for a critical acceleration of k c = 0.4, Table 17 shows a significant statistical improvement between R-3F and R-3F+1 with respect to the original network for pairs of nodes over long distances, for a critical acceleration of k c = 0.4.
In Fig. 9 and Fig. 10 again when we disaggregate with the distance, the impact of the proposed model can be quantified. Reliability is significantly improved for those users that are more distant. HÉCTOR CANCELA (Senior Member, IEEE) received the degree in computer systems engineering from Universidad de la República, Uruguay, in 1990, and the Ph.D. degree in computer science from the University of Rennes 1, INRIA, Rennes, France, in 1996. He was the Dean of the Engineering School, Universidad de la República, from 2010 to 2015. He is currently the Head of the Computer Science Institute, Engineering School, Universidad de la República. He is also a Researcher with the National Program for the Development of Basic Sciences (PEDECIBA), Uruguay. His research interests include operations research techniques, especially in network models and stochastic models, applied jointly with optimization methods for solving problems in different areas. He was the President of Centro Latinoamericano de Estudios en Informática (CLEI) and a member of the Uruguayan Engineering Academy (ANIU). He was a former President of Asociación Latino Ibero Americana de Investigación Operativa (ALIO). He has participated in the IEEE CS Education Award Selection Committee.