Transportation Resilience Modeling and Bridge Reconstruction Planning Based on Time-Evolving Travel Demand during Post-Earthquake Recovery Period

Abstract

After major earthquakes, communities may experience time-evolving population in terms of size and distribution, and varying travel demands, along with the displacement and recovery of residents caused by the damage and restoration of dwelling units. Community transportation can be significantly affected if the changes in population size and distribution are considerable. As a result, the post-earthquake infrastructure reconstruction process is essentially like urban replanning to meet the realistic traffic needs of the remaining and recovering residents and further maximize the sustainability of the community. To fill the gap in existing studies that considered the travel demand as fixed during the long-term recovery stage, it is important to investigate the effects of time-evolving travel demand on transportation resilience modeling and bridge reconstruction planning during the post-earthquake recovery period. A new methodology is proposed to analyze such impact by assessing the time-dependent resilience performance of transportation networks during the post-earthquake recovery stage. Traffic efficiency and safety are the two resilience performance indicators used to evaluate the transportation network. Post-earthquake infrastructure restoration planning is conducted using a heuristic algorithm based on the time-dependent resilience performance indicator. A demonstrative case study is carried out at Shelby County, Tennessee.

1. Introduction

1.1. Characteristics of the Post-Earthquake Recovery Stage

Earthquakes may have significant and long-term impacts on the population and infrastructure of modern communities, both of which substantially affect the performance of the transportation networks. Major earthquakes often cause changes in population size and distribution in the communities, which may affect the travel demand and origin-destination (OD) conditions accordingly. For example, after an earthquake, some residents may have to be relocated to other cities because of the damaged dwelling units and the lack of necessary living resources (Hazus 2020) [1]. These people may move back to the community along with the recovery of the lodging units, water, power, etc. The change in the population size can significantly affect the total travel demand of the community (Chang et al., 2012) [2]. In addition to those residents relocated to other cities, some residents may be displaced to the short-term or long-term shelters located at different zones in the community and move back to their pre-earthquake zones when their dwelling units are repaired. As a result, the community population is redistributed extensively during the post-earthquake recovery stage, which may considerably affect the traffic OD conditions of the community. So, it is important to incorporate the time-evolving feature of the travel demand when simulating the resilience performance of transportation networks during the post-earthquake recovery stage.

However, most of the existing studies assumed the travel demand and the OD condition during the post-earthquake recovery stage to be a fixed ratio of the pre-earthquake level (e.g., Alipour and Shafei, 2016; Wu et al., 2021) [3,4]. Such an assumption may be acceptable for minor hazards with no or slight effect on a community’s population displacement and redistribution. It, however, cannot reflect the reality when simulating transportation resilience during the recovery stage after disasters like major earthquakes. If the resilience performance simulation is conducted based on the fixed travel demand, the needs of the pre-disaster population, instead of the remaining and future residents during the recovery stage, will become the focus, which may deviate significantly from the reality. Thus, a more rational simulation of the performance and resilience of transportation networks subjected to hazards requires the consideration of the time-evolving population and travel demand of the community during the post-disaster recovery period.

The post-earthquake reconstruction process of damaged infrastructures can be deemed as an urban replanning if there are significant changes in population and travel demand (Di Ludovico et al., 2020) [5]. Earthquakes can cause severe damage to transportation infrastructures and considerably affect the traffic efficiency and travel safety of the community during the post-earthquake recovery stage (Wu et al., 2021) [4]. To make it even worse, the repair of the damaged transportation infrastructures may take up to several years to finish, making the transportation network in the community remain disrupted over an extended period during the post-earthquake recovery stage. Timely recovery of a disrupted transportation network is of utmost importance because it can not only help improve the transportation performance of the affected community but also, more importantly, expedite the long-term recovery efforts of other critical infrastructures that are heavily dependent on accessible, efficient, and safe transportation. Like urban planning, the reconstruction planning of damaged transportation infrastructures should focus on meeting the realistic needs of the current/remaining and the future/recovering residents. After major earthquakes, the time-evolving population size and distribution of the affected community resulting from the resident displacement and recovery may be at the similar magnitude of the population change of a fast-developing city. So, the reconstruction planning cannot be conducted based on a fixed travel demand like what has been done in most existing studies (e.g., Alipour and Shafei, 2016; Wu et al., 2021) [3,4]. Therefore, it is significant to develop rational methodologies to predict the time-dependent performance of transportation networks with time-evolving travel demand so that the stakeholders can make risk-informed decisions regarding optimal reconstruction planning and further improve the resilience of the whole community during the post-disaster recovery period.

1.2. Existing Studies on Post-Earthquake Transportation Resilience Modeling

Numerous studies have been conducted to investigate the performance and resilience modeling of transportation networks in the context of earthquakes. Some research efforts were focused on the connectivity of transportation networks or the accessibility to various resources (e.g., Aydin et al., 2018; Guo et al., 2015; Jenelius et al., 2006; Kilanitis and Sextos, 2019a; Kondo et al., 2012; Li and Zhou, 2020; Mahmassani et al., 2013; Qian et al., 2012; Viriyasitavat et al., 2011; Wu and Chen, 2019; and Wu and Chen, 2023a) [6,7,8,9,10,11,12,13,14,15,16]. Boakye et al. (2022) [17] introduced a novel framework to assess the influence of transportation networks on communities during natural disasters. The researchers measured both the well-being of the community and the equitable distribution of social justice. To evaluate the network’s ability to preserve the health of affected individuals, they established a performance indicator, which was derived from the concept of connectivity theory, and considered the shortest path between injured individuals and medical resources. Wu and Chen (2023a) [16] presented a methodology that comprehensively incorporated the redundancy of transportation systems during emergency response periods. They utilized the equivalent resistance theory and employed the emergency medical demand of various zones as the weighing factor. The researchers examined the effectiveness of private-vehicle-based emergency response and applied the developed metric to aid in pre-disaster mitigation planning. However, most of these studies were related to the transportation resilience modeling during the post-disaster emergency response stage when the traffic is moderate. The connectivity-based methodologies tend to be less sufficient for transportation resilience modeling during the post-earthquake recovery stage because of the considerable travel demand and the potential traffic congestion and delays.

Many studies have integrated traffic flow and travel time modeling methods, devising dependable approaches to simulate transportation resilience in scenarios with substantial traffic volumes (Chang et al., 2012; Feng et al., 2020; Wu and Chen, 2023b; Zou and Chen, 2019) [2,18,19,20]. Wu et al. (2021) [4] introduced a resilience indicator by considering both travel time and safety for assessing transportation networks during the extended recovery phase following an earthquake. The researchers examined the impact of partially closed links resulting from severe bridge damage on travel time and traffic safety. They employed static traffic assignment and the Poisson regression model to derive the travel time and traffic accident frequency. The study assumed that the travel demand during the long-term recovery stage is same as the pre-earthquake level, which is apparently different from reality. Some researchers suggested deriving the travel demand during the post-earthquake recovery stage by applying modification factors based on seismic intensity, building damage, or network capacity (Kiremidjian et al., 2007; Shinozuka et al., 2008; and Zhou et al., 2010) [21,22,23]. These travel demand modification methods have been followed by some studies when simulating transportation resilience. Alipour and Shafei (2016) [3] examined the influence of deteriorating components on network resilience. They developed a framework that considered the degradation of components over time and its consequential impact on the network functionality. By applying this framework to the transportation network in the San Francisco Bay Area, they utilized historical seismic activity data along with current estimates of component deterioration. They used reduction ratios for different trip purposes as functions of peak ground accelerations.

Chang et al. (2012) [2] introduced an in-depth exploration of modeling approaches for assessing transportation networks in the aftermath of earthquakes. The authors delved into various methodologies for damage assessment, network performance evaluation, and recovery planning. They introduced modification factors for trip production and attraction based on different functionalities and hazard scenarios of the zones. They also adopted the classic gravity model to estimate the changed travel demand during the post-earthquake stage. However, most of these studies were based on the modified but fixed travel demand of the network, failing to capture the time-evolving nature of the travel demand during the recovery stage. Some studies considered the time-evolving feature of travel demand during the post-earthquake recovery stage (Zhou et al., 2010; Kilanitis and Sextos, 2019b) [23,24]. Kilanitis and Sextos (2019b) [24] investigated the effects of earthquake-induced bridge damage and time-varying traffic demand on the resilience of road networks when simulating post-earthquake transportation network resilience. The study focused on the interaction between the network capacity and the trip generation during post-earthquake scenarios, based on which they introduced a trip reduction factor to simulate the time-evolving travel demand. However, post-earthquake traffic demand is affected not only by the network capacity, but also by the time-evolving population size and distribution.

1.3. Contributions Made by This Study

All existing studies summarized above contributed to the post-earthquake resilience modeling of transportation networks and reconstruction planning with some limitations. For example, most of these studies either used the pre-earthquake travel demand or applied a reduction or modification factor to quantify the number of trips during the post-earthquake recovery stage. These simplifications may differ from reality when simulating the performance of transportation networks, especially during the post-earthquake recovery stage after disasters. After major earthquakes, the time-evolving community population and travel demand may change drastically, which may differ significantly from the pre-earthquake community. Thus, transportation resilience modeling and reconstruction planning based on a specific ratio of the pre-earthquake travel demand may not reflect the unique characteristics of the transportation during the recovery stage. It would be desirable to incorporate the time-dependent OD demand into the simulation to rationally assess the transportation performance and resilience during the post-earthquake recovery stage. Moreover, the reconstruction planning of damaged transportation infrastructures was not properly conducted. As mentioned above, the reconstruction planning should focus on meeting the actual demands of the current/remaining and future/recovering residents like urban planning. However, the reconstruction plans in existing studies were based on approximated travel demands as a certain ratio of the pre-earthquake OD demand instead of using reasonably estimated time-evolving data.

This study is to investigate the time-dependent transportation resilience modeling and reconstruction planning during the post-earthquake recovery stage to close existing research gaps. The contributions of this study include (1) modeling more realistic time-dependent travel demand of the community during the post-earthquake recovery stage; (2) conducting the reconstruction planning with a focus on meeting the realistic demands of the current/remaining and future/recovering residents. Specifically, both the time-evolving population size and distribution are simulated based on earthquake-induced building damages, public shelter information, community demographics, and building recovery information. The time-evolving travel demand is further estimated based on the varying zone population, community demographics, and a verified travel demand model in the case study area.

2. Methodology

This study proposes a methodology to quantify the resilience performance of transportation networks during the post-earthquake recovery stage by incorporating the time-dependent travel demand of the community and support the reconstruction planning of damaged transportation infrastructure. As presented in Figure 1, the model begins with the given earthquake information (location, intensity, etc.) to generate the bridge and building damage scenarios using the seismic attenuation function (Atkinson and Boore, 1995) [25] and fragility curves of structures (Nielson and DesRoches, 2007) [26]. Building damages and community demographics (U.S. Census Bureau 2019) [27] are used to estimate the total displaced population (Hazus, 2020) [1] and the initial travel demand using the trip generation model (Kimley-Horn et al., 2007) [28]. Building recovery reflects the restoration progress of dwelling units, which is used to estimate the recovery of the relocated residents and the time-dependent population size and distribution, and travel demand (Hazus 2020) [1].

In Figure 1, bridge damage scenarios are related to link functionalities based on existing studies (Zhang et al., 2019) [29] and are generated using the Latin Hypercube sampling method. The time-dependent link travel time and traffic flow can be estimated using traffic assignment methods, the sampled link functionality scenarios, and the time-evolving travel demand (Wardrop, 1952) [30]. The traffic flow and traffic accident models based on the community’s historical traffic crash data can be used to derive the expected traffic accident frequency (Washington et al., 2010) [31]. The expected time-dependent resilience index is defined based on system travel time, traffic accident frequency, and bridge damage uncertainty. Finally, the optimal bridge reconstruction planning can be derived using heuristic algorithms based on the proposed resilience index. Figure 1 shows the detailed process of the methodology and

Table 1. List of notations.

Variables	Description
$AD{T}_i$	Average daily traffic per lane of road segment $i$
BPR	Bureau of Public Roads
$b_{res}$	Resource constraint for reconstruction planning
$c_i$	Traffic capacity of road $i$
$D{H}_i$	Number of total relocated households of zone $i$
$H{H}_i$	Total number of households for zone $i$
$H{H}_i^{t=0}$	Initial total number of households for zone $i$ at the beginning of recovery
$H{H}_i^t$	Number of households of zone $i$ at recovery time $t$
$L_i$	Length of link/road segment $i$
$m$	Weight parameter
$M$	The set of all sampled network damage scenarios
$me{d}_i$	Median peak ground acceleration value of damage state $i$
$\%M{F}_i$	Parameter related to multi-family dwelling units damage states
$\%MF{C}_i$	Possibility for multi-family residential buildings suffering complete damage
$\%MF{E}_i$	Possibility for multi-family residential buildings suffering extensive damage
$\%MF{M}_i$	Possibility for multi-family residential buildings suffering moderate damage
$MF{U}_i$	Total numbers of multi-family dwelling units in zone $i$
$M_w$	Moment magnitude
$n_b$	Total number of damaged bridges
$N{D}_i$	Number of displaced households move to the public shelter in zone $i$
$n_s$	Number of shelters in the community
$n_z$	Number of zones in the community
$\mathrm{PGA}$	Mean peak ground acceleration
$PO{P}_{total}^0$	Community population at the pre-earthquake stage
$PO{P}_{total}^t$	Community population at recovery time $t$
$q_i$	Traffic volume on road $i$
$\mathrm{R}$	Epicentral distance
$\%RD{U}_i^t$	Percentage of recovered dwelling units of zone $i$ at time $t$
$R{I}^t$	Time-dependent resilience index
$\%S{F}_i$	Parameter related to single-family dwelling units damage states
$\%SF{C}_i$	Possibility for single-family residential buildings suffering complete damage
$\%SF{E}_i$	Possibility for single-family residential buildings suffering extensive damage
$\%SF{M}_i$	Possibility for single-family residential buildings suffering moderate damage
$SF{U}_i$	Total numbers of single-family dwelling units in zone $i$
$ST{P}_{ij}$	Number of people from zone $j$ that are assigned to the public shelters in zone $i$
$ST{P}_j$	Number of people moving from zone $j$ to public shelters
T	Total recovery time of the system
$TAF$	Expected total accident frequency of the system
$TA{F}_0$	Expected total accident frequency of the system at the pre-earthquake stage
$TA{F}_t$	Expected total accident frequency of the system at recovery time $t$
$t_i$	Travel time of the road $i$
$t_i^0$	Free-flow travel time of road $i$
$t_i^s$	Travel time of link $i$ in sampled damage scenario $s$
$T{P}_i^H$	Average household trip production of zone $i$
$T{P}_i^t$	Total trip production of zone $i$ at recovery time $t$
$TTT$	Total travel time of all links
$TT{T}_0$	Total travel time of all links at the pre-earthquake stage
$TT{T}_t$	Total travel time of all links at recovery time $t$
$W_{MFC}$	Default weight for completely damaged multi-family dwelling units
$W_{MFE}$	Default weight for extensively damaged multi-family dwelling units
$W_{MFM}$	Default weight for moderately damaged multi-family dwelling units
$W_{SFC}$	Default weight for completely damaged single-family dwelling units
$W_{SFE}$	Default weight for extensively damaged single-family dwelling units
$W_{SFM}$	Default weight for moderately damaged single-family dwelling units
$x_i$	Vector of explanatory variables for traffic accident regression
$x_i^t$	Binomial variable specifying if bridge $i$ is being repaired at time $t$
$\alpha$	BPR function parameter
$\beta$	Vector of estimable parameters for traffic accident regression
${\beta }_c$	BPR function parameter
${\epsilon }_i$	Error term for traffic accident regression
${\zeta }_i$	Dispersion of damage state $i$
${\lambda }_i$	Expected number of accidents per unit time on road $i$
${\lambda }_i^s$	The expected number of traffic accidents of link $i$ in sampled damage scenario $s$

outlines the list of notations used in this paper.

2.1. Seismic Vulnerability and Household Displacement

2.1.1. Seismic Intensity Measure and Structural Damage

As a popular earthquake intensity measure, the peak ground acceleration, PGA, of a given location can be estimated using the attenuation laws in existing studies based on the fault information such as moment magnitude ($M_w$) and epicentral distance ($R$) (Campbell 1981) [32]. Typical attenuation laws are expressed as:

$$PGA=f\left(M_w,R\right)$$

(1)

The seismic fragility of structures (buildings and bridges, etc.) can be defined as the probability of occurrence of a given infrastructure damage state under given intensities of earthquakes, e.g., PGA. The failure in terms of a limit state is defined as the excess of the limiting values of the structural performance indicators. Based on the failure severity level, the damage states for bridges and buildings are usually divided into five types: no damage, slight/minor damage, moderate damage, extensive damage, and complete damage (Nielson and DesRoches, 2007; Hazus, 2020) [1,26]. The damage states of buildings and bridges can be used to estimate building habitability and link functionality. The population size and distribution and the initial travel demand at the beginning of the recovery stage can be calculated based on the building damage states.

2.1.2. Household Displacement

The damaged dwelling units may result in an estimable number of displaced households (e.g., Wang and van de Lindt, 2021; Wang et al., 2021; Wang and van de Lindt, 2022) [33,34,35]. These households may need immediate shelters such as family or friends’ houses, rental homes, and public shelters. For the convenience of estimation, two types of displaced households are assumed in this study: those moving to temporary homes out of the community and those moving to public shelters within the community. Households moving out of the community are considered to relocate to and live in their new communities/cities, while households moving to shelters within the community are assumed to have the same travel habits as the residents in those zones where the public shelters locate. The public shelters may be in different zones from those before the disaster. For both types of displaced households, it is assumed they will move back when their pre-disaster dwelling units are repaired. For each zone, the number of total displaced households is estimated as (Hazus 2020) [1]:

$$D{H}_i=\left(SF{U}_i\times \%S{F}_i+MF{U}_i\times \%M{F}_i\right)\times \left(\frac{H{H}_i}{SF{U}_i+MF{U}_i}\right)$$

(2)

where for zone $i$, $D{H}_i$ is the number of total displaced households; $SF{U}_i$ is the total number of single-family dwelling units; $MF{U}_i$ is the total number of multi-family dwelling units; and $H{H}_i$ is the total number of households. All these data are included in the demographics of the community (U.S. Census Bureau 2019) [27]. $\%S{F}_i$ and $\%M{F}_i$ are parameters related to the building damage states (Hazus 2020) [1], which are defined as:

$$\%S{F}_i=W_{SFM}\times \%SF{M}_i+W_{SFE}\times \%SF{E}_i+W_{SFC}\times \%SF{C}_i$$

(3)

$$\%M{F}_i=W_{MFM}\times \%MF{M}_i+W_{MFE}\times \%MF{E}_i+W_{MFC}\times \%MF{C}_i$$

(4)

where $\%SF{M}_i$, $\%SF{E}_i$, and $\%SF{C}_i$ are the probabilities of single-family residential buildings suffering moderate, extensive, and complete damage states, respectively. Similarly, $\%MF{M}_i$, $\%MF{E}_i$, and $\%MF{C}_i$ are the probabilities for multi-family residential buildings suffering moderate, extensive, and complete damage states, respectively. $W_{SFM}$, $W_{SFE}$, $W_{SFC}$, $W_{MFM}$, $W_{MFE}$, and $W_{MFC}$ are the default weights for moderately, extensively, and completely damaged single- and multi-family dwelling units with values of 0, 0, 1, 0, 0.9, and 1, respectively.

According to Hazus (2020) [1], the number of households displaced to public shelters within the community for a given zone can be estimated based on its population size and the income, ethnicity, home ownership, and age distribution of the residents in the zone. Thus, the number of households displaced to other cities can be estimated by deducting the households displaced to public shelters from the total displaced households. The estimated household displacement information can be used to derive the initial population size and distribution of the community at the beginning of the post-earthquake recovery stage.

2.2. Time-Evolving Population and Travel Demand

2.2.1. Time-Evolving Population

To estimate the time-evolving population size and distribution of the community and further derive the travel demand, the reconstruction of damaged dwelling units needs to be introduced to simulate the recovery of the displaced residents. After an earthquake, restoring damaged buildings is part of the community’s recovery process. Given the distributions of the recovery time of buildings with different damage states, the population of each zone at any given recovery time $t$ can be estimated by assuming residents will move back immediately after their damaged dwelling units are restored. Based on Equation (2), the initial total households for zone $i$ after dislocation is:

$$H{H}_i^{t=0}=\{\begin{array}{c}H{H}_i-D{H}_i, Zone i has 0 public shelter\\ H{H}_i-D{H}_i+N{D}_i, Zone i has public shelter\left(s\right)\end{array}$$

(5)

$${{\displaystyle \sum }}_1^{n_z}ST{P}_j={{\displaystyle \sum }}_1^{n_s}N{D}_i$$

(6)

where $N{D}_i$ is the number of displaced households that move to the public shelter in zone $i$; $n_z$ is the number of zones in the community; $ST{P}_j$ is the number of people moving from zone $j$ to public shelters; and $n_s$ is the number of shelters in the community. Considering that distance can be a significant factor for the shelter choice of displaced residents, this study reasonably assumes that the residents who need public shelters will seek temporary housing at the closest public shelter until it is full. Based on this assumption, $N{D}_i$ is calculated using linear programming.

According to Hazus (2020) [1], the mean recovery time of single-family and multi-family dwelling units is listed in

Table 2. Building recovery times in days (Hazus 2020) [1].

Occupancy Class	Recovery Time
	Structural Damage State
	None	Slight	Moderate	Extensive	Complete
Single-family dwelling	0	5	120	360	720
Multi-family dwelling	0	10	120	480	960

. Following existing studies, this study assumes the recovery time of dwelling units follows a lognormal distribution with a covariance of 0.4 (Lin and Wang 2017) [36].

As assumed above, the displaced residents will return to their pre-earthquake dwelling units soon after the buildings being restored to meet the habitability criterion. The total households of zone $i$ at time recovery time $t$ can be estimated as:

$$H{H}_i^t=\{\begin{array}{c}H{H}_i^{t=0}+{\%RDU}_i^t\ast {DH}_i, Zone i has 0 public shelters\\ H{H}_i^{t=0}+{\%RDU}_i^t\ast {DH}_i-{\sum }_j{\%RDU}_i^t\ast {STP}_{ij}, Zone i has public shelter\left(s\right)\end{array}$$

(7)

where $H{H}_i^t$ is the number of households of zone $i$ at recovery time $t$; $\%RD{U}_i^t$ is the percentage of recovered dwelling units of zone $i$ at time $t$; and $ST{P}_{ij}$ means the number of people from zone $j$ that are assigned to the public shelters in zone $i$. $\%RD{U}_i^t$ can be derived based on the recovery time (Table 2) and the specified distribution information. The time-dependent number of households, $H{H}_i^t$, can be used to estimate the travel demands of the community at any time during the recovery stage.

2.2.2. Time-Dependent Travel Demand

The recovery stage can last up to several years after major earthquakes. Unlike the usually sudden, short, and intense post-earthquake emergency response stage that can easily disrupt almost all regular travel plans (Wu and Chen 2023b) [19], the relatively long-lasting recovery stage may allow the remaining residents to keep up with their daily activities despite the ongoing reconstruction of the damaged infrastructures (Chang et al., 2012) [2]. The Memphis Metropolitan Planning Organization (MPO) developed a travel demand model in the Long-Range Transportation Plan for the Memphis region. They used trip generation, trip attraction, and the trip distribution matrix based on the gravity model to derive the travel demand of the region, the model settings of which have been verified to be reliable (Kimley-Horn et al., 2007) [28].

According to the Memphis MPO travel demand model (Kimley-Horn et al., 2007) [28], there are nine typical trip purposes, including journey to work, home-based school, home-based university, home-based shopping, home-based social-recreational, home-based pickup/drop-off, home-based other, non-home-based work, and non-home-based non-work. For household trip production, the numbers of trips generated for different purposes are related to at least two of the following parameters: number of persons in a household, number of persons aged 0–17 in a household, number of workers in a household, and number of vehicles in a household (Kimley-Horn et al., 2007) [28]. Based on the Memphis MPO model, the total trip production of a given zone at any recovery time can be estimated as:

$$T{P}_i^t=T{P}_i^H\ast H{H}_i^t$$

(8)

where $T{P}_i^t$ is the total trip production of zone $i$ at recovery time $t$; and $T{P}_i^H$ is the average household trip production of zone $i$ derived using the Memphis MPO model (Kimley-Horn et al., 2007) [28] with the demographics of the zone (U.S. Census Bureau 2019) [27].

For the trip attraction of each zone, the number of attracted trips for each purpose is related to one or more of the following parameters of the zone: total employment, school enrollment, university enrollment, retail employment, service employment, total households, and office employment (Kimley-Horn et al., 2007) [28]. All the above-mentioned zone information is available in the demographics of the Memphis area (U.S. Census Bureau 2019) [27]. This study follows the Memphis MPO model for travel demand estimation. A commonly used method, the gravity model (Erlander and Stewart 1990) [37], is adopted for trip distribution. The factors used to determine the trip distribution matrix are travel distances and production-attraction dummies.

Based on the estimated daily travel demand, the origin-destination (OD) matrix will be derived during the study period according to the historical time-of-day travel data of the Memphis area. Similarly, a mode choice logit model is adopted to estimate the percentage of vehicle-based trips. The factors considered in the mode choice logit model are origin-destination highway distance, production-attraction zone types (e.g., CBD, urban, suburban, and rural), and household auto ownership. The derived time-dependent OD demand during the study period can be used with various traffic assignment methods (e.g., Wardrop 1952 [30]) to get the link travel time and traffic accident frequency of the community.

2.3. Travel Time and Safety for Roads with Work Zones

2.3.1. Modified Bureau of Public Roads Function

Damaged bridges may disrupt road links following earthquakes and create work zones. According to Padgett and DesRoches (2007) [38], bridges with extensive damage can recover 50% of their capacities within one week after the disaster. Considering that the recovery stage usually does not start immediately after earthquakes, the results from their research will be applied in this study. We focus on two-way two-lane roads (two lanes in each direction), which are common in the US. In this study, a two-way two-lane bridge is deemed fully closed if it is completely damaged; and deemed at 50% capacity if the bridge suffers extensive damage but one lane in each direction remains open. With extensive damage, since one lane in each direction on the bridge will be closed, the link will be degraded into a road essentially with “work zones” in both directions (Figure 2). As a result, there are three typical traffic scenarios for two-way double-lane roads: normal four-lane traffic, disrupted traffic with work zones on the bridge (Figure 2), and fully closed to traffic.

Significant differences exist between the travel patterns of normal roads and those with work zones (Hou and Chen, 2019) [39]. To consider the impact of work zones on travel time, a modified Bureau of Public Roads (BPR) function is adopted from an existing study (Zhang et al., 2019) [29] to estimate the travel time of roads $i$ in Figure 2. The travel time function for road $i$ can be expressed as:

$$t_i=t_i^0\left[1+\alpha {\left(\frac{q_i}{c_i}\right)}^{{\beta }_c}\right]$$

(9)

where $t_i$ is the travel time of the road $i$; $t_i^0$ is the free-flow travel time of the road; $q_i$ is the traffic volume on the road; $c_i$ is the traffic capacity of the road; and $\alpha$ and ${\beta }_c$ are BPR function parameters. Different from the default values of 0.15 and 4.0 for $\alpha$ and ${\beta }_c$ for normal roads, these parameters have values of 1.429 and 4.293, respectively, for road $i$ with the traffic pattern as specified in Figure 2 (Zhang et al., 2019) [29]. According to the US Highway Capacity Manual (TRB 2010) [40], the ideal capacity of one intact lane is usually assumed to be 2000 vehicles. Zhang et al. (2019) [29] have confirmed that Equation (9) works when the ratio of heavy vehicles is not higher than 10%. Since trucks and heavy vehicles only account for 4.2% of the highway vehicles (Sprung 2018) [41], the BPR function can be applied in this study.

2.3.2. Work Zone Traffic Accident Estimation

For the road scenario in Figure 2, the changed traffic pattern will affect traffic accident frequency (Garber and Zhao, 2002; La Torre et al., 2017; Pigman and Agent, 1990; Weng and Meng, 2011; and Yang et al., 2015) [42,43,44,45,46]. As an important topic that has been widely investigated, work zone traffic safety can be as important as traffic efficiency to travelers. The Poisson regression model has been used for traffic accident counts but with a significant limitation: the mean and standard deviation of accidents need to be close. To handle the over-dispersed accident data, negative binomial regression is often used as an alternative for discrete and nonnegative events. According to Washington et al. (2010) [31], for the negative binomial distribution, the expectation of the number of accidents on road $i$ is:

$${\lambda }_i=e^{\beta x_i+{\epsilon }_i}$$

(10)

where ${\lambda }_i$ is the Poisson parameter for road $i$, which is equal to the expected number of accidents per unit time on road $i$; $x_i$ is a vector of explanatory variables; $\beta$ is a vector of estimable parameters; and $EXP\left({\epsilon }_i\right)$ is a gamma-distributed error term (Washington et al., 2010) [31]. For work zone traffic accidents, the explanatory variables may include link length, light conditions (daytime or nighttime), adjusted traffic volume, posted speed limit, speed reduction to the normal speed limit, the number of open lanes, the number of lane closures, road type, the number of ramps, and the number of intersections within the work zone (Ozturk et al., 2013) [47]. It should be noted that the explanatory variables for normal roads are usually different from those of work zones.

2.4. Resilience Index and Reconstruction Planning

2.4.1. Time-Dependent Resilience Index

The expected total travel time (TTT) is one of the commonly used system performance metrics in transportation engineering to characterize the efficiency of a transportation network:

$$TTT=\frac{{{\displaystyle \sum }}_s{{\displaystyle \sum }}_i{t}_i^s}{\left|M\right|}$$

(11)

where $t_i^s$ is the travel time of link $i$ derived from Equation (9) and traffic assignment methods in sampled damage scenario $s$; and $M$ is the set of all sampled network damage scenarios. Similarly, the expected total accident frequency (TAF) of the system is used to quantify the traffic safety risk of the network, which is defined as:

$$TAF=\frac{{{\displaystyle \sum }}_s{{\displaystyle \sum }}_i{\lambda }_i^s}{\left|M\right|}$$

(12)

where ${\lambda }_i^s$ is the expected number of traffic accidents of link $i$ derived from Equation (10) in sampled damage scenario $s$. The integrated resilience index (RI) at time $t$ is proposed to consider both traffic efficiency and travel safety performances. Considering that travel demand and community population are both functions of the recovery time, a time-dependent integrated resilience index of the transportation system to offset the effect of population change can be defined as:

$$R{I}^t=m·\frac{TT{T}_0}{PO{P}_{total}^0{}^{\delta }}/\frac{TT{T}_t}{PO{P}_{total}^t{}^{\delta }}+\left(1-m\right)·\frac{TA{F}_0}{PO{P}_{total}^0{}^{\delta }}/\frac{TA{F}_t}{PO{P}_{total}^t{}^{\delta }}$$

(13)

where $m$ is a weight parameter between [0, 1] usually defined by decision-makers based on the specific scenarios and preferences; $TT{T}_0$, $TA{F}_0$, $PO{P}_{total}^0,$ and $TT{T}_t$, $TA{F}_t$, $PO{P}_{total}^t$ are the TTT, TAF and community populations at the pre-disaster period and recovery time $t$respectively; and $\delta$ is an adjustment factor for RI. Considering that both TTT and TAF may change exponentially with recovery time, a proper $\delta$ will keep RI between 0 to 1. For demonstrative purposes, this study let $m=0.5$ to evenly distribute the importance to travel time and traffic safety and use $\delta =2$.

2.4.2. Reconstruction Planning of Transportation Infrastructure

According to the above-mentioned characteristics of the post-earthquake recovery stage, the post-earthquake reconstruction planning can be deemed as an urban replanning if the population size and distribution changes are significant. Compared to the reconstruction planning of the damaged infrastructures with a fixed travel demand that is usually derived based on the pre-earthquake population in existing studies, the infrastructure restoration in this study is based on the time-evolving travel demand of the community. Therefore, the reconstruction planning in this study focuses on serving the remaining and recovering residents instead of the pre-earthquake residents in most existing studies.

Considering the common situation of limited resources that would impede the repairing efforts of all damaged infrastructures simultaneously, this study focuses on developing the optimal reconstruction planning of the damaged bridges to improve transportation resilience. For long-term recovery, minor damages are not critical and there are two types of bridge damage scenarios according to Section 2.3: a completely damaged bridge that is fully closed to traffic, and an extensively damaged bridge that partially opens to traffic with work zones. Other than the bridge damage state, the reconstruction time of a bridge is also an important factor affecting its repair priority and effect on system performance. In this study, the performance indicator is based on the proposed time-dependent resilience index in Section 2.4.1. The optimal bridge reconstruction plan is to maximize the system resilience. The objective function of the reconstruction optimization is defined as:

$$\mathrm{Objective}Z=\max \left({{\displaystyle \int }}_0^{T}R{I}^{t}dt\right)$$

(14)

$$\mathrm{Subject}\mathrm{to}:{{\displaystyle \sum }}_{i=1}^{n_b}{x}_i^t\le b_{res}$$

(15)

where T is the total recovery time of the system; $x_i^t$ is a binomial variable specifying whether bridge $i$ is being repaired at time $t$; $n_b$ is the total number of damaged bridges; and $b_{res}$ is the resource constraint.

3. Case Study

3.1. Shelby County, Tennessee

Shelby County, Tennessee is in an earthquake-prone area. There are 37 zip codes in Shelby County and the demographics of each zip code area are available from the US Census Bureau (Figure 3 and Figure 4) (U.S. Census Bureau 2019) [27]. Among all the zip codes, 38054, 38058, 38131, and 38132 are not considered in this study because of their special functionalities (e.g., PO box, local airport). To better fit Shelby County in the model with a reasonable scope, every main road network intersection is considered as a single node/zone (Figure 5). The basic information for the zip codes and zones is presented in

Table 3. Zip code information of Shelby County.

Zip Code	TDP 1	Avg 2	Households
38002	25,802	1.84	14,053
38016	33,445	1.78	18,801
38017	48,012	2.67	17,980
38018	29,209	2.04	14,317
38028	4696	1.77	2646
38053	23,981	2.29	10,490
38103	35,566	5.08	7005
38104	37,758	3.02	12,517
38105	17,336	6.48	2676
38106	28,019	2.69	10,429
38107	15,159	2.14	7073
38108	19,164	2.85	6729
38109	43,959	2.59	16,989
38111	38,634	2.16	17,860
38112	18,764	2.75	6826
38114	25,507	2.4	10,624
38115	37,249	2.32	16,053
38116	43,802	2.96	14,793
38117	37,887	3.15	12,029
38118	101,321	7.21	14,062
38119	28,514	2.83	10,065
38120	42,874	6.28	6827
38122	23,456	2.34	10,044
38125	33,446	2.27	14,722
38126	10,891	4.18	2607
38127	35,673	2.33	15,320
38128	35,595	2.21	16,080
38133	35,919	4.67	7689
38134	42,981	2.86	15,053
38135	18,557	1.66	11,170
38138	34,024	3.37	10,082
38139	10,335	1.84	5629
38141	19,847	2.58	7678

Note: ¹ TDP means total daily time population, including daytime workers and daytime residents; ² Avg means average daily household size.

and

Table 4. Zonal information of Shelby County.

Zones	Longitude	Latitude	Households	Zone Type
1	35.343022	−89.788899	3513	Rural
2	35.308079	−89.822665	2098	Rural
3	35.288887	−89.686639	3513	Rural
4	35.2561	−89.676892	3513	Rural
5	35.210145	−89.653792	3513	Rural
6	35.158731	−89.694493	1323	Rural
7	35.195068	−89.792187	9401	Suburban
8	35.20436	−89.77236	2563	Suburban
9	35.156498	−89.774194	9401	Suburban
10	35.156446	−89.638166	3596	Rural
11	35.044817	−89.637635	3596	Rural
12	35.046272	−89.689285	3596	Rural
13	35.033613	−89.689272	3596	Rural
14	35.020894	−89.689143	3596	Rural
15	35.157881	−89.794974	14,317	Suburban
16	35.211229	−89.63495	1323	Rural
17	35.332123	−90.015529	2098	Rural
18	35.34534	−89.906004	2098	Rural
19	35.331266	−89.920981	2098	Rural
20	35.319164	−89.873442	2098	Rural
21	35.153188	−90.059146	7005	Suburban
22	35.137466	−90.024084	6259	CBD
23	35.130376	−90.023955	6259	Urban
24	35.155747	−90.037473	3537	CBD
25	35.151606	−90.040735	2676	CBD
26	35.124652	−90.054897	2086	Suburban
27	35.10264	−90.029534	2086	Suburban
28	35.103061	−90.019985	2086	Suburban
29	35.071525	−90.023332	2086	Suburban
30	35.158338	−90.020242	3537	CBD
31	35.182929	−89.932609	6729	Suburban
32	35.07939	−90.057321	8495	Rural
33	35.020647	−90.08822	8495	Rural
34	35.124294	−90.067707	2086	Suburban
35	35.131173	−89.956728	5953	Urban
36	35.133559	−89.956127	5953	Urban
37	35.135243	−89.968143	1707	Urban
38	35.079982	−89.930549	5953	Urban
39	35.080578	−89.955869	5312	Urban
40	35.133766	−89.984064	1707	Urban
41	35.150891	−89.98136	1707	Urban
42	35.158259	−89.980416	1707	Urban
43	35.108035	−89.98827	5312	Urban
44	35.057674	−89.833646	16,053	Urban
45	35.021095	−90.025306	7397	Suburban
46	35.020989	−90.002561	7397	Suburban
47	35.131269	−89.901452	6015	Suburban
48	35.11449	−89.901795	6015	Suburban
49	35.057248	−89.932866	4687	Suburban
50	35.021407	−89.937372	4687	Suburban
51	35.020528	−89.898062	4687	Suburban
52	35.052222	−89.812832	5033	Suburban
53	35.084536	−89.878106	5033	Suburban
54	35.101952	−89.870081	3414	Suburban
55	35.130947	−89.867721	3414	Suburban
56	35.171825	−89.943552	3348	Urban
57	35.149089	−89.927158	3348	Urban
58	35.152387	−89.884629	3348	Urban
59	35.037425	−89.795708	7361	Suburban
60	35.020556	−89.795966	7361	Suburban
61	35.142383	−90.051141	2607	CBD
62	35.236202	−90.014406	5107	Rural
63	35.199458	−90.034662	5107	Rural
64	35.191005	−90.016037	5107	Rural
65	35.235393	−89.895316	4020	Suburban
66	35.205276	−89.923297	4020	Suburban
67	35.192791	−89.926172	4020	Suburban
68	35.205485	−89.834333	11,170	Suburban
69	35.237388	−89.793563	2563	Suburban
70	35.204362	−89.791932	2563	Suburban
71	35.204642	−89.867163	3763	Suburban
72	35.180442	−89.860168	3763	Suburban
73	35.168234	−89.85416	3763	Suburban
74	35.093018	−89.808541	10,082	Suburban
75	35.075809	−89.774767	5629	Suburban
76	34.994567	−89.84841	7678	Suburban
77	35.207011	−89.894558	4020	Suburban
78	35.169095	−89.90233	3763	Suburban

respectively.

The main road network of Shelby County comprises 78 nodes and 115 links. This study evenly distributed the population of each zip code to the nodes inside it. It is assumed there are eight different types of buildings in the community according to the building classification in Hazus (2020) [1]: 1–2 stories wood frame, other wood frame, 1–3 stories steel moment frame, 4–7 stories steel moment frame, 1–3 stories concrete moment frame, 4–7 stories concrete moment frame, 1–3 stories concrete shear wall, and 4–7 stories concrete shear wall. Due to the lack of detailed and reliable data, the percentages of different building types are randomly generated for each zone for demonstrative purposes. Among all the links, it is assumed that 31 of them have bridges on them (Figure 6). There are five different types of bridges in the network: multi-span continuous (MSC) concrete, slab, and steel girder bridges, and multi-span simply supported (MSSS) concrete girder and steel girder bridges. In this study, it is assumed that only bridge damage may cause the failure or degradation of a link. Of the 115 links, 31 are, therefore, vulnerable to earthquakes due to the existence of bridges (B1–B31) on them (Figure 6). The geographic information of the nodes and bridges is available on Google Maps and through ArcGIS database. The households of each zone are divided into single-family and multi-family households. The percentage of single-family households is 70%, according to the demographics of Shelby County (Census profile 2018) [48].

3.2. Data Analysis

In this study, an earthquake with an epicentral depth of 10 km is assumed to happen in the community being investigated. Atkinson and Boore (1995) [25] developed the seismic attenuation law for earthquakes with an epicentral depth of 10 km. It is applied in this study to derive the ground motion at different zones and bridges. According to Atkinson and Boore (1995) [25], the mean peak ground acceleration at a certain distance from the epicenter can be estimated as:

$$\log \left(\mathrm{PGA}\right)={\mathrm{c}}_1+{\mathrm{c}}_2\left({\mathrm{M}}_w-6\right)+{\mathrm{c}}_3{\left({\mathrm{M}}_w-6\right)}^2-\log \mathrm{R}-{\mathrm{c}}_4\mathrm{R}$$

(16)

where ${\mathrm{M}}_w$ is the moment magnitude of the earthquake; $\mathrm{R}$ is the epicentral distance (km) of the location; and ${\mathrm{c}}_1$, ${\mathrm{c}}_2$, ${\mathrm{c}}_3$, and${\mathrm{c}}_4$ are the parameters from regression analyses. According to the study by Atkinson and Boore (1995) [25], the following parameters are used: ${\mathrm{c}}_1=3.79$, ${\mathrm{c}}_2=0.298$, ${\mathrm{c}}_3=-0.0536$, and ${\mathrm{c}}_4=0.00135$. Attenuation laws have many uncertainties, including magnitude, epicentral distance, ground motion, and limited knowledge about earthquake phenomena, and it is crucial to consider these uncertainties when deriving ground motion intensity measurements. According to Adachi and Ellingwood (2007) [49], PGA is assumed with a coefficient of variance of 0.6 to reflect the uncertainties.

As mentioned in Section 2.1, the peak ground accelerations derived from Equation (16) can be used to evaluate the building and bridge damage. According to Nielson and DesRoches (2007) [26], the probability of a given infrastructure damage state is defined as:

$$P\left[Damage State i or greater|PGA\right]=\mathsf{\Phi }\left[\frac{\ln \left(PGA\right)-\ln \left(me{d}_i\right)}{{\zeta }_i}\right]$$

(17)

where $me{d}_i$ and ${\zeta }_i$ are the median PGA value and the dispersion of damage state $i$, respectively. For certain damage states of buildings and bridges, the $me{d}_i$ and ${\zeta }_i$ are shown in

Table 5. Seismic fragility parameters of five bridge types (Nielson and DesRoches 2007) [26].

Bridge Type	Median PGA Values (g)		${\mathit{\zeta }}_{\mathit{i}}$
	$\mathbf{Extensive}\mathit{m}\mathit{e}{\mathit{d}}_{\mathit{i}}$	$\mathbf{Complete}\mathit{m}\mathit{e}{\mathit{d}}_{\mathit{i}}$
MSC concrete	0.75	1.03	0.7
MSC slab	0.78	1.73	0.7
MSC steel	0.39	0.5	0.55
MSSS concrete	0.83	1.17	0.65
MSSS steel	0.56	0.82	0.5

and

Table 6. Median and dispersion values for seismic fragility of eight building types (Hazus 2020) [1].

Building Type	Stories	$\mathbf{Moderate}\mathit{m}\mathit{e}{\mathit{d}}_{\mathit{i}}$	$\mathbf{Extensive}\mathit{m}\mathit{e}{\mathit{d}}_{\mathit{i}}$	$\mathbf{Complete}\mathit{m}\mathit{e}{\mathit{d}}_{\mathit{i}}$	${\mathit{\zeta }}_{\mathit{i}}$
Wood frame	1–2	0.43	0.91	1.34	0.64
	all other	0.35	0.64	1.13	0.64
Steel moment frame	1–3	0.22	0.42	0.8	0.64
	4–7	0.21	0.44	0.82	0.64
Concrete moment frame	1–3	0.23	0.41	0.77	0.64
	4–7	0.21	0.49	0.89	0.64
Concrete shear wall	1–3	0.3	0.49	0.87	0.64
	4–7	0.26	0.55	1.02	0.64

. Moderate $me{d}_i$, extensive $me{d}_i$, and complete $me{d}_i$ in Table 5 and Table 6 are the mean PGA values for moderate, extensive, and complete damages, respectively. Per Hazus (2020) [1], typical buildings are designed based on four seismic standards: high-code, moderate-code, low-code, or not seismically designed (also called pre-code). Buildings designed based on codes with higher seismic standards are more robust under earthquakes. This study assumes all the buildings in Table 6 were designed based on the moderate-code seismic standards.

According to the historical traffic data of the Memphis/Shelby County area, the morning rush hour (7:00 a.m.–8:00 a.m.) has the heaviest traffic of the day (Kimley-Horn et al., 2007) [28]. The resilience assessment of this study is based on the traffic performance of the morning rush hour during the recovery stage. The percentages of different trip purposes as introduced in Section 2.2.2 of the morning rush hour over total daily trips are presented in

Table 7. Percentages of trips by purpose.

	Percentages of Trips by Purpose over Daily Trips				All Purposes
Time Period	Journey to Work	Home-Based School/Home Based University	Other Home-Based Purposes	Non-Home Based
7:00–8:00 a.m.	16.7%	23.6%	7%	3.8%	12.52%

(Kimley-Horn et al., 2007) [28]. In this study, due to insufficient information, the percentages of different trip purposes during the morning rush hour are assumed to be the same as the pre-earthquake period throughout the recovery stage.

This study assumes the earthquake has an epicenter latitude and longitude of (35.107757, −89.940678). For comparison purposes, the data analysis also explored earthquakes with four moment magnitudes: 6.0, 6.5, 7.0, and 7.25. Given the demographic data, earthquake information, building fragility, and the community population during the recovery stage can be estimated using Equations (2)–(7). The expected total auto-mode OD demand during the morning rush hour at any recovery time can be derived from the household information, community population, and zone information using the Memphis MPO model (Kimley-Horn et al., 2007) [28]. The expected community population and OD demand are shown in Figure 7 and Figure 8 as functions of recovery time. The figures show both the community population and OD demand increase with the recovery time, which conforms to the reality that the displaced residents keep moving back to their original zones due to the restoration of their pre-earthquake dwelling units. The percentage increments in the community population and OD demand are similar because the community OD demand is a rough product of the community household number and the average household trips, considering the latter may only vary slightly among different zones within the community. As shown in Figure 7 and Figure 8, the higher the earthquake magnitude, the smaller the community population and OD demand, and the faster they recover. The lower starting point is consistent with the common understanding that higher magnitude generally causes more damaged housing units, more displaced residents, and greater reduced travel demand. The faster increase is caused by the lognormal distribution of the building recovery time. In this case, most damaged houses will recover within the mean recovery time (Table 2). Figure 7 and Figure 8 also show the changes in the community population and OD demand can be around 10%, which may significantly affect the transportation performance of the community. The considerable difference underscores the necessity to consider the time-dependent population and OD demand when simulating transportation resilience during the post-earthquake recovery stage.

The displaced population was estimated based on the information of the buildings (Table 6) designed according to the moderate-code seismic standards. However, there are many countries and regions where some or even most dwelling units were designed based on low-code or even pre-code seismic standards (Bilham, 2010; Ilki and Celep, 2012; Wang, 2008) [50,51,52]. Under major earthquakes, the affected communities may experience extensive infrastructure damages, a significantly displaced population, substantially decreased travel demands, and a long-lasting recovery stage (Naddaf, 2023) [53].

Considering the lack of the specific work zone traffic accident data in the Memphis area, this study uses the work zone traffic accident model by Ozturk et al. (2013) [47] to model the traffic safety in work zones. Their model uses the negative binomial regression based on the observed data in the State of New Jersey. The traffic accident model of fully functioning roads adopted in this study is based on the observed data in Tennessee to estimate the travel safety of roads with moderate or lower damages (Kiattikomol, 2005) [54]. Based on the existing studies, the case study assumes the speed limits are 55 mph and 45 mph for fully functional roads and work zones, respectively (Mannering, 2007; Mekker et al., 2016) [55,56]. Considering the trivial effects of ramps and intersections on traffic accidents according to the studies summarized by Ozturk et al. (2013) [47], this case study assumes no ramps and intersections on the links. Since all roads are typical four-lane roads, the numbers of operating and closed lanes are both two for the road situation in Figure 2. Based on the available information, the traffic accident frequency of road $i$ is expressed as:

$${\lambda }_i=\{\begin{array}{c}2\times {10}^{-10}·L_i{}^{0.7856}AD{T}_i{}^{2.26}, Fully functional roads\\ 0.0799·L_i{}^{0.477}AD{T}_i{}^{0.512}, Road with work zones\end{array}$$

(18)

where $L_i$ is the length of link/road segment $i$; and $AD{T}_i$ is the average daily traffic per lane of road segment $i$. $AD{T}_i$ can be derived from the community travel demand and the traffic assignment method. The link lengths can be estimated from the longitudes and latitudes (Table 4) of the links’ starting and ending nodes.

3.3. Results

This study assumes that only bridge reconstruction will generate work zones during the post-earthquake recovery period. Assuming that the above-mentioned earthquake has a moment magnitude of 7.25, Latin hypercube sampling is used to simulate bridge damage scenarios with 10,000 samples. Based on the estimated travel demand (Figure 8), the traffic simulation for the morning rush hour is conducted during the post-earthquake recovery stage through static traffic assignment. The expected system total travel time, total accident frequency, and resilience index are calculated using Equations (11)–(13) and (18). They are presented in Figure 9, Figure 10 and Figure 11 as functions of recovery time.

The box plots in Figure 9 and Figure 10 suggest both the expected total travel time and total accident frequency increase with time. Such phenomena are caused by the increasing travel demand along with the recovery process of the population, and thus are in accordance with reality. The increased percentages of the expected TTT and TAF are 3.9% and 14.07%, respectively. The smaller increase in TTT may be caused by the fact that the traffic volumes on most links are lower than their capacities, which reflects a trivial increase in their travel time. Only a few links with high traffic volumes have significant traffic delays and contribute to the increase in TTT. The larger increase in TAF may result from the increased traffic volume. Unlike link travel time, traffic accident frequency is more sensitive to link traffic volume, whether it is below the link traffic capacity or not (Equation (18)). Figure 11 shows that the expected RI also increases with the recovery time, despite the increase in TTT and TAF. It is because the resilience index is related not only to the expected total travel time and accidents on the transportation network, but also to the community population it serves. As a result, the system resilience index of the transportation network gradually increases during the recovery period.

The pre-earthquake TTT and TAF are 1389 and 4305, respectively, assuming no link damage. Figure 9 and Figure 10 show that TTT and TAF for most sampled network damage scenarios fall in ranges of [1533, 1694] and [5464, 7160], respectively. The results suggest that TAF is more sensitive to bridge damages with larger variations as compared with TTT. The less-deviated TTT may be caused by the existence of many redundant links, so there are always alternative paths with similar travel time in case of link damage. The larger variations in TAF may result from partially functional bridges with work zones. This corroborates the traffic accident model in Equation (18), which shows that roads with work zones tend to have higher traffic accident frequencies than fully functional roads, unless the traffic is very heavy. This observation also proves the significance of considering work zone traffic accidents when simulating transportation resilience during the post-earthquake recovery stage, especially when the recovery takes a very long time. Although the expected TTT in this case study is relatively stable, the expected TAF presents a time-evolving difference of 14.07%, which is non-negligible. Like traffic efficiency, traffic safety is also critical to the travel experience of the remaining and recovering residents. The considerable difference in TAF throughout the recovery stage confirms the significance of simulating the time-evolving community population and travel demand. Also, as mentioned above, many countries and regions have many older buildings designed using codes with lower standards or even pre-code seismic standards (Bilham, 2010; Ilki and Celep, 2012; and Wang, 2008) [50,51,52]. For those areas, the impact of earthquakes on population, travel demand, and traffic patterns can be more significant, so that both TTT and TAF may exhibit considerable changes during the recovery stage.

3.4. Post-Earthquake Bridge Reconstruction Planning

Due to the usually limited resources, including monetary budget, personnel, material supply, etc., it is not realistic to reconstruct all the damaged bridges simultaneously. To investigate the effect of time-evolving travel demand and the restoration priorities of bridges with different damage states, this study assumes that only one bridge can be repaired at a time for demonstrative purposes, without losing generality. So, it is meaningful to repair the most critical bridges first. In this study, the most critical bridges will be identified based on their impacts on the system resilience index. Considering the computational complexity to include the bridge damage uncertainties in the optimization, a specific bridge damage scenario was assumed in

Table 8. Bridge repair time by damage state.

Bridge #	Damage State	Repair Time (Week)
3	Extensive	9
7	Extensive	7
8	Complete	17
9	Extensive	6
10	Complete	22
11	Extensive	9
13	Complete	18
14	Extensive	5
16	Extensive	7
18	Complete	22
24	Complete	16
25	Complete	14
26	Complete	14
28	Extensive	7

to demonstrate the reconstruction planning. All bridges to be repaired have suffered either extensive or complete damage; however, moderately or less-damaged bridges will fully recover their traffic-carrying capacities before the recovery stage (Padgett and DesRoches, 2007) [38]. According to Twumasi-Boakye and Sobanjo (2018) [57], the estimated mean recovery time of general highway bridges with extensive and complete damages is 48.5 and 133 days, respectively. Assuming the repair time obeys a normal distribution with a coefficient of variation of 0.2, the repair time of the damaged bridges in this study is sampled and presented in Table 8. The post-earthquake time-evolving travel demand after the moment magnitude 7.25 earthquake (Figure 8) is used to carry out the reconstruction planning optimization. The effect of an intuitive reconstruction sequence (simply following the order of bridge numbers) on the system resilience index is also investigated for comparison purposes. The goal here is to identify which plan is better in terms of meeting the demands of the remaining and recovering residents in the community.

Based on Equations (14) and (15), the bridge reconstruction optimization is derived using the Greedy Algorithm. The Greedy Algorithm solves an optimization problem by selecting the best option available at every step. For the bridge reconstruction optimization in this study, the repair of every bridge is deemed as a step. To achieve the optimization objective specified in Equation (14), the algorithm will repair the bridge that can increase the resilience index the most quickly at every single step. The optimal reconstruction sequence of the bridges is 9, 14, 3, 28, 7, 11, 16, 25, 10, 18, 24, 13, 26, and 8. For both the optimal and intuitive reconstruction plans, their corresponding resilience indexes as functions of recovery time are presented in Figure 12. Figure 12 shows that both curves start and end at the same values, which conforms to the reality that the two reconstruction plans have the same starting and ending bridge damage situations. The curve of the resilience index for the optimal plan increases fast at the beginning and slows down at $t=120 \mathrm{weeks}$. The curve of the resilience index for the intuitive plan rises slowly during the first 60 weeks, then increases similarly to the optimal plan, and then rises sharply at $t=173 \mathrm{weeks}$. The curve of the resilience index for the optimal plan has higher values than the intuitive plan throughout the recovery stage, which means the optimal reconstruction plan provides a better transportation experience to the remaining and recovery residents. This phenomenon is likely because the critical bridges have received priority in the optimal plan. It is also found that all seven bridges with “extensive” damage states are fixed first. Although bridges losing full functionality may have more significant impact on the capacity of the transportation network than those only losing partial functionality, extensively damaged bridges have larger impacts on the overall system resilience, since TAF is very sensitive to the existence of work zones. In contrast, TTT is more stable than TAF due to many redundant links according to the above analysis. Also, the shorter repair time of extensively damaged bridges means that their restoration can increase the resilience index faster.

4. Conclusions

Transportation networks may remain disrupted for years, significantly affecting the traffic efficiency and safety of the community’s remaining and recovering residents during the post-earthquake recovery stage. A framework to model the time-dependent resilience performance of transportation networks during the post-earthquake recovery stage was developed, with a focus on the transportation experience of the current and future residents of the attacked community. To simulate the uncertainties in bridge damage scenarios, the seismic attenuation law, bridge fragility functions, and Latin hypercube sampling method were adopted in this study to generate numerous network disruption scenarios. The building fragility functions, building recovery information, community demographics, and the MPO model were applied to estimate the time-evolving community population and travel demand. The modified BPR function for work zones, static traffic assignment method, and work zone traffic accident model based on negative binomial regression were applied to simulate the traffic efficiency and travel safety during the morning rush hour throughout the recovery stage.

This study defined a new time-dependent integrated resilience index to assess the transportation performance experienced by the time-evolving community residents. The methodology incorporates uncertainties into link degradation to make the resilience modeling more reliable. The proposed resilience index was used to derive the optimal bridge construction plan to maximize the overall benefits of the time-evolving community residents. Finally, the proposed methodology was demonstrated in Shelby County, Tennessee, which belongs to an earthquake-prone area. A comparative analysis was conducted in terms of reconstruction planning to investigate the effects of the optimal reconstruction plan on transportation resilience during the post-earthquake recovery stage.

It is found that earthquakes may significantly affect the community population and travel demand and, in turn, impact the time-dependent travel time, traffic accident, and system resilience index. Therefore, the model is more realistic and reliable than traditional methodologies that did not reasonably consider the time-evolving population and travel demand. The reconstruction planning emphasized the significance of focusing on the needs of adopting a more realistic population group instead of a fixed resident group derived from the pre-earthquake population. This study still has some limitations: due to the lack of reliable and site-specific data, some assumptions had to be made, including some specific parameters related to buildings and bridges, the bridge reconstruction time, the travel behavior of residents, etc.; the influence of sequential earthquakes on the infrastructure fragility curves was not considered, which can be significant to both the travel demand and bridge damage (Yaghmaei-Sabegh & Mahdipour-Moghanni, 2019) [58]; some factors that may affect the traffic simulation were not incorporated in this study, including traffic signal lights and the correlation between traffic delays and accidents; and the reconstruction planning was scenario-based due to the computational complexity when both the time-evolving feature and uncertainty are involved at the same time. Despite these assumptions and simplifications in the demonstrative study, the proposed methodology is general enough to be applied to conduct a more advanced study once more reliable and site-specific data, as well as more powerful computational technologies, become available in the future. Finally, although this methodology was proposed based on earthquakes, it can be easily adapted and applied to the long-term recovery stage involving other disasters or incidents causing similar disruptions and time-evolving population and travel demand.