#### Abstract

Earthquakes can have significant impacts on transportation networks because of the physical damage they can cause to bridges. Hence, it is essential to assess the seismic risk of a bridge transportation network accurately. However, this is a challenging task because it requires estimating the performance of a bridge transportation network at the system level. Moreover, it is necessary to deal with various possible earthquake scenarios and the associated damage states of component bridges considering the uncertainty of earthquake locations and magnitudes. To overcome these challenges, this study proposes a new method of system-level seismic risk assessment for bridge transportation networks employing probabilistic seismic hazard analysis (PSHA). The proposed method consists of three steps: (1) seismic fragility estimation of the bridges based on PSHA; (2) system-level performance estimation using a matrix-based framework; and (3) seismic risk assessment based on the total probability theorem. In the proposed method, PSHA enables the seismic fragility estimation of the component bridges considering the uncertainty of earthquake locations and magnitudes, and it is systemically used to carry out a posthazard bridge network flow capacity analysis by employing the matrix-based framework. The proposed method provides statistical moments of the network performance and component importance measures, which can be used by decision makers to reduce the seismic risk of a target area. To test the proposed method, it is applied to a numerical example of an actual transportation network in South Korea. In the seismic risk assessment of the example, PSHA is successfully integrated with the matrix-based framework to perform system reliability analysis in a computationally efficient manner.

#### 1. Introduction

Natural disasters can have significant impacts on infrastructure, such as transportation, electricity, gas, and water distribution networks, causing structural damage and economic losses in commercial and residential activities. In particular, earthquakes that occurred in recent years caused serious physical damage and disruption to transportation networks. Damage to transportation system is a major concern because it can impose extra burden on the other lifelines [1, 2]. One of the most significant impacts of earthquakes is the disconnection in bridge transportation networks, which can impede posthazard emergency responses, such as the movement of emergency vehicles. This is because bridges are one of the most critical components of transportation networks acting as “bottlenecks”: the structural failure of a bridge can interfere with traffic flow and decrease network performance [3]. Hence, it is essential to assess the seismic risk of a bridge transportation system and predict the posthazard performance accurately.

The objective of a seismic risk assessment is to obtain useful information for risk-informed decision making on seismic hazard mitigation and disaster management. For this reason, seismic risk assessment of critical infrastructure has been conducted extensively. Nuti et al. [4] proposed a methodology for the reliability assessment of electric power, water, and road systems, not considering the interdependence between the networks, whereas Poljnašek et al. [5] proposed a method for gas and electricity transmission networks considering the increased vulnerability due to interdependency. Dueñas-Osorio et al. [6] evaluated seismic responses considering the interdependency of the water and power networks in Shelby County, Tennessee 56 area, and proposed a method to apply mitigation action efficiently. With regard to transportation networks, Kiremidjian et al. [7] evaluated the risk from earthquakes to a transportation system in terms of direct loss from damage to bridges by studying the San Francisco Bay Area. Moreover, various studies have proposed postearthquake flow capacity models for evaluating the impact of seismic events and the functionality of the networks [8, 9].

For accurate seismic risk assessment, it is necessary to deal with the uncertainties associated with earthquakes and infrastructure responses [10], and sampling-based approaches are often used for this purpose. For example, Sánchez-Silva et al. [11] proposed a model for the efficient assignment of resources for the reliability of a transportation network, and Rokneddin et al. [12] conducted a system reliability analysis and determined the importance of bridges in a transportation network through Origin-Destination (O-D) connectivity based on Markov Chain Monte Carlo (MCMC) simulations. Shinozuka et al. [13] proposed a probabilistic model to determine the effect of repair of bridge damage on the improvement of the transportation network and measure the system performance in terms of a “Drivers Delay” index. Significant research efforts have been dedicated to estimating the seismic performance using sampling-based approaches; however, the approaches may increase the computational costs and can be time-consuming tasks to deal with possible component failure scenarios. In particular, when a sampling-based approach is used for the seismic risk assessment, in which the earthquake uncertainty needs to be considered, owing to its sampling nature, the analysis may require a large number of simulations to achieve an acceptable level of convergence in the results [14].

To overcome these challenges, a few nonsampling-based approaches have been developed. Li and He [15] proposed a recursive decomposition algorithm for seismic reliability evaluation to compute the probabilities of disconnections in a network, while Kang et al. [16] and Kang and Song [17] proposed a new nonsampling-based system reliability analysis method, the matrix-based system reliability (MSR) method. The matrix-based framework of the MSR method enables rapid calculation of multiple probability scenarios and separation of network and vulnerability analyses. Using the MSR method, Lee et al. [18] estimated the posthazard flow capacity of a bridge transportation network considering ten seismically vulnerable bridges, and the MSR method was successfully applied for evaluating the system reliability within a network.

Meanwhile, most of the previous studies conducted seismic risk assessments without probabilistic seismic hazard analysis (PSHA). PSHA can play a crucial role in considering earthquake uncertainty because it is useful for determining the uncertainty in earthquake locations and magnitudes in a target region [19]. However, it is a challenging task to introduce PSHA into seismic risk assessment because it requires dealing with a large number of possible earthquake scenarios and quantifying the performance of a bridge transportation network.

For this reason, this paper proposes a new method of system-level seismic risk assessment for bridge transportation networks employing PSHA. The proposed method consists of three steps: (1) seismic fragility estimation of the bridges based on PSHA; (2) system-level performance estimation using the matrix-based framework of the MSR method; and (3) seismic risk assessment based on the total probability theorem. In the proposed method, PSHA enables the seismic fragility estimation of the component bridges considering the uncertainty of earthquake locations and magnitudes, and it is systemically used to carry out a posthazard bridge network flow capacity analysis employing the matrix-based framework of the MSR method. The proposed method provides statistical moments of the network performance and component importance measures, which are useful for decision makers to reduce the seismic risk of the target area.

#### 2. Theoretical Background

##### 2.1. Probabilistic Seismic Hazard Analysis (PSHA)

In PSHA, a seismic hazard is defined as a physical phenomenon such as ground shaking or failure affected by an earthquake which can cause serious effects on human activities [20]. The main goal of seismic hazard analysis is to refine the understanding of earthquake magnitudes and the corresponding intensity of ground shaking. In particular, PSHA allows a better insight into earthquake generation and seismic effects on a region by quantifying the uncertainties and estimating the distribution of earthquake occurrences [11]. For the precise estimation of the risk caused by a seismic event at a particular area, seismic hazard should be analyzed probabilistically, by considering uncertainties in earthquake locations and magnitudes.

In the proposed method, the uncertainty of earthquake locations and magnitudes are identified using PSHA based on past earthquake records. The uncertainty of earthquake locations can be considered by using the epicenters of past earthquakes for the seismic risk assessment. Meanwhile, to account for the uncertainty of earthquake magnitudes, a series of modeling procedures is required. In this paper, one such modeling procedure is briefly introduced, whereas more details on PSHA can be found in Baker [19].

###### 2.1.1. Earthquake Magnitude Uncertainty Modeling

To account for the uncertainty of earthquake magnitudes, as the first step, the occurrence rate of earthquakes in a target region is assumed to follow the Gutenberg-Richter (G-R) recurrence law [21] based on the following logarithm using base 10:where *m* is the specific earthquake magnitude of interest, *λ*_{m} is the occurrence rate of earthquakes of magnitudes greater than *m*, and *a* and *b* are the constants which are referred to as G-R recurrence parameters and can be determined from past earthquake records. This G-R recurrence law is based on the Poisson process assumptions [22]: (1) an earthquake can occur at any random time, (2) the occurrence of an earthquake in a given time is statistically independent of that at any other time, (3) the occurrence probability of an earthquake in a small interval is proportional to the interval, and (4) the probability of experiencing more than one earthquake in a small interval is negligible.

When the minimum and maximum magnitudes are determined using equation (1), the cumulative distribution function (CDF) of the earthquake magnitude can be derived as [19]where denotes the CDF of a random variable, *M* is the earthquake magnitude, and *m*_{min} and *m*_{max} are the minimum and maximum of earthquake magnitude, respectively. By differentiating equation (2), the probability density function (PDF) can be obtained aswhere denotes the PDF of a random variable. This bounded PDF is termed the *bounded Gutenberg-Richter recurrence law* [19], and it can be obtained based on the frequency of past earthquakes with varying magnitudes.

Once the bounded PDF is obtained, to generate possible earthquake scenarios with varying magnitudes as an input of seismic risk assessment, the continuous distribution of earthquake magnitudes needs to be converted into a discrete set of magnitudes of interest. The probabilities of occurrence, according to this discrete set of magnitudes, represent only a partial distribution of magnitude at a site. Then, a normalizing process that divides all of the cumulated values by their sum is needed so that the sum of the probability distribution in the partial magnitudes becomes 1.0.

###### 2.1.2. Seismic Attenuation Model

After modeling the probability distribution of earthquake magnitudes, the associated distances from the earthquake source to the bridges (i.e., source-to-site distances) of the target transportation network and ground motion intensities need to be analyzed. Given the earthquake magnitude and location, the ground motion intensities at different bridges must be analyzed based on a seismic attenuation model. A representative model is the ground motion prediction equation (GMPE). In this equation, the ground motion intensity is expressed as a function of several parameters including earthquake magnitude, distance, and local site effects [23]. For seismic event *i* recorded at the site *j*, the general form of GMPE considering the total variability of the ground motion is given as [24]where *Y*_{ij} represents the response variable, such as peak ground acceleration (PGA), peak ground velocity (PGV), and spectral acceleration (SA), which often corresponds to the logarithm (natural or common), *M*_{i} is the earthquake magnitude of the event, *R*_{ij} is the distance between the epicenter of event *i* and the site *j*, *ξ*_{ij} is the geomorphic factor affecting the ground motion, is the mean of the response variable, and *η*_{i} and *ε*_{ij} are the inter- and intraevent parameters representing the uncertainty of the ground motion. The parameter *η*_{i} represents the uncertainty of the ground motion due to the inherent earthquake itself, termed the *Earthquake-to-Earthquake* variability, whereas the parameter *ε*_{ij} denotes the uncertainty of the ground motion because of the energy paths and geological characteristics, termed the *Site-to-Site* variability [24].

In this study, as in HAZUS-MH [25], the structural vulnerabilities of a bridge are described by the probability conditioned to the ground motion intensity which is expressed in terms of SA. In addition, the GMPE proposed by Emolo et al. [24] is introduced. The GMPE was derived using 222 earthquakes recorded at 132 stations in South Korea using the nonlinear mixed effects regression analysis. Moreover, the equation includes both fixed and random effects accounting for inter- and intraevent residual values. Using the equation, the mean of the response variable in equation (4) can be given aswhere SA_{ij} is the spectral acceleration caused by an earthquake event *i* at site *j*, *h* is the focal depth, *c*_{k} are the regression coefficients, and *s* is the station dummy variable, which assumes a value of −1, 0, or 1. The dummy variable depends on the sign of the mean residual (negative, zero, or positive, respectively), and it is generally chosen by the seismological observatory.

In addition, the inter- and intraevents parameters in equation (4) can be expressed by the following equations [26–29]:where *σ*_{η}^{2} and *σ*_{ε}^{2} are the inter- and intrastandard residuals, respectively, Δ_{ij} is the distance between the epicenter of event *i* and the site *j*, and is the spatial correlation equation. For the special correlation equation, in this study, the following equation suggested by Goda and Hong [26] is introduced:

As described so far, PSHA enables the consideration of earthquake uncertainty. For the given location and magnitude of an earthquake, the ground shaking intensity of individual bridges is expressed in terms of SA using the GMPE in equation (4). Then, the probabilities of several damage states of bridges can be provided by fragility curves. Subsequently, it becomes possible to estimate the performance of the bridge transportation network, which requires system reliability analysis.

##### 2.2. Matrix-Based System Reliability (MSR) Method

As previously mentioned, sampling-based approaches are often used for the performance estimation of lifeline networks subjected to earthquakes. However, these approaches may prevent rapid risk assessment and make providing useful probabilistic measures difficult [18]. To overcome these, the MSR method was recently proposed and successfully applied to evaluate the postearthquake performance of a bridge transportation network in terms of the disconnection probability [16] and the maximum flow capacity [18]. The MSR method provides a matrix-based framework of system reliability analysis in which two tasks of “system event description” and “probability calculation” are performed separately. This enables efficient evaluations of the posthazard capacity under various failure scenarios without performing deterministic analyses of the maximum flow capacity repeatedly. In this paper, the MSR method is briefly introduced, and more details can be found in Lee et al. [18].

A bridge transportation network is considered, consisting of *n*_{b} bridges, each of which has *n*_{d} distinct damage states. When assuming all of the component bridges are statistically independent, there is a total of damage scenarios, which are mutually exclusive and collectively exhaustive (MECE) events. Due to the mutual exclusiveness, the probability of the system event , i.e., is the sum of the probabilities that belong to the system event. Therefore, can be easily computed by the inner product of the two vectors, probability and event vector [16, 18].

Let , , , indicate the probability that the bridge enters the damage state. Then, the probability of dealing with multiple damage states can be expressed by the following sequential matrix calculations:where **p** denotes the probability vector for all possible damage scenario of the system, *d*_{i} is the damage state of the bridge (where ), and denotes the probability of a damage state system with the numbers in the subscript. The subscript signifies the component damage states. For example, the second row, , means all of the components are in the first damage state except the first component which is in the second damage state.

In the MSR method, to carry out a posthazard flow capacity analysis, a certain corresponding quantity can be estimated by using the matrix-based framework. Therefore, a new column vector **q**, termed the “quantity vector,” is constructed, which has same size as **p** in equation (8). For all damage states of the system, the quantities are generalized towhere **q** denotes the quantity vector, denotes the performance quantity of the system with damage states as the subscript, and is the postearthquake capacity of the network. In this paper, the maximum flow capacity, i.e., the maximum number of vehicles that can pass per unit time, is determined as a measure of network flow capacity [30, 31]. For each maximum flow capacity corresponding to damage states of the network, the MATLAB® version of Boost Graph Library [32, 33] is used.

As results of the probability vector **p** in equation (8) and the quantity vector **q** in equation (9), the statistical parameters such as the mean, variance, and coefficient of variation (c.o.v.) can be obtained as follows:where , , and are the mean, variance, and c.o.v. of the maximum flow capacity of the network *Q*, respectively. In equation (10), “” denotes the element-by-element multiplication.

A potential drawback of the MSR method is that the sizes of the vectors and matrices increase exponentially as the number of components (e.g., the number of bridges in a transportation network) increases [14, 16, 34, 35]. However, it was also reported that this drawback could be overcome by the following: (1) transforming a large MSR problem into multiple smaller problems using a multiscale approach [14, 36] or (2) integrating the branch-and-bound method with the MSR method to prioritize important system states and evaluate only their probabilities [35].

Furthermore, an importance measure of components can also be evaluated using a matrix-based formulation. For regional authorities who need to make decisions for allocating budgets and other resources, it is important to identify the important locations of a transportation network [37]. To evaluate the relative importance of components, several importance measures have been developed in system engineering. In this study, the reduction factor (RF), a new version of importance measure, which was proposed by Lee et al. [18], is adopted. RF computes the reduction of the expected mean of the maximum flow capacity by the observed event *E*_{obs}, e.g., the failure of a bridge, which can quantify the relative importance of component bridges. The computation of the proposed RF is as follows:where represents the conditional mean of the maximum flow capacity given an observed event, *E*_{obs}, and is the probability vector constructed using the conditional probabilities of components given by *E*_{obs}.

#### 3. Proposed Method

##### 3.1. Matrix-Based Framework in the Proposed Method

The proposed method aims at estimating the performance of a bridge transportation network after a seismic event considering earthquake uncertainty. For this purpose, the proposed method introduces PSHA to deal with the uncertainty of earthquake magnitudes and locations. The performance of a bridge transportation network which is subjected to earthquakes also becomes uncertain, and the proposed method employs the matrix-based framework of the MSR method to evaluate this uncertainty at the system level.

The MSR method was originally developed to perform system reliability analyses of various structures [16, 38]. However, its convenient matrix-based framework enables the system reliability analysis of lifelines, and it was successfully applied to evaluate the system-level performance of bridge transportation networks [18, 35]. In particular, it was applied to estimate the posthazard flow capacity of a bridge transportation network considering the structural deterioration of bridges [18]. In the previous research, the time-varying bridge fragilities causing structural deterioration could be computed efficiently using the MSR method, while the corresponding flow capacities were estimated using a maximum flow capacity analysis algorithm. The matrix-based framework separates probability calculation and network flow capacity analysis, which enabled performing extensive parametric studies and time-varying posthazard flow capacity analyses without repeated network flow capacity analyses. Similarly, in this research, the matrix-based framework of the MSR method is introduced to deal with the earthquake uncertainty.

As stated previously, by assuming that all the components are statistically independent, the basic MECE events can be easily computed by using the matrix calculation proposed in equation (8). Nevertheless, many reliability problems of network systems are composed of various components which are statistically dependent. However, using the concept of *common source random variable* (CSRV), a system event can be described by the combination of component MECE events which are conditionally independent.

In the proposed method, the earthquake magnitude and location are introduced as CSRVs that influence the damage states of individual bridges. When the uncertainty of earthquake magnitudes and locations is characterized by PSHA and defined as CSRVs, using the total probability theorem, the probability of a system can be expressed aswhere is the system event of interest, is the earthquake magnitude, *L* is the earthquake location, and is the joint PDF of the earthquake magnitude and location when and . Using the conditional PDF, the joint PDF in equation (12) can be changed towhere is the marginal PDF of the earthquake location and is the marginal PDF of the earthquake magnitude conditioned to .

Since the damage states of individual bridges become conditionally statistically independent given the earthquake magnitude and location, the conditional probability vector can be constructed using equations (8), (12), and (13):where denotes the conditional probability of a component’s damage state with the numbers in the subscript. At the same time, it is noticeable that, unlike the probability vector, the quantity vector **q** can remain the same as in equation (9) because the matrix-based framework of the MSR method enables the separation of the construction of the probability and quantity vectors.

Hence, using the total probability theorem, the statistical parameters can be obtained as

Similarly, the reduction factor RF can be calculated aswhere is the conditional probability vector which can be constructed using equation (11).

##### 3.2. Seismic Risk Assessment Employing the Proposed Method

The proposed method consists of three steps: (1) seismic fragility estimation of the bridges based on PSHA, (2) system-level performance estimation using the matrix-based framework of the MSR method, and (3) seismic risk assessment based on the total probability theorem. In the proposed method, PSHA enables the seismic fragility estimation of the component bridges considering the uncertainty of earthquake locations and magnitudes, and it is systemically used to carry out a postearthquake bridge network flow capacity analysis employing the MSR method.

The MSR method provides an efficient framework for seismic risk assessment which performs separate calculations for seismic hazard and network flow capacity analyses. PSHA contains the investigation of earthquake generation using several proposed relations provided in equations (1)–(7), in order to finally obtain the probabilities of structural damage scenarios considering bridge fragilities given the earthquake magnitude and location. When the corresponding network flow capacities are identified and the probability vector and quantity vector are constructed, the MSR method enables the calculation of various risk-informed measures and statistical use of the total probability theorem. Figure 1 shows a flow chart of the seismic risk assessment proposed in this research.

To start the analysis employing the proposed method, it is first necessary to collect input data in the form of information related to the study area. The input data can be classified into three groups: exposure data, hazard data, and structural vulnerability data. The exposure data contains the topology information of the target transportation network such as the number of nodes and links in the target region. The past earthquake data are also necessary to identify the earthquake uncertainty using PSHA. The hazard data contain the mean occurrence rate of earthquakes with varying magnitudes and the seismic attenuation law (i.e., GMPE) with the spatial correlations described in equations (4)–(7). Lastly, the structural vulnerability data include the determination of damage states of individual bridges. For this task, the fragility curve parameters provided in HAZUS-MH [25] are adopted in this study.

The next step is to perform seismic risk assessment employing PSHA and the matrix-based framework of the MSR method. First, from the past earthquake data, earthquake sources are analyzed using the G-R relation law. The purpose of PSHA is to identify uncertainties related to the earthquake itself and calculate the resulting intensity of ground motion. Therefore, final results are the PDF of earthquake magnitude and the probability of the different damage states of the component bridges. After performing the seismic hazard analysis, the proposed approach predicts the posthazard flow capacity of a transportation network for given magnitudes and locations of earthquake. In this process, the quantity vector **q** is constructed for all possible combinations of bridge damage states using the maximum flow capacity analysis. Next, the conditional probability vector in equation (14) is constructed based on the probabilities of bridge damage states.

Lastly, the conditional probability vector construction is repeated for various earthquake scenarios with different earthquake magnitudes (*M*) and locations (*L*). Then, the performance of the transportation network can be estimated in terms of the statistical moments of the maximum flow capacity using equation (15). In addition, the importance measure of RF can be estimated for all component bridges using equation (16). It is noteworthy that only the tasks in the dotted boxes in Figure 1 need to be repeated but the maximum flow capacity analysis which is computationally expensive does not.

#### 4. Numerical Example

##### 4.1. Problem Description

###### 4.1.1. Target Bridge Transportation Network

To test the proposed method, it is applied to an actual transportation network around Pohang city, South Korea. The study area is located in the southeast of the Korean Peninsula, and it experienced a 5.4-magnitude earthquake in 2017 [39]. Although South Korea was known to have relatively low seismic risk, this earthquake and its aftershocks raised lasting concerns which this study aims to help address.

Figure 2 illustrates the topology of the target network, consisting of 37 nodes (red solid circles) and 46 links (blue lines). It is a network of expressways and national routes in and around Pohang city and includes ten (i.e., *n*_{b} = 10) relatively long bridges (colored black) in the area. In this example, the objective is to measure the maximum flow capacity of the network after a seismic event, and Nodes 3 and 30 represent the origin and destination, respectively. Table 1 shows the locations (i.e., connecting nodes) of the ten bridges, and Table 2 shows their structural information. Table 3 shows the maximum flow capacities of the links (given in the number of passing vehicles per hour). These assumptions enable the performance assessment of this network in terms of the maximum flow capacity using a node-to-node flow capacity analysis.

To account for uncertainty in the seismic damage states of bridges, seismic fragility curves are introduced. Seismic fragility is defined as the conditional probability that the demand of a structure exceeds a specified threshold for a given earthquake intensity [40–42], and seismic fragility curves are often used for setting retrofit and repair priorities of bridges after an unexpected and disastrous event [43]. In this study, SA, which can be calculated by the GMPE provided in equation (4), is introduced as the earthquake intensity. In addition, seismic fragility curves are obtained from HAZUS-MH [25] in which bridges are classified by several factors including their length, type, and seismic design methods, and the corresponding seismic fragility curves are provided. Similarly, the seismic fragility curves of the ten considered bridges are determined based on the structural information presented in Table 3.

In this example, five damage states of no, slight, moderate, extensive, and complete damage are assumed (i.e., *n*_{d} = 5). For the maximum flow capacity of a damaged bridge, quarter-based maximum flow capacities, which have been used in previous studies [44, 45], are adopted, as described in Table 4. In the table, 100%, 75%, 50%, 25%, and 0% of the original flow capacities represent the remaining maximum flow capacity of the bridges under the five damage states. However, the accuracy of the analysis can be improved if more realistic data on the maximum flow capacity of a damaged bridge are made available. For each combination of bridge damage states in the probability vector, in the MSR method, these values of maximum flow capacity are assigned to the corresponding bridges during the maximum flow capacity analysis so that the quantity vector can be constructed.

###### 4.1.2. Earthquake Uncertainty Identification

To identify the earthquake uncertainty in the target region, past earthquake data (with magnitude, *M*_{L}, greater than or equal to 3.0) were collected from the Korea Meteorological Association (KMA) website from January 1st, 1918, to August 22nd, 2018. As a result, 20 earthquake records were collected. Table 5 presents the earthquake information, and Figure 3 shows the locations of the 20 epicenters and the ten bridges. In this example, it is assumed that all of the past earthquake epicenters have the same likelihood of suffering an earthquake again, and the distances between the epicenters and the bridge locations are calculated based on their location information.

Based on these earthquake records, the earthquake uncertainty is identified using the G-R law given in equation (1). For the task, first, the range of earthquake magnitudes is divided by an interval of 0.1, and the rate at which each earthquake magnitude exceeds is calculated based on the past earthquake records. Next, the parameters of *a* and *b* in the G-R law (i.e., equation (1)) are calculated using a regression analysis. In this study, a MATLAB function, *lsqcurvefit*, which can solve nonlinear curve-fitting problems using the least-squares method, is employed, and the parameters *a* and *b* in equation (1) are obtained to be 2.167 and 0.699, respectively. Figure 4 shows the exceedance rate of the earthquake magnitude and the regression function. In the figure, it is noteworthy that the exceedance rate is plotted using a log scale, which makes the regression function linear.

Then, using equation (3), the bounded PDF of the earthquake magnitude can be constructed, and the corresponding occurrence probabilities with varying earthquake magnitudes in the range can be obtained, which is shown in Figure 5.

Once the occurrence probability of earthquake magnitude is constructed, the earthquake intensities at the ten bridge locations are calculated using equation (4). First of all, the mean of the response variable can be calculated using equation (5). For the calculation, the regression coefficients *c*_{k} are assumed to be −5.15, 0.95, −0.92, 6.8, −0.0003, and 0.208, respectively [24], and the station dummy variable *s* is assumed to be −1, which was recommended by the seismological observatory (station code: PHA2) of the region.

##### 4.2. Analysis Results

First, for each set of earthquake magnitude and location, the proposed method provides the probability mass function (PMF) and cumulative distribution function (CDF) of the maximum flow capacity of the network. For example, for EQ8’s epicenter, Figure 6 shows the PMF of the maximum flow capacity for a magnitude of 7.5. In addition, for the same epicenter, Figure 7 shows the CDFs with magnitudes of 5.0, 6.0, and 7.5 and an uncertain magnitude (*M*), and it can be seen that the CDF values increase as the earthquake magnitude increases. Furthermore, the probability that the maximum flow capacity of the network is smaller than a given threshold value can be obtained. For example, for the given earthquake magnitudes of 5.0, 6.0, and 7.5 and an uncertain magnitude, it can be seen from Figure 7 that the probabilities that the maximum flow capacities are smaller than 3310 (i.e., vertical line in red), which is the 75% level of the original maximum capacity, are estimated to be 1.4 × 10^{−4}, 6.1 × 10^{−2}, 9.9 × 10^{−1}, and 3.0 × 10^{−2}, respectively.

To verify the analysis results of the proposed method, an MCS with 3 × 10^{5} samples was conducted for *M* = 7.5. From the proposed method, the mean and standard deviation of the maximum flow capacity were calculated to be 471.45 and 742.97, respectively. Figure 8 shows how the mean and standard deviation of the maximum flow capacity obtained from the MCS converge to those from the proposed method as the number of samples increases. In the figure, it can be seen that the convergence of the mean and standard deviation is achieved with 2 × 10^{5} to 3 × 10^{5} samples, which takes approximately 12 minutes using a general-purpose personal computer with a 4.0 GHz CPU and 16 GB of RAM. However, it should be noted that the time cost is for only one set of earthquake magnitude and epicenter. To calculate the mean for an uncertain earthquake magnitude and location using the MCS, similar analyses should be conducted for many different sets of earthquake magnitudes and epicenters. In this numerical example, the range of earthquake magnitude (i.e., 4.5 to 7.5) is divided by 0.1 intervals and there are 20 epicenters. This means that the same MCS analysis should be conducted a total of 620 times, which will take approximately 124 hours. However, for the proposed method, it takes approximately 2 hours to obtain all of the analysis results, and this is possible because the matrix-based framework of the proposed method separates the system probability calculation and maximum flow capacity analysis, which enables a repeated analysis of the maximum flow capacity to be avoided.

**(a)**

**(b)**

For earthquake magnitudes between 4.5 and 7.5, the statistical moments of the maximum flow capacity according to the twenty earthquake sources are also obtained from the proposed method. Figure 9 shows the mean of the maximum flow capacity with varying earthquake magnitudes for the twenty earthquake locations presented in Table 5. From the results, it can be seen that the mean flow capacity decreases with increasing earthquake magnitude. When the earthquake magnitudes are relatively small, the mean of the maximum flow capacity is close to the original maximum flow capacity of the network, 4440, but it decreases with increasing earthquake magnitude. In addition, the rates of decrease are different in the twenty earthquake sources, since the associated site-source distances and focal depths are different.

**(a)**

**(b)**

**(c)**

**(d)**

For a better comparison, Figure 10 presents the mean flow capacities of the twenty earthquake scenarios and their average (colored blue) which means the mean of the maximum flow capacity for uncertain earthquake location (*L*). It can be clearly seen that the mean of the maximum flow capacity for uncertain location also decreases as the earthquake magnitude increases. In addition, it is found that, among the twenty earthquake scenarios, Earthquake 8 is the most critical one followed by Earthquake 20.

Similarly, the standard deviation and the c.o.v. of the flow capacities of the twenty earthquake scenarios and their averages (colored blue) for uncertain earthquake location can be calculated, which are presented in Figures 11 and 12, respectively. Overall, the standard deviation increases as the earthquake magnitude increases. In some scenarios, however, it decreases after certain magnitude. In addition, it is observed in Figure 12 that the c.o.v., i.e., a standardized measure of dispersion, increases with the increasing magnitude, which means the stronger a seismic event is, the more uncertainty the maximum flow capacity of the network has.

Furthermore, the mean, standard deviation, and c.o.v. of the maximum flow capacity of the network for an uncertain earthquake (i.e., uncertain magnitudes and locations of earthquake) can be calculated using equation (15), which are shown in Table 6.

It is noteworthy that sampling-based approaches would be inefficient for this sort of parametric study because network flow capacity analysis should be conducted for all of the individual magnitude values and locations of earthquake. On the contrary, the proposed method enables to perform this parametric study efficiently.

The proposed method also enables the computation of the reduction factor RF using equation (16). RF contains the performance measure; hence, the results can be used to investigate the relative importance of bridges in the network. Figure 13 shows RFs of all bridges for the two serious earthquake scenarios (i.e., Earthquake 8 and Earthquake 20) and for uncertain earthquake. As a result, it is observed that Bridges 3, 5, 6, 9, and 10 are relatively important in the region with the assumed origin and destination of Node 30 and Node 3, respectively.

Furthermore, one important advantage of the proposed method is its efficiency, which makes it feasible to calculate the sensitivity of the statistical moments with respect to each parameter of the bridge fragility curves. Although sampling-based approaches such as the MCS may not be feasible for such a sensitivity analysis, the proposed method enables a sensitivity calculation, by applying the method with the finite difference method. A fragility curve is determined by a median and a logstandard deviation [25], and Table 7 shows the sensitivity of the mean of the maximum flow capacity of the network for an uncertain earthquake with respect to the median of each damage fragility curve. Among the ten bridges considered in this numerical example, Bridges 3 and 4, which are found to be important and less important in Figure 13, respectively, are selected for a comparison purpose. Each sensitivity value in the table indicates how the mean of the maximum flow capacity of the network will increase when the median of each damage fragility curve of each bridge increases. The increase in the median of a bridge fragility curve means the bridge becomes more robust, and it can be seen in the table that the sensitivities of Bridge 3 are larger than those of Bridge 4. This shows Bridge 3 is more important than Bridge 4, which also means that it is more efficient to reinforce Bridge 3.

#### 5. Conclusions

In this study, a new method has been proposed for system-level seismic risk assessment of bridge transportation networks considering earthquake uncertainty. The proposed method consists of three steps: (1) component failure probability calculation of bridges based on PSHA; (2) system-level performance estimation of the transportation network using the matrix-based framework of the MSR method; and (3) seismic risk assessment based on the total probability theorem. In the proposed method, PSHA enables the seismic fragility estimation of the component bridges considering the uncertainty of earthquake locations and magnitudes, and it is systemically used to carry out the estimation of the postearthquake performance (i.e., maximum flow capacity) of the target bridge network by employing the matrix-based framework. The matrix-based framework enables efficient evaluations of the network performance with various magnitudes and locations of earthquake, without performing deterministic analyses of maximum flow capacity repeatedly.

The proposed method has been tested through its application to a numerical example of an actual transportation network in Korea, considering various earthquake scenarios with twenty locations and a range of earthquake magnitude from 4.5 to 7.5, and ten bridges with five damage states. As a result, the probabilistic distributions and statistical moments of the maximum flow capacity are obtained. It was observed that as the earthquake magnitude increases, the mean value of the maximum flow capacity of the network decreases but the uncertainty of the maximum flow capacity represented by the c.o.v. increases. In addition, the mean, standard deviation, and c.o.v. of the maximum flow capacity of the network for an uncertain earthquake are computed. The proposed method also enables the computation of the reduction factors for all bridges, and it was observed that Bridges 3, 5, 6, 9, and 10 are relatively important in the target transportation network. The proposed method can also provide the mean sensitivity of the maximum flow capacity of the network with respect to the median of each fragility curve. Based on the analysis results, it has been successfully demonstrated that the proposed method enables the system-level seismic risk assessment of bridge transportation networks considering the earthquake uncertainty.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This research was supported by a grant (19SCIP-B146946-02) from the Construction Technology Research Program funded by the Ministry of Land, Infrastructure and Transport of Korean government. This work was also supported by the 2019 Research Fund (1.190011.01) of UNIST (Ulsan National Institute of Science and Technology).