Integrated seismic risk and resilience assessment of roadway networks in earthquake prone areas

Abstract

Intercity networks constitute a highly important civil infrastructure in developed countries, as they contribute to the prosperity and development of the connected communities. This was evident after recent strong earthquakes that caused extensive structural damage to key transportation components, such as bridges, overpasses, tunnels and geotechnical works, that in turn led to a significant additional loss associated with the prolonged traffic disruption. In cases of seismic events in developed societies with complex and coupled intercity transportation systems, the interdependency between citizens’ life and road functionality has further amplified the seismically-induced loss. Quantifying therefore, the resilience of road networks, defined as their ability to withstand, adapt to, and rapidly recover after a disruptive event, is a challenging issue of paramount importance towards holistic disaster risk mitigation and management. This study takes into account the above aspects of network resilience to earthquake loading and establishes a comprehensive, multi-criterion framework for mitigating the overall loss expected to be experienced by the community due to future earthquake events. The latter is decoupled into the direct structural damage-related loss and the indirect loss associated with the travel delays of the network users, as well as the wider socio-economic consequences in the affected area. In order to reflect the multi-dimensional nature of loss, a set of novel, time-variant indicators is herein introduced, while cumulative indicators are proposed for assessing the total loss incurred throughout the entire recovery period. This probabilistic risk management framework is implemented into a software to facilitate informed decisions of the stakeholders, both before and after a major earthquake event, thus prioritizing the pre-disruption strengthening schemes and accelerating the inspection and recovery measures, respectively.

1 Introduction

Intercity transportation networks constitute a vital component of prosperity in modern, dense populated societies by facilitating the mobility of people, goods and services. Their smooth and undisruptive operation is crucial for ensuring sustainability after extreme natural disasters. Earthquakes for instance, have caused extensive damage worldwide primarily to seismically sub-standard road or highway components (i.e., bridges, overpasses, tunnels etc.) (Kawashima and Buckle 2013; Shen et al. 2014), as well as to adjacent geotechnical works that influence road functionality (e.g. slopes). These damages have led to enormous direct and indirect loss to the affected areas (Zhou et al. 2010). Direct loss is related to the repair of the damaged components, if one for the sake of quantification, neglects the priceless loss of human life, while indirect loss refers to the reduced functionality of the road network and the subsequent increase of travel time, the disturbance to social and professional life, business interruption, additional transportation cost and environmental implications (Kiremidjian et al. 2007; Dorra et al. 2013).

Direct and indirect loss associated with future seismic events affecting highway networks and their secondary roadways is assessed probabilistically (Chang et al. 2000; Werner 2000; Günneç and Sibel Salman 2011; Dong et al. 2014a) by a seismic risk assessment that convolutes hazard at the site(s) of interest (Bommer and Crowley 2006; Sokolov and Wenzel 2011; Han and Davidson 2012) with the structural/geotechnical vulnerability (DesRoches et al. 2008; Kwon and Elnashai 2010; Shafieezadeh et al. 2014), as well as the altered traffic flow (Zhou et al. 2014) and the wider economic, social and environmental consequences of both physical damages and traffic diversion. Primary purpose of risk assessment is to support decision making related to pre-event loss mitigation strategies, such as, seismic retrofit and emergency planning (Fragiadakis et al. 2013; Tesfamariam and Goda 2013). Key concept to the reliable risk assessment and to the effective decision-making is network resilience that encompasses the dimensions of network capacity, redundancy and recovery time to express its ability to withstand and adapt to a natural disaster, while being able to recover and restore the services offered quickly. Note that several alternative definitions exist as per different systems and contexts as reviewed by Hosseini et al. (2016). Due to the complexity of roadway networks and their extension in large geographical areas as well as the multidimensional nature of network resilience concept itself (Cavalieri et al. 2012), incorporating resilience metrics to risk-informed decision-making procedures faces several challenges that hinder the practical application of such innovative approaches.

One major difficulty of doing so is that, first, damage in one network component (i.e. bridge, tunnel etc.) influences the functionality of the entire network, which subsequently may impede the emergency response, the recovery activities, the rehabilitation process and the accessibility of critical transportation (ports, airports, train) or operation of energy facilities. The interdependence between modern life activities and the difficulty in assessing loss pertaining different sectors (e.g. economy, society, and environment) adds further to the complexity to the indirect loss estimation (Lounis and McAllister 2016; Tapia and Padgett 2016) a fact that was highlighted by some studies that identified multiple dimensions of resilience (Bruneau et al. 2003; Cimellaro et al. 2016) but has not yet been addressed in a reliable, computationally efficient and quantitative manner. Moreover, a roadway consists of a number of components that, after an earthquake, experience different levels of damage. As the spatially variable damage is restored gradually but not simultaneously for the entire network, its spatial distribution also varies in time and subsequently, the post-earthquake traffic and the associated indirect loss are ultimately recovery-dependent. This time dimension of resilience in conjunction with the traffic assignment that is necessary for simulating traffic conditions, further imply that the Origin–Destination (OD) matrix that controls the traffic generation is dynamic, in other words, the drivers’ behavior is event-dependent (i.e., they tend to alter their destination after a disaster and divert towards home, shelter or emergency facilities). For these reasons, only a wider network perspective that includes the refined network simulation and detailed traffic flow estimation for consecutive recovery phases can capture the time variation of the indirect loss and provide a reliable approach for its overall estimate, throughout the entire recovery period.

Time has an additional effect on infrastructure resilience as ageing of transportation infrastructure (Alipour and Shafei 2016) affects not only the initial state of functionality of the pre-disaster roadway network (as seen in Fig. 1) but effectively determines the immediate resilience drop of functionality to the residual (robustness) level $Q\left( {t_{0E1} } \right)$ right after the natural disaster, a fact that is commonly overlooked.

Fig. 1

Schematic representation of disaster resilience. (Reproduced with permission from Cimellaro et al. 2010)

Another difficulty is that the roadway risk associated with seismic hazard cannot be evaluated using conventional approaches. The latter is because, probabilistic seismic hazard analysis (PSHA), typically used for single structure analysis is site-specific, and hence not applicable for the case of a road network spanning over a wide region, where a number of geographically distributed network components (i.e. bridges, tunnels etc.) are exposed to different levels of seismic hazard. Most importantly, the spatial distribution of seismic intensity (i.e., peak ground acceleration, spectral acceleration among others) after an earthquake event attenuates from the source to the sites of interest, resulting to damage distribution within the network, which varies in space. The problem is usually tackled on the basis of deterministic specific scenario events (Padgett et al. 2010), or, if a stochastic analysis is foreseen, by using a set of probabilistic earthquake scenarios. While the former approach is more tractable, the latter can better account for the prevailing uncertainties and lead to an unbiased risk assessment. More precisely, Werner et al. (2000) developed a probabilistic risk-based methodology that incorporated a walk-through analysis of roadway networks by simulating earthquakes over a certain period of time (e.g. 1000 years). In some other studies conventional Monte Carlo simulations were applied to integrate the uncertainties associated with seismic hazard (Chang et al. 2012). The computational demand of a risk assessment that includes probabilistically modeled earthquakes deriving from the application of the above approaches is considerably high, and it can be hardly addressed, especially in cases where network resilience is taken into account and refined network analyses are employed for multiple phases of network operation. Recognizing, the high computational cost of the probabilistic management of seismic hazard at a regional basis, as opposed to conventional site-specific hazard, several researchers have proposed methods for generating a limited number of hazard-consistent seismic maps (Chang et al. 2000; Shinozuka et al. 2003). These approaches are based on the adjustment of the probabilities of a number of earthquake scenarios so that the hazard associated to the whole range of possible earthquake events in the region is approximated. However, this approximation is achieved in terms of ground motion intensity values on a few arbitrary selected network sites and thus the whole range of possible intensity measure distributions is not taken into account.

In light of these challenges, a number of studies proposed frameworks for assessing seismic risk of roadway networks. In the most comprehensive of these frameworks, risk is expressed in the form of risk curves which give the annual probability of exceeding various loss levels in terms of travel delay (Shiraki et al. 2007), monetary units (Werner et al. 2001), accessibility (Miller and Baker 2016), etc. These efforts, however, do not explicitly deal with network resilience aspects since traffic rerouting and the consequent indirect costs are only approximately estimated (Shiraki et al. 2007) or the evolution of the network functionality throughout the recovery period is neglected (Du and Peeta 2014; Miller and Baker 2016). On the other hand, much of existing literature in the field of resilience assessment is rather conceptual (Zobel and Khansa 2014; Cimellaro et al. 2016). In fact, the few existing approaches for quantifying seismic resilience of infrastructure networks provide formulas for calculating resilience metrics (Franchin and Cavalieri 2015) but do not integrate these formulas into a comprehensive risk framework. For instance, Alipour and Shafei (2016) presented a framework for the analysis of the resilience of highway networks with deteriorating components due to aging but therein seismic hazard is treated deterministically. Similarly, Ouyang et al. (2012) introduced a three-stage resilience analysis framework for urban infrastructure systems that, among others, deals with the uncertainty in the seismic hazard, however, without including procedures for estimating system performance (i.e. the performance metrics needed for assessing resilience metrics are assumed to be known).

Based on the above, there is a need for a framework that treats uncertainties related to hazard, damage and consequences in a balanced way, while at the same time explicitly accounts for time-dimension of recovery, network functionality and multidimensionality of seismic losses, three key aspects for quantifying resilience.

Along these lines, this paper seeks to develop a robust, multi-dimensional, applicable, framework for the quantification of road network risk and resilience in earthquake regions. More precisely, it aims to:

• (a) develop a probabilistic, multi-event, seismic hazard approach, specifically tailored to network resilience analysis. The proposed approach is similar to the standard PSHA but in this case hazard is computed independently for every possible seismic source without aggregating the multiple seismic sources to account, in a more refined manner, for the impact of seismic source geometry and orientation with respect to the roadway network topology. To this end, the active tectonic faults and seismic zones are discretized into segments and seismic maps with specific annual exceedance probability are generated. This serves the purposes of event-driven, pre- and post-quake traffic and consequence analysis.

• (b) include multiple phase, traffic scenarios, resulting from a Monte Carlo simulation of post-earthquake damage samples. Each damage sample is associated to several recovery phases following the gradual restoration of the damaged components, thus introducing the time dimension of roadway network resilience. The above disaggregation of the recovery time into phases is utilized for the time variant analysis of the post-earthquake traffic.

• (c) introduce time-variant and cumulative roadway indicators to quantify resilience and seismic risk respectively through indicators that capture the multidimensional nature of resilience at a community level.

With the aid of an open GIS-based software developed herein and an interconnected dynamic traffic assignment engine (Zhou et al. 2014), the framework permits roadway stakeholders to form alternative risk management strategies to minimize loss and optimize functionality after an earthquake disaster.

2 Overview of the proposed resilience framework

The framework proposed herein (Fig. 2) initiates with the Network definition of its topology and key components (i.e., bridges, tunnels, slopes etc.), as well as network operation and pre-disaster traffic data. Then, regional seismic hazard is analyzed through an extension of the conventional PSHA. As it was mentioned, this approach involves the discretization of the nearby active seismic sources into segments, the definition of seismic hazard levels with specific annual exceedance probabilities and ultimately the generation of a number of probabilistic seismic maps. Subsequently, key components are classified into fragility classes of similar structures that can be classified as having identical probability of failure expressed in terms of fragility curves. Seismic maps and key component fragilities are combined to formulate probabilistic damage distributions. The latter are then used for Monte Carlo simulations that generate deterministic damage samples. The damage samples together with a restoration matrix that describes the estimated repair times of every damaged key component are converted to multiphase traffic scenarios that cover the entire recovery period. Traffic cost deriving from the traffic assignment performed for the pre-earthquake network state, as well as for all the post-quake phases, is coupled with structural cost associated to the probabilistic damage distributions. Several cumulative and time-variant indicators are defined to assess the impact of probable earthquake events with certain annual rates of exceedance on the resilience of the network. These indicators can then be used by the stakeholders to assess the resilience of the “as built” network and to identify candidates of possible (pre-earthquake), effective risk management strategies that are expected to significantly improve the resilience of the roadway network in case of a future earthquake event.

Fig. 2

General workflow of the proposed framework

3 Network definition

3.1 Network topology

Every network intersection and every network location from which the drivers are originated or are heading to, is considered a node, while two successive nodes are connected with a link. Α unique ID number and the corresponding coordinates are assigned to every node. Road links are defined by their length and the IDs of their edge nodes. For visualization purposes, this is reflected on a GIS platform; in what concerns the software developed herein, shape files are loaded on a Matlab-based script as discussed in more detail in Sect. 9.

3.2 Pre-earthquake traffic conditions

Traffic flow estimation process requires a reference (initial, pre-earthquake) origin–destination (OD) matrix as an input. This matrix describes the traffic demands for all the possible travel origins and destinations of the network and it is typically provided by the stakeholder. Because of the relatively stable nature of the pre-earthquake traffic conditions, the drivers’ behavior is rather consistent and hence, it can be easily expressed in the form of a static OD matrix format, which is then assumed constant over the entire pre-earthquake period. Based on that matrix and also on the additional input of the traffic capacity and speed limit of every network link, traffic flows and times over the whole network are calculated according to Zhou et al. (2014). This involves the assumption that drivers will always choose travel paths that minimize their travel times, given the links capacity and speed limit, which is a rather realistic assumption at least under normal (everyday) traffic patterns.

3.3 Key components

All the bridges, overpasses, tunnels and slopes affecting the network (Fig. 3), herein referred as key components (KC), are considered to be potential sources of seismic loss. Damage to other network components (e.g. straight road segments, pavements, small technical works etc.) is neglected.

Fig. 3

Examples or traffic disruption due to key network component failures (bridges, tunnels and slopes)

The I key components of the network are firstly identified and then defined through the description of their IDs and coordinates. Site-specific soil type according to Eurocode 8 (EC8 2004), total rebuilding cost (TBC), repair cost ratios (RCR) for every possible damage state (i.e. damage state-dependent ratio of the repair cost to the total rebuilding cost) are also provided for each key component at this stage. Moreover, the IDs of the network link affected functionally by the damages of each key component has to be defined.

3.4 Network operation

Specific network points of interest (POIs) that need to remain accessible after an earthquake event are defined. The same applies to links whose failure or excessive traffic may have an undesirable environmental impact. The network region is then divided into activity zones with similar organizational structure (e.g. industrial zones, agricultural areas, etc.) and each POI and network link is assigned to one zone. All the above data is utilized in the resilience analysis as discussed in Sect. 7.

4 Seismic hazard analysis

The spatial variation of the intensity measure (IM) of interest along the road network is derived from a seismic hazard analysis. The conventional probabilistic seismic hazard assessment (PSHA) is based on a method introduced by Cornell (1968) that integrates the temporal and spatial distribution of earthquake events of different magnitudes with a Ground Motion Prediction Equation (GMPE). The integration is accomplished on the basis of the Total Probability Theorem as (Baker 2013):

$$\lambda (IM > x) = \sum\limits_{i = 1}^{{n_{sources} }} {\lambda (M_{i} > m_{\rm min} )} \cdot \int\limits_{{m_{\rm min} }}^{{m_{\rm max} }} {\int\limits_{0}^{{r_{\rm max} }} {} P(IM > x/m,} r) \cdot f_{{M_{i} }} (m) \cdot f_{{R_{i} }} (r) \cdot dr \cdot dm$$

(1)

where $\lambda (IM > x)$ is the annual rate λ of exceeding an IM value x, n_sources is the number of sources considered, $P(IM > x/m,r)$ is the probability that the IM of interest will exceed a value x given the earthquake magnitude m and the source-to-site distance r; $\lambda (M_{i} > m_{\hbox{min} } )$, is the mean annual rate of occurrence of an event with magnitude greater than m_min for a source i; $f_{{M_{i} }} (m)$ and $f_{{R_{i} }} (r)$ are the probability density functions for the magnitude and the source-to-site distance of source i, respectively.

Notably, conventional seismic hazard maps that combine IMs to several sites within a region of interest, estimated on the basis of the above integration, are not compatible with the problem of post-earthquake traffic distribution studied herein. This is because post-event traffic depends on the individual post-quake operation of each network key component (bridge, tunnel etc.). The latter depends on the specific earthquake rupture scenario (i.e., triggering event) examined, that is, a single earthquake event with specific magnitude and rupture location (Bommer and Crowley 2006; Sokolov et al. 2010). As discussed in the introduction, to tackle this problem several methods employ Monte Carlo simulations of earthquake occurrences that lead to an excessive number of seismic scenarios, while some other generate a limited set of hazard-consistent scenarios but do not explicitly treat the uncertainty in the rupture location. To address this issue, an efficient treatment of the uncertainty in seismic hazard is proposed, which is based on the discretization of the active seismic sources within the region of interest into segments. Moreover, specific hazard levels (HL) corresponding to predefined annual rates of exceedance that are needed for the resilience-based network analysis and decision-making are defined by the stakeholder (e.g. HL1, HL2, HL3 and HL4 that have 2%, 0.20%, 0.10% and 0.05% annual rate of exceedance, respectively). Assuming that each segment corresponds to a possible earthquake rupture, seismic maps associated with each hazard level are generated for each source segment. The generation of the seismic maps is based on a modification of Eq. (1) introduced herein to account for the source discretization. Hence, the uncertainty in the rupture locations is explicitly captured through the generation of seismic maps that correspond to all the source segments, while the uncertainties in the earthquake rates and magnitudes are reflected by the hazard levels and the relevant annual rate of exceedance associated to each map.

The application of the above method involves the following steps:

• 1. U seismic sources (e.g. tectonic faults, seismic zones) capable of producing damage to the network under consideration are first identified. Similarly to the conventional PSHA, a seismic source is represented by either a line or an area of uniformly distributed seismicity (i.e., an earthquake rupture initiation at any location within a source is equiprobable).

• 2. The annual rate of occurrence of earthquakes greater than m_min for each identified source u, $\lambda (M_{u} > m_{\hbox{min} } )$, is defined. It is noted that $\lambda (M_{u} > m_{\hbox{min} } )$ corresponds to $\lambda (M_{i} > m_{\hbox{min} } )$ parameter needed for a conventional PSHA and can be extracted by existing seismotectonic studies of the region of interest.

• 3. The U seismic sources are divided into segments of equal seismicity. For a source u a total of $N_{u}^{{}}$ segments with equal length or area are defined:

  $$N_{u} = \frac{{\lambda (M_{u} > m_{\hbox{min} } )}}{{\lambda_{ref} }}$$

  (2)

  where $\lambda_{ref}$ is a reference value that is equal to the annual rate of occurrence of earthquakes greater than m_min to any of the $L = \sum\limits_{u}^{U} {N_{u} }$ discrete segments of the U seismic sources.

• 4. The segment-to-site distances $R_{l,i}^{{}}$ between all the pairs of source segments derived after the dicretization of the U sources and network key components are defined (Fig. 4).

  Fig. 4

  

  Epicenter to site distances associated to the segment l

• 5. The annual rate of exceeding an IM value at the location of a key component i, due to the seismicity of segment l ∈ {1,2,…,L} is given by Eq. (3), which is essentially an adjustment of Eq. (1) to account for source discretization:

  $$\lambda_{l,i}^{{}} (IM > x) = \lambda_{ref} \cdot \int\limits_{{m_{\hbox{min} } }}^{{m_{\hbox{max} } }} {P(IM > x/m_{{}} ,} R_{{l,i_{{}} }} ) \cdot f_{{M_{u} }} (m) \cdot dm$$

  (3)

  where $f_{{M_{u} }} (m)$ is the probability density function of earthquake magnitudes related to seismic source u, that is, the source to which segment l belongs to (the estimation of this function is also a prerequisite in a conventional PSHA and is based on the assumption of an earthquake recurrence model such as Gutenberg-Richter law).

  P(IM > x/m, R_l,i) is the probability of exceeding a certain IM value given the earthquake magnitude and the segment-to-site distance R_l,i which is estimated on the basis of a ground motion prediction equation.

• 6. Κ number of hazard levels (HL) are defined, each one associated with an annual rate of exceedance λ_k (e.g. HL1, HL2, HL3 and HL4 that have 0.02, 0.002, 0.001 and 0.0005 annual rate of exceedance, respectively). Then, for a given hazard level k and seismic segment l, Eq. (3) is employed and intensity measure values to all the key component sites are derived as shown in Fig. 5.

  Fig. 5

  

  Indicative derivation of the (unamplified) PGA values corresponding to a hazard level k and a seismic segment l

• 7. All the intensity measure values derived from the previous step are multiplied with an amplification factor to account for the effects of local soil and site conditions on the ground motion. Given the soil class at each the site of eack network key component (i.e., bridges, tunnels etc.), this factor is approximately taken equal to the Eurocode 8 (EC8 2004) soil factor (S) that is typically used for adjusting the elastic response spectrum.

• 8. Given a hazard level k, L ground motion maps, with the same annual rate of exceedance λ_k, are generated out the L possible earthquake ruptures assumed after the segment discetization of U faults identified. Thus, in total K × L maps are generated.

Note that several source-to-site distance definitions were introduced by the different GMPEs proposed in the literature depending on the assumed rupture location (e.g. epicenter, hypocenter, closest point on the rupture surface etc.). Segment-to-site distances $R_{l,i}^{{}}$ needed for the application of the proposed hazard analysis method must be defined in accordance to the source-to-site distance definition (i.e. by treating each segment as a seismic source) of the GMPE utilized for estimating $P(IM > x/m,R_{l,i} )$ (Eq. 3).

It is acknowledged that intra-event (between sites) and inter-event (between earthquake events) IM correlation is not explicitly taken into account in the proposed hazard approach, as mentioned for reasons of computational efficiency. This issue has been presented in detail elsewhere (Jayaram and Baker 2010) but it is deemed beyond the scope of this study.

5 Fragility analysis

5.1 Probabilistic damage distribution

Fragility curves reflect the probability of a structure to experience damage that exceeds predefined Damage States (DS) given a seismic intensity level, expressed by means of an appropriate intensity measure (IM), typically peak ground acceleration (PGA) or spectral acceleration at the fundamental period (S_a(T₁)), among others. Fragility curves are usually plotted assuming a lognormal distribution function as follows:

$$P_{{S > DS_{t} /IM = im}} = \varPhi \left[ {\frac{1}{{\beta_{t} }}\ln \left( {\frac{im}{{im_{mt} }}} \right)} \right]$$

(4)

where Φ is the standard normal cumulative distribution function. im_mt is the median threshold value of the IM associated with damage state t. β_t the lognormal standard deviation of the IM associated with damage state t

A set of four Damage States (DS1 to DS4) is assumed to define the fragility of all key component classes, corresponding to minor, moderate, extensive damage and collapse. Given the network component (bridge, tunnel, slope etc.) classification, the appropriate fragility curves are assigned to every single key component of the network. Because the value of the Intensity Measure at the location of each key component is different, the probability that the component will experience damage corresponding to each Damage State (DS1-4, Fig. 6) is also different and can be directly computed from the limit-state exceedance probabilities:

$$\begin{aligned} P_{{DS_{0} /IM}} & = 1 - P_{{S \ge DS_{1} /IM}} ,\;P_{{DS_{1} /IM}} = P_{{S \ge DS_{1} /IM}} - P_{{S \ge DS_{2} /IM}} \\ P_{{DS_{2} /IM}} & = P_{{S \ge DS_{2} /IM}} - P_{{S \ge DS_{3} /IM}} ,\;P_{{DS_{3} /IM}} = P_{{S \ge DS_{3} /IM}} - P_{{S \ge DS_{4} /IM}} ,\;P_{{DS_{4} /IM}} = P_{{S \ge DS_{4} /IM}} \\ \end{aligned}$$

(5)

Fig. 6

Sample set of fragility curves corresponding to the four damage states (minor, moderate, extensive and complete damages)

Wherever possible, classes of key components with similar structural, geometrical or geotechnical characteristics (hence similar vulnerability) are used. It is noted that the potential decision to group different network components into classes of similar fragility (for instance, a number of overpasses with the same structural system and approximately equal dimensions and sections) is computationally very practical but can introduce additional uncertainty to the problem and has to be made with caution. However, it can be a reasonably good simplification for the case of extended networks where ad-hoc, structure-specific fragility analysis (Stefanidou and Kappos 2017) is not possible. Overall, it is a good strategy to group components of minor importance in classes and perform more detailed fragility assessment for potentially critical network components, such as those associated with high traffic loads or increased travel times in case of failure. Refinement of fragility could be further increased if the selection of earthquakes needed to analyze the structures’ response takes into account the discretization of the seismic sources proposed in the previous section.

After the key components are organized into classes of equal fragility, the probabilities of Eq. (5) are defined for each one of the K x L seismic maps derived from the seismic hazard analysis and each one of the I key components of the network. The above 5 × I probabilities associated to a seismic map k, l (i.e. the seismic map that corresponds to hazard level k and seismic segment l) express the uncertainty in the damage estimation given the hazard and are cumulatively referred as the probabilistic damage distribution k, l.

5.2 Seismic fragility of structural components

The structural components of a road network that are usually damaged after an earthquake are bridges and overpasses. Several alternative approaches have been proposed for the estimation of their fragility, mainly being either empirical, analytical or numerical. Empirical fragility curves are derived by statistically processing of actual damage data collected after major earthquake events (e.g., Northridge and Kobe earthquakes (Basoz et al. 1999; Elnashai et al. 2004). However, due to the uncertainty associated with visual inspection and the inherent subjectivity in defining damage states on site, analytical methods have also been developed based on non-linear numerical models (Mackie and Stojadinovic 2004; Nielson and DesRoches 2007). In the methodology presented herein, fragility curves of the structural components of the network can be either derived numerically, or existing (empirical or analytical) expressions can be retrieved for specific classes directly from the literature (Choi and Desroches 2004; Padgett and DesRoches 2007).

5.3 Seismic fragility of geotechnical components

In contrast to the extensive research conducted on the seismic fragility of structural road network components, fragility assessment methods for geotechnical components is very limited. To the best of the authors knowledge, up to now there is only one procedure for deriving fragility curves for shallow (Debiasi et al. 2012; Argyroudis and Pitilakis 2012) or underground (Andreotti and Lai 2015) tunnels that actually cannot be directly applied to deep tunnels (e.g. mountain tunnels). Moreover, even though landslides are often triggered by strong earthquakes, the seismic fragility of slopes affecting roads has not been sufficiently studied yet (Liu et al. 2017). Regardless the lack of previous research in the area, it is recommended that the vulnerability of geotechnical components is taken into consideration even if the parameter values required in Eqs. (4–5) are based on engineering judgment rather than on refined fragility analyses. Note that as the Intensity Measures typically employed for the fragility analysis of tunnels (i.e., PGD, PGV) are generally different than those used for bridges (PGA, S_a(T₁)), an IM transformation is required in order to ensure compatibility between the fragility estimates of the structural and geotechnical components for the same seismic map.

6 Recovery analysis

6.1 Generation of recovery phases

All the probabilistic damage distributions derived within Fragility analysis are simulated employing Monte Carlo method. Assuming the size of sampling S, a total number of S damage samples is formed, including schemes of damage states (DS0, DS1, DS2, DS3, DS4, DS5) for all the network key components and each one of the KxL probabilistic damage distributions.

Having identified samples of spatial distribution of damage within the network, the time associated with the gradual recovery of the network is defined. This is made feasible by defining a restoration matrix of dimensions I×4, which quantifies the time required to repair different damage levels (DS1 to DS4) of each different key component i.

As explained in Sect. 3.3, each key component (bridge, tunnel etc.) affects the functionality of a certain roadway network link. Reversely, each network link j may be affected by zero, one or several key components. As a result, the above damage and repair time of key components is effectively reflected on the recovery of the links as well. Therefore, depending on the sample examined, the damage state identified for each key component and the restoration matrix, a set of repair time values $\left\{ {rt_{j1} ,rt_{j2} , \ldots ,rt_{jYj} } \right\}$ is assigned for each network link j. The number Y_j of the repair times associated is equal to the number of the key components affecting the functionality of the respective link. The maximum repair time $ct_{j} = \hbox{max} \left\{ {rt_{j1} ,rt_{j2} , \ldots ,rt_{jYj} } \right\}$ of the key components within a link will define the closure time of the link itself, assuming that link j cannot be used by any vehicle until all the key components that affect it are fully restored.

After identifying closure times for all the J network links $\left\{ {ct_{1} ,ct_{2} , \ldots ,ct_{J} } \right\}$ the durations are sorted in ascending order, defining distinct points in time $\left\{ {ct_{1} ,ct_{2} , \ldots ,ct_{J} } \right\}$ when each one of the links opens. Because there is the chance that different links recover (open) at the same time, the characteristic opening time vector has an equal or lower number of terms $\left\{ {ct_{{_{1} }}^{\prime} ,ct_{{_{2} }}^{\prime} , \ldots ,ct_{{_{C} }}^{\prime} } \right\}$ where C < J. Taking also into consideration that a number of network links may be non-associated with the damage of any key component, it is also possible that some of the first terms are zero. By eliminating the zero recovery time terms, the number P of the post-earthquake recovery phases is defined and the characteristic opening times vector has p terms $\left\{ {ct_{{_{1} }}^{\prime\prime} ,ct_{{_{2} }}^{\prime\prime} , \ldots ,ct_{{_{P} }}^{\prime\prime} } \right\}$.

Based on the above, the duration t_p of each phase p (i.e., until the opening of the next closed link) is obtained from:

$$t_{p} = \left\{ {\begin{array}{ll} {ct_{1}^{\prime\prime} } & {{\text{if}}\;{\text{p = 1}}} \\ {ct_{p}^{\prime\prime} - ct_{p - 1}^{\prime\prime} } & {{\text{if}}\;{\text{p > 1}}} \\ \end{array} } \right.$$

(6)

where $ct_{1}^{\prime\prime}$ and $ct_{{_{2} }}^{\prime\prime} , \ldots ,ct_{{_{P} }}^{\prime\prime}$ are the lower non-zero link closure time and the unique link closure times sorted in ascending order, respectively, expressed in the same units (i.e., days, weeks, months) that the restoration matrix was initially defined. The functionality of the network is then described by a binary network state vector ${NS}_{p}$ taking values of (1) for each network link that is open and (0) when closed.

The above procedure is repeated for the defined seismic hazard levels, leading to KxLxS damage samples as shown in Fig. 7 for which the number of post-earthquake phases is to be found. This is because the latter depends on the closure times of the network links and varies from damage sample to another. For example, a higher hazard level is expected to generate more extensive in space and significant damage and hence, naturally, a higher number of recovery phases. Overall, the total number of the recovery phases to be analyzed is equal to $P_{tot} = \sum\nolimits_{k}^{K} {\sum\nolimits_{l}^{L} {\sum\nolimits_{s}^{S} {P_{k,l,s} } } }$, where P_k,l,s is the number of phases derived from the recovery analysis of each damage sample and each seismic map (for K hazard levels of interest and L discretized segments of all sources).

Fig. 7

Development of the recovery phases associated to the K × L × S damage samples

6.2 Traffic assignment

With the aid of a traffic assignment model (Zhou et al. 2014) traffic flows and travel times are estimated for the pre-earthquake traffic scenario and the P_tot post-earthquake traffic scenarios that correspond to the P_tot phases of the K × L × S damage samples. Pre-earthquake traffic scenario is defined by the Network topology and the Pre-earthquake traffic data (Sect. 3) while the traffic scenarios corresponding to post-earthquake phases are derived by setting to zero the traffic capacity of the links that are closed to traffic according to the corresponding network state vector ${NS}_{p}$.

At this point it is noted that travel demand after an earthquake event can be significantly altered due to the emergency situation impact on the drivers’ choices. For instance, considering the increased travel time and cost some drivers may chose to postpone or cancel their trip, head to emergency facilities, change destination or simply return home. Moreover, as everyday activities rebound gradually and life comes back to normality, this altered post-earthquake demand may also evolve in time. In the proposed framework, this dynamic aspect of the post-earthquake traffic demand can be taken into account by assigning a new OD matrix to every recovery phase. As the prediction of this OD transformation involves social and behavioral aspects that introduce a level of additional uncertainly, it is not explicitly addressed herein, however, the methodology and software presented in Sect. 9 permits such matrix update in case of existing data.

Another interesting point is that unlike to the pre-earthquake traffic conditions where drivers usually know which is the fastest route for their trip and tend to select it, this is not always the case after an earthquake event. However, in this study it is assumed that emergency management authorities and (real time traffic) navigation tools can aid travelers in choosing the fastest route for a certain trip and hence the fastest-path algorithm is again applied for calculating traffic flows over the network for every post-earthquake phase. The above assumption holds true only if communication networks (e.g. 3G/4G internet services) are resilient enough to ensure the smooth traffic data transmission after an earthquake event.

7 Resilience assessment

7.1 Quantification of resilience

To assess the resilience of the road network and the earthquake impact on the affected area, a set of Time-variant and Cumulative (total consequence throughout the recovery period) indicators are proposed, their main difference being that the former reflect the evolution of wider consequences over the recovery phases attributed to a damage sample, while the latter quantify the cumulative structural and traffic cost per seismic map. It is deemed that both (time-variant and cumulative) proxies are useful for informed decision-making of the stakeholders as per the necessary measures needed to be taken before an earthquake to minimize post-quake loss and improve network resilience.

7.2 Time-variant indicators

An earthquake event results in an abrupt disorder to the network functionality that is gradually restored following the repair of the network damage. Time-variant indicators, focus on the description of the time evolution of losses incurred by the community from the onset of the earthquake throughout the recovery period, and are introduced herein for the resilience-based assessment of the network ability to survive, adapt and recover after an earthquake. In order to explicitly reflect the time-dimension aspect, these indicators are distinctively evaluated for every damage sample arising from MC simulation.

7.2.1 Network functionality

A Network functionality-time relationship is an indicator that reflects the network functionality as a function of time. In this relationship, network functionality is given as a percentage of full network operation, which is the reference functionality before the earthquake, weighted by the importance of the links that remain operative after the earthquake. It is assumed that a link is non-functional either when it is closed (due to damage) or if it is intact but does not serve to any travels.

The importance factor (γ_j) of every network link is defined as the ratio between the link initial traffic load and the sum of initial (pre-quake) traffic loads of all the network links:

$$\gamma_{j} = \frac{{V_{j0} }}{{\sum {\text{V}}_{\text{j0}} }}$$

(7)

Network functionality$f_{p}$ is then given by the sum of the importance factors of the links that are functional:

$$f_{p} = 100 \cdot \sum\limits_{j = 1}^{J} {f_{pj} }$$

(8)

where

$$f_{pj} = \left\{ {\begin{array}{ll} {\gamma_{j} } & {if\;j\;link\;is\;functional\;during\;phase\;p} \\ 0 & {if\;j\;link\;is\;non{ - }functional\;during\;phase\;p} \\ \end{array} } \right.$$

Network functionality is calculated for every phase of the recovery period (Fig. 8). Every horizontal branch of this graph stands for one phase of the recovery period. At the end of the last recovery phase, network functionality is restored back to 100%.

Fig. 8

“Functionality-time” correlated to a sample damage sample

It is noted that the calculation of this indicator does not include quantities that are explicitly related to traffic such as flows and travel times. However, f_p values are useful for a preliminary (i.e. before the execution of the traffic assignment) assessment of the possible earthquake consequences (e.g. they may be used for modifying the pre-quake OD matrix to generate post-quake OD matrices which in turn may be used in the subsequent traffic assignment of the recovery phases).

7.2.2 Additional traffic cost

By developing the complete set of P post-earthquake traffic phases associated to a damage sample s, the increase in travel time is derived along with the associated additional (travel) cost. The latter is expressed in monetary terms per time unit (e.g. euros per day) and evolves in time, being gradually decreased until it is diminished upon re-establishment of the full network functionality (as seen in Fig. 8). In this study, this additional cost is defined as a function of the travel delay cost which is derived according to Goodwin (Goodwin 2004), plus a penalty cost assigned to every trip that is canceled after the earthquake occurrence either due to absence of an alternative route or due to a decrease in the travel demand (in case that a new OD matrix is used for the traffic assignment of the recovery phases):

$$EC_{p} = \sum\limits_{j = 1}^{J} {(V_{jp} t_{jp} - V_{j0} t_{j0} )} \cdot VOT + (TR_{0} - TR_{p} ) \cdot PC$$

(9)

where EC_p is the additional daily cost during phase p (per time unit). VOT is the value of time meaning the mean cost of 1-h delay per vehicle. V_jp is the traffic load in network link j during phase p.t_jp is the travel time in network link j during phase p.V_j0 is the traffic load in network link j before the earthquake occurrence. t_j0 is the travel time in network link j before the earthquake occurrence. J is the total number of network links. TR₀ is the total number of network trips completed before earthquake occurrence. TR_p is the total number of network trips completed during phase p.PC a penalty cost assigned to every trip canceled after the earthquake occurrence (Fig. 9).

Fig. 9

Evolution in time of the additional traffic cost associated with a sample damage scenario

7.2.3 Consequences vector

To further quantify the impact of the network disruption in the region of interest, a qualitative, also time-dependent, Consequence vector$\left\{ {ECO_{p} ,CON_{p} ,ENV_{p} } \right\}$ is introduced. This is a three-component vector used for assessing earthquake loss to the wider financial life of the affected area, the connectivity among various Points of Interest (POIs) and the environment.

Τhe vector is composed by three loss factors which vary between 0 and 1. The lower the factor (i.e., closer to zero), the higher the consequences. As anticipated, immediately before the earthquake the values of all the three factors are 1, drop suddenly when the earthquake event occurs and are gradually restored to 1 by the end of the recovery period. As mentioned, the Consequence Vector is populated by three distinct factors:

The wider economic loss factor$ECO_{p}$ is used to express the increase in the transportation and travel cost, the disturbance in the productive activities, the decrease in the tourist business and generally the earthquake consequences that have an economic impact to the society (Zhou et al. 2010; Cavalieri et al. 2012). It is wider than the additional travel cost of Eq. (9) but is again based on Goodwin’s approach and on the assumption that the economic consequences are proportional to the increase of travel time. The factor value for every phase of the recovery period is given as:

$$ECO_{p} = \frac{{\sum\nolimits_{j = 1}^{J} {(V_{j0} \cdot t_{j0} )} }}{{\sum\nolimits_{j = 1}^{J} {(V_{jp} \cdot t_{jp} )} }}$$

(10)

It is noted that in case of a link that is closed, traffic time tends to infinity hence, the relevant terms are taken equal to zero in the denominator sum.

The connectivity loss factor$CON_{p}$ refers to the consequences due to loss of access to points of interest. The POIs are structures or areas that it is crucial to be accessible after an earthquake such as hospitals, power stations, administrative buildings, ports and border checkpoints among others. It is noted that this index solely considers the accessibility and not the traffic load or travel time to each destination (Du and Peeta 2014). The calculation of the connectivity loss factor is based on the assumption that network functionality loss in a region where POIs are located affects the accessibility to the particular POIs. To derive the $CON_{p}$ it is necessary to define the number of POIs per activity zone (Fig. 2), so that a zone with a higher number of POIs can take higher weighting in the equation. Only zones that contain at least one important location are taken into account:

$$CON_{p} = \frac{{\sum\nolimits_{z = 1}^{Z} {(\gamma_{zp} \cdot N_{z} )} }}{{\sum\nolimits_{z = 1}^{Z} {N_{z} } }}$$

(11)

where $\gamma_{zp} = \frac{{L_{z0} - L_{zp} }}{{L_{z0} }}$ is the accessibility factor in activity zone z during phase p of the recovery period. L_z0 is total network length in activity zone z. L_zp is the total length of the network links that belong to activity zone z and experience a functionality loss during phase p. Ν_z is the total number of POIs located at activity zone z. Z is the total number of activity zones that include POIs.

The environmental loss factor$ENV_{p}$ assesses the consequences to environmental-sensitive areas due to an increase of traffic after an earthquake and the associated gas emissions, noise pollution and heavy vehicle transport (Dong et al. 2014b). Such areas may be national parks, protected forests, cultural heritage sites and regions with wild life passing routes. The calculation of the factor relies on the definition of the links that are environmentally sensitive (defined in the initial setup of the network, Fig. 2). Environmentally sensitive links may be the links that cross environmentally sensitive areas but also other links that have a high environmental impact themselves.

$$ENVp = \frac{{\sum\nolimits_{j = 1}^{V} {l_{j} } }}{{\sum\nolimits_{j = 1}^{V} {a_{jp} \cdot l_{j} } }}$$

(12)

where α_jp is the ratio between the traffic load at environmentally sensitive network link j during phase p to the pre-earthquake traffic load at that link. l_j the length of the environmentally sensitive network link j. V the total number of the environmentally sensitive network links.

The graphical representation of the Consequences vector$\left\{ {ECO_{p} ,\;CON_{p} ,\;ENV_{p} } \right\}$ as function of time illustrates the evolution of earthquake wider financial, connectivity and environmental impact until the functionality of the road network is fully restored (Fig. 10).

Fig. 10

“Consequences vector-time” correlated to a damage sample

7.3 Cumulative indicators

Εvery damage state that a key component i (bridge, overpass, slope etc.) may experience is associated to a structural repair cost ratio (Mackie and Stojadinovic 2006) defined as the repair cost over the total re-construction cost of that specific network component. It is necessary to provide such repair cost ratio indices, during the Network definition stage, for all four damage states prior to the analysis for each component class. Based on the repair cost ratio and the probability of attaining every damage state for each key component i, the estimated structural cost (ESC) due to the l^th IM distribution (seismic map) of hazard level k is derived as the sum of the cumulative direct cost of structural damage within the network due to that specific IM distribution:

$$ESC_{k,l} = \sum\limits_{i = 1}^{I} {TBC_{i} \cdot \sum\limits_{\theta = 1}^{4} ( RCR_{\theta }^{i} \cdot P_{{DS_{\theta } }}^{k,l,i} )}$$

(13)

where $TBC_{i}$ is the total cost of re-constructing key component i.$RCR_{1}^{\theta }$ are the repair cost ratios corresponding to DS_θ (θ ranges from 1 to 4 for slight to collapse DS respectively). $P_{{DS_{\theta } }}^{k,l,i}$ is the probability that the damage of the key component i exceeds DS_θ due to the seismic map k, l. I the total number of key network components.

On the other side, an earthquake-induced traffic cost (TC) can be calculated for every Monte Carlo simulated damage sample. This cost is the additional cost during the entire recovery period of that damage sample, and as such, it is the sum of the product of each phase duration, times the corresponding additional travel cost.

$$TC_{k,l,s} = \sum\limits_{p = 1}^{{P_{k,l,s} }} {EC_{{k,l,s_{,} p}} } \cdot t_{{k,l,s_{,} p}}$$

(14)

where $t_{k,l,s,p}$ is the duration of phase p of the sth damage sample derived from the seismic map k, l (i.e. the map associated to hazard level k and seismic segment l). $P_{k,l,s}$ is the total number of recovery phases associated with the sth damage sample derived from the seismic map k, l.

Subsequently, the estimated traffic cost (ETC) can be associated to every seismic map, as the mean of the costs calculated for the S Monte Carlo samples simulated from that map:

$$ETC_{k,l} = \frac{{\sum\nolimits_{s = 1}^{S} {TC_{k,l,s} } }}{S}$$

(15)

Overall, L number of total network costs (TNC) are computed as the sum of the estimated structural and traffic cost out of the L seismic maps associated to every hazard level k. This information is key in assessing and improving the resilience of the network as discussed in the following.

8 Decision-making

8.1 Available data for informed-decision making

Based on the above methodology and formulations, for each hazard level k that corresponds to the selected four annual rates of exceedance (e.g. 2%, 0.20%, 0.10% and 0.05%), the stakeholder has a set of indicators, namely:

• (a) L sets of cumulative, structural (ESC), traffic (ETC) and total network (TNC) costs, corresponding to the L seismic maps of each hazard level k. Notably every set of these costs has the same rate of exceedance, given that the L maps, are also equally probable.

• (b) A set of L×5 charts, plotting the time evolution of the five time-variant indicators presented in Sect. 7.2, namely, the additional network functionality $f_{p}$, traffic cost, as well as the qualitative financial, social and environmental loss factors as expressed by the Consequence Vector $\left\{ {ECO_{p} ,CON_{p} ,ENV_{p} } \right\}$. Each of these charts consists of S curves. It is recalled that S represents the size of the MC sampling resulting from the probabilistic assessment of structural/geotechnical damage of the key network components given the distribution of IM in space for a seismic segment l.

Notably, the above set of computed results constitutes a very detailed prediction of network resilience. However, its size is such that does not facilitate easy decision making. For instance, in case of L = 10 seismic segments (possible epicenters) that generated 10 seismic maps of IM distribution, 50 different charts (10 × 5 time-variant indicators) would be made available. For this reason, the following decision-making strategy is adopted.

8.2 Strategy for improving road network resilience

The proposed strategy involves a number of steps that need to be taken into account in order to condense the above information into an applicable set of decision-making criteria. More precisely, the stakeholders, with the aid of the adhoc software developed, shall be able to:

Step 1: First, identify, for each hazard level k (annual rates of exceedance of 2%, 0.20%, 0.10% and 0.05%), amongst the L seismic segments (and the subsequent, equally probable, L maps derived according to Sect. 4), the critical seismic segment that leads to the highest, cumulative, total network cost (TNC_k) or, in other words, has the highest contribution to the overall loss resulting from equiprobable possible losses. The total cost is defined as the sum of the structural/geotechnical damage (ESC_k,l, Eq. 13) of the key components (bridge, slopes, tunnels etc.) and the additional traffic delay cost (ETC_k,l Eq. 15) resulting from the partial operation of the network after an earthquake event:

$$TNC_{k} = \max_{l \in [1,L]}^{{}} \left( {ESC_{k,l} + ETC_{k,l} } \right)$$

(16)

It is noted that adoption of the critical seismic segment as the criterion for network resilience assessment is one of the possible approaches that can be made for the stakeholder, whose decision depends on its overall attitude toward risk. For instance, some stakeholders may be more interested in adopting a strategy that faces in a weighted way the whole range of possible future earthquakes (e.g. the L seismic maps) rather than is based on the worst case.

Step 2: Plot for every hazard level the charts for the time-variant indicators (similar to Fig. 11) but just associated with the corresponding critical seismic segment in order to identify the time that takes for the functionality of the “as built” roadway network to be fully restored, as well as plot the evolution of the Consequences vector.

Fig. 11

Evolution of one of the time-variant indicators (network functionality) for a sample seismic map

Step 3: Convert, the above 5 charts (inclusive S curves each) into easy to compare form. This is needed because, due to the form of the time-variant indicators, it is not possible to be directly used to comparatively assess alternative risk mitigation strategies. In order to facilitate comparison therefore, a normalised Resilience Index Vector (RI) is introduced, namely {RI^fun, RI^eco, RI^con and RI^env} matching to network functionality and the Consequence Vector terms $\left\{ {f_{p} ,ECO_{p} ,CON_{p} ,ENV_{p} } \right\}$. It is noted that the fifth time-variant indicator which is the additional cost $EC_{p}$, is cumulatively accounted for in the estimated traffic cost and as such, it is neglected herein. Assuming a reference period of recovery t_ref (e.g., 500 days, Fig. 12) within which integration of Resilience is made, the Resilience Index Vector terms are derived as (Cimellaro et al. 2016):

$$RI_{k} = \frac{{\sum\nolimits_{s = 1}^{S} {} (t_{ref} - \sum\nolimits_{p = 1}^{{P_{{k,l_{wk} ,s}} }} {t_{{k,l_{wk} ,s,p}} \cdot (1 - QI_{{k,l_{wk} ,s,p}} )} )}}{{S \cdot t_{ref} }}$$

(17)

where $QI_{{k,l_{wk} ,s,p}}$ is the value of the time-variant indicator $\left\{ {f_{p} ,ECO_{p} ,CON_{p} ,ENV_{p} } \right\}$ associated to the pth phase of the sth damage sample, derived from the map associated with the critical seismic segment l_wk for hazard level k. As shown in Fig. 12, the grey area quantifies the integration of the loss over the target recovery period, while the green one reflects the (desirable) resilience.

Fig. 12

Schematic representation of the areas that define the Resilience Index term for the economic cost (ECO)

Step 4: Assess the Total Network Cost TNC of Step 1 and the time evolution of the Resilience Index Vector of Step 3 versus pre-defined Resilience Objectives (for instance, against is the tolerable total cost, downtime or consequence for the particular network after an earthquake with a given rate of exceedance λ_k), as shown in Fig. 13.

Fig. 13

Detailed workflow for the proposed decision making process

Step 5: In case that the Resilience Objectives are not met, the stakeholder needs to revert to the evolution charts of Step 2 to identify the critical network components (bridges, overpasses, slopes, etc.) whose failure has the highest impact on the total cost and resilience of the overall network. This can be made for instance by focusing on the components which appear to be critical in the highest number of MC samples for the (critical) seismic segment examined (Fig. 13).

Step 6: Form a set of Resilience Improvement Strategies that consist of alternative loss mitigating measures (retrofit schemes, planning for improved post-earthquake response, construction of alternative roads for detours and network robustness, among others) applied to identified critical components and estimate their investment (e.g. initial) cost. Given that there are several alternative measures that can be taken and also that each resilience improvement strategy may consider different sets (combinations) of critical components, a number of different strategies can be formed by the stakeholder.

Step 7: Rerun the analysis of the network for each one of the resilience improvement strategies formed and comparatively illustrate the total network cost and Resilience Indices per scenario (Fig. 13).

Step 8: Reassess versus the Resilience Objectives and select the appropriate resilience improvement strategy.

A disclaimer is made that, as noted already, this procedure refers to the most critical seismic segment but can be repeated for additional segments with high impact as well, as per the stakeholders approach to risk. A more refined optimization scheme for aggregating the impact of all seismic segments and back-identifying the critical network components with the highest contribution to this overall risk would be more rigorous, however, it is deemed that it further increases the computational cost and hence it is out of the scope of this study. It is also noted that the software developed and presented in the following pilot study, facilitates at a great extent the execution and visualization of steps 1–8 thus manual repetition of this approach for additional seismic segments is quite handy. It is acknowledged that this process can be utilized for the identification of the best Resilience Improvement Strategy within the group of candidate strategies formed in Step 6 rather for the definition of the overall optimum Resilience Improvement Strategy which would require the application of an optimization algorithm and lead to a prohibitive computational cost.

9 Software development and case study

The above holistic probabilistic framework is materialized as a standalone, GIS-based interactive freeware that is made available to the engineering community (www.retis-risk.eu). This software implements all framework stages in an efficient Matlab code and incorporates open-source traffic assignment engine DTALite (Zhou et al. 2014) for computing the pre- and post- earthquake traffic flows by analyzing the normal traffic conditions and the generated traffic scenarios, respectively.

The same software and framework is been applied herein for demonstration purposes and the case of a simple, idealized, road network that consists of 12 uni-directional links as shown in Fig. 14. A total number of 9 key components are identified, classified into three general classes of identical vulnerability (namely, class 1: four span R/C bridges with continuous deck and low damping rubber bearings at the deck-abutment connections, class 2: two span R/C bridges with monolithic deck-abutment connections and class 3: typical shallow R/C tunnels with circular cross section of 5 m radius). A traffic demand of 20,000 cars per day and destination is assumed to be originating from nodes 9 and 8 towards nodes 3 and 4, while the drivers tend to choose the fastest route, i.e., 9-7-6-3 and 8-5-4, which correspond to the shortest paths with the higher speed limits. Four hazard levels are examined (1, 2, 3 and 4 with annual exceedance rate of 2%, 0.20%, 0.10% and 0.05%, respectively), each one defined by ten seismic maps, (i.e., ten seismic segments are considered). For Damage States DS1 to DS4, a repair time of 0, 7, 150 and 450 days, respectively is defined for all the 9 key components and the corresponding restoration matrix with 9 identical rows was created according to Sect. 6. Similarly, Damage State-specific repair cost ratios identical for the 9 key components are defined according to Werner et al. (2006) (i.e., RCR index value is taken 0.03, 0.25, 0.75 and 1 for the four Damage States).

Fig. 14

Case study: network topology

As described in the methodology, the probability of each key network component to experience a certain level of damage given the Intensity Measure at the site(s) of interest for each distinct seismic map and each hazard level, leads to a Monte Carlo analysis of open/closed components coupled with traffic and consequence analysis. Post-earthquake damage samples derived in this way are then analyzed into traffic phases to account for the gradual restoration of network functionality. Figure 15 shows such a decomposition of an indicative damage sample into phases (in red the links that are closed due to closed key network components).

Fig. 15

Recovery analysis of a damage sample

It is noted that during Phase 1 and 2, network functionality is zero because both the two pre-earthquake traffic paths remain closed. Moreover during the same period, additional traffic cost remains unchanged since the opening of link 7 at the end of Phase 1 does not have any beneficial effect to the network in terms of travel paths availability. The restoration of path 8-5-4 at the end of phase 2 leads to a stepwise network functionality increase and to a corresponding drop to the additional traffic cost.

The above calculation for all possible samples permits the quantification of the critical seismic segment, the subsequent map, the Total Network Cost, the evolution of the five resilience indicators and finally the Resilience Index Vector for the worst case scenario as described in Steps 1–6 of the previous section.

The Resilience assessment is then repeated (Step 7) assuming that specific measures are taken to improve the post-earthquake recovery response. Here, the adopted resilience improvement strategy involves better organized recovery teams that can reduce the time needed to recover all the damaged key network components. In this case, the improved restoration matrix corresponds to repair times of 0, 4, 100 and 300 days for Damage States DS1 to DS4, respectively for all the damaged key components, as opposed to the repair times of 0, 7, 150 and 450 days of the “as-built” network. Repeating all steps of the analysis, the above resilience improvement strategy is found to lead to a significantly lower Total Network Cost, due to the reduction of the traffic cost. More precisely, comparing the TNC between the two decision making diagrams and for the four annual exceedance rates of 2%, 0.20%, 0.10% and 0.05%, it is reduced from {110, 260, 290, 340} to {75, 175, 185, 240} million euros, respectively, which roughly corresponds to a reduction of approximately 30% overall (Figs. 16, 17). On the other hand, structural cost (ESC) remains the same in the two cases for all four hazard levels considered. This is due to the fact that only recovery times are improved, while damage probabilities and repair/cost ratios remain the same. Apparently, several other “nominee” mitigation strategies can be examined as per Step 8, such as retrofit of key network components for instance, until the optimum approach for improving the network resilience is identified.

Fig. 16

Case 1: Decision making diagram for the reference case of the “as built” roadway network and the most critical seismic segment

Fig. 17

Case 2: Decision making diagram for the a resilience improvement strategy that consists of a faster recovery of all the damaged key network components

10 Conclusions

In this study, a holistic framework is presented for the multi-criterion assessment and management of the seismic risk and resilience of roadway networks. Different sources of uncertainty that contribute to the overall network seismic risk, namely, hazard and vulnerability, coupled with consequences analysis are accounted for and integral aspects of resilience such as network functionality and post-earthquake time-dimension are integrated into the overall process. Several, novel, time-variant and cumulative indicators are introduced for reflecting time evolution of losses and providing a compact and more easily to interpret estimate of the extent of the losses incurred by the community after a major seismic event as well as throughout the entire recovery period. The holistic estimate of the direct (structural) and indirect (traffic) monetary loss as well the wider financial, network connectivity and environmental impact of hazard levels with different annual rate of exceedance is used as a means to describe the evolution of network functionality and resilience taking explicitly into account its gradual recovery. The proposed framework is implemented into a GIS-based software and constitutes a useful decision-making tool for the stakeholders, to quantify and improve the resilience of their roadway network, by testing and adopting the most appropriate alternative resilience improvement strategy.