Life-cycle seismic resilience of deteriorating highway bridges under mainshock-aftershock sequences

Abstract

Reinforced concrete highway bridges often experience seismic events coupled with progressive deterioration due to corrosion over their lifespan. In many cases, the seismic mainshocks are followed by significant aftershocks, causing cumulative damage that disrupts normal operations and delays the restoration of bridges. Current research however lacks detailed assessment of seismic resilience considering the effects of aftershocks and corrosion deterioration over the lifetime of bridges. This study proposes a methodology for evaluating the life-cycle seismic resilience of deteriorating structures under mainshock and aftershock (MS-AS) sequences. Three multi-span reinforced concrete highway bridges with different geometries are used as benchmarks. A suite of 80 pairs of ground motion sequences is selected for undertaking the resilience evaluations based on the seismic scenarios considered. The Park-Ang damage index is adopted for the purpose of quantifying the cumulative damage. Nonlinear dynamic analysis is used to provide detailed insights into the mechanisms through which the aftershocks affect the cumulative damage. Based on the results, time-dependent system fragility curves under MS-AS sequences are developed in conjunction with a cumulative damage capacity model for the bridge piers. The seismic resilience of the bridges is subsequently assessed under ground motion sequences at different service times, and the effects of aftershocks and corrosion-induced deterioration on the resilience are examined. Finally, the life-cycle seismic resilience of the deteriorating benchmark bridges under MS-AS sequences is evaluated using the suggested framework. It is shown that the influence of aftershocks on the cumulative damage depends on a number of inter-related factors, including the relative mainshock-aftershock intensity as well as the dynamic characteristics of the bridges. The findings highlight the merits of the proposed framework in evaluating the life-cycle seismic resilience of bridges for different hazard scenarios and deterioration conditions.

1 Introduction

Highway bridges are a crucial component of the transportation network, facilitating economic exchange and social development between regions. Past notable seismic events worldwide have shown that earthquakes typically occur in sequences, whereby a mainshock (MS) is generally followed by aftershocks (AS) (Goda 2012; Di Sarno 2013; Li et al. 2014; Raghunandan et al. 2015; Yu et al. 2021). Numerous post-earthquake reconnaissance reports have indicated that bridges are highly susceptible to severe damage during ground motion sequences (Priestley 1988; Ruiz-Garcia et al. 2008; Di Sarno et al. 2019; Ozkula et al. 2023). For instance, several reinforced concrete (RC) piers of the I-5/I-605 separator, a nine-span highway overpass, sustained damage during the October 1, 1987, Whittier Narrows earthquake (M_w = 5.9) and experienced further structural damage after the main aftershock (M_w = 5.3) three days later. During the August 2016 earthquake sequence in Central Italy (M_w = 6.1), several masonry bridges located at transportation hubs suffered severe damage to components such as the partial collapse of wing walls. Other bridges were also heavily damaged during the February 6, 2023, Turkey earthquake sequence (M_w = 7.8). Consequently, the seismic performance of highway bridges during mainshock-aftershock (MS-AS) sequences has attracted considerable attention from both practitioners and researchers, prompting extensive studies on the topic.

The assessment of seismic behavior of highway bridges under MS-AS sequences primarily includes deterministic (Ruiz-Garcia et al. 2008; Fakharifar et al. 2015; Shin and Kim 2017) and probabilistic analysis methods (Alessandri et al. 2013; Jeon et al. 2016; Omranian et al. 2018; Abdollahzadeh et al. 2021; Chen et al. 2022, 2024; Sedaghati et al. 2022), as well as experimental shake table studies (Di Sarno et al. 2017; Zhou et al. 2019). Using deterministic approaches, Ruiz-García et al. (2008) assessed the nonlinear response of typical Mexican RC highway bridges exposed to MS-AS sequences, and found that aftershocks could increase the lateral drift demands. Shin and Kim (2017) emphasized the significant influence of the frequency content of aftershocks on the seismic response of bridge piers by comparing the peak and residual displacements under spectrally-matched and real aftershock records. Fakharifar et al. (2015) compared the performance of three types of repair jackets on RC bridge columns damaged by a mainshock and subjected to aftershocks, concluding that the three repair jackets can effectively improve the collapse capacity of the bridge. With the development of the performance-based earthquake engineering (PBEE) framework (J Moehle and GG Deierlein 2004), researchers have incorporated the uncertainties of ground motions, structural parameters, and material properties to provide a more comprehensive probabilistic assessment of the performance of highway bridges under MS-AS sequences. For example, Alessandri et al. (2013) proposed an effective tool to expedite the decision-making process for assessing the seismic risk of bridges damaged by mainshocks, and subjected to aftershocks. Jeon et al. (2016) introduced a methodology within the PBEE framework to develop aftershock fragilities, enhancing post-earthquake evaluation of damaged bridges and ensuring the functionality of transportation networks. Chen et al. (2022) examined the seismic fragility and direct financial losses of bridges with tall piers under MS-AS sequences, using a vector-based intensity measure that accounts for both mainshocks and aftershocks. Their findings suggested that completely ignoring the impact of aftershocks could lead to underestimating the probability of failure and financial losses for such bridges, especially in cases where strong aftershocks follow mainshocks. Chen et al. (2024) established a framework for MS-AS fragility analysis using an input-output Hidden Markov Model, providing a novel approach and insights into seismic fragility under multiple aftershocks. Zhou et al. (2019) investigated the low-cycle fatigue performance of transverse steel dampers (TSDs) under ground motion sequences using a 1/35-scale model of a cable-stayed bridge in a shake table test. The test results indicated that the TSD exhibits excellent low-cycle fatigue performance, maintaining its original performance without any fatigue failure even under multiple strong earthquakes.

Although existing research has made progress in exploring the seismic response and risk assessment of bridges under MS-AS sequences, previous studies primarily used the peak response of bridge components as the engineering demand parameter (EDP). Given that MS and AS constitute a continuous loading process, where bridge structures may undergo multiple loading cycles resulting in cumulative damage, there is a need for further research on these cumulative damage effects under ground motion sequences. Within the few available studies, Jung and Andrawes (2018) evaluated the effectiveness of external shape memory alloy (SMA) spirals in retrofitting multiple-frame bridges during strong MS-AS sequences by capturing the cumulative damage of reinforcing bars, and concluding that SMA confinement helps reduce fatigue damage. Ghosh et al. (2015) also proposed a framework to predict damage accumulation in structures subjected to multiple shock scenarios by developing damage index prediction models accounting for the probabilistic nature of the hazard. It was shown that repeated shocks within the time window of interest significantly increased the damage index exceedance probabilities. However, the aforementioned studies did not delve in detail into the mechanisms by which aftershocks contribute to the cumulative damage of bridge piers.

Beyond assessing the seismic behavior of structures, understanding their functionality loss during disasters, and the associated post-disaster recovery costs, is essential for making informed decisions regarding the operation, maintenance, and management of civil infrastructure (Zhou et al.2024). To this end, resilience assessment has become a crucial tool, attracting significant attention from the academic and engineering communities (Bruneau et al. 2003; Bruneau and Reinhorn 2007; Cimellaro et al. 2010). Dong and Frangopol (2015) and Gidaris et al. (2022) independently developed probabilistic methods for assessing the seismic resilience of bridges under MS-AS sequences. However, RC bridges are subjected to the detrimental effects of aging processes throughout their service life due to aggressive chemical attacks, such as chlorides (Frangopol 2011; Biondini et al. 2015). This leads to corrosion of steel reinforcement and deterioration of concrete, which in turn reduces the capacity of the structure to withstand seismic loading and negatively impacts the resilience of the system (Titi and Biondini 2013). Other previous studies have demonstrated that corrosion-induced deterioration significantly affects the cumulative damage to highway bridges (Panchireddi and Ghosh 2019; Otárola et al. 2022), highlighting the necessity for evaluating the life-cycle seismic resilience of bridges from a cumulative damage perspective. However, there is currently a lack of research on the seismic resilience of deteriorating individual bridges over their service life under MS-AS sequences.

In addition to the above considerations, since earthquakes occur randomly throughout the service life of a bridge, it is necessary to incorporate the stochastic nature of hazard occurrence in life-cycle resilience assessments (Li et al. 2020b). Among the limited available studies, Yang et al. (2019) introduced a renewal theory to model lifetime hazards and proposed a general approach for the life-cycle management of civil infrastructure under natural and/or manmade hazards. On the other hand, Li et al. (2020b) developed an integrated framework for assessing the long-term resilience and losses of highway bridges exposed to multiple hazards and subsequently proposed a long-term loss higher-order analysis method based on the moment-generating function (Li et al. 2020d). However, the resilience loss due to deterioration from aggressive environments has not been incorporated in these studies. Moreover, existing studies have not considered the influence of aftershocks within the lifetime resilience assessment framework.

To address the aforementioned research needs, this study deals with the life-cycle seismic resilience of deteriorating highway bridges subjected to mainshock and aftershock (MS-AS) sequences. An integrated framework to quantify the long-term seismic resilience of deteriorating structures is developed. The seismic behavior of bridges is assessed using cumulative damage, with the Park-Ang damage index serving as the quantification metric. Three continuous multi-span RC highway girder bridges with different geometric configurations are adopted as the benchmark structures. Based on the hazard scenarios considered, 80 pairs of MS-AS sequences are selected. Nonlinear time history analyses are conducted on the benchmark bridges. Drawing from the nonlinear dynamic analysis results, the mechanisms through which aftershocks affect the cumulative damage of bridges are elucidated, and the influence of corrosion-induced deterioration on cumulative damage is examined. By integrating the capacity model of the structure, the time-dependent system fragility curves are established. Finally, the proposed framework is applied to estimate the life-cycle resilience of bridges under MS-AS sequences.

2 Framework for long-term seismic resilience of deteriorating structures

2.1 Assessment of seismic resilience

Seismic resilience refers to the ability of the system to mitigate hazards, manage disaster effects when they occur, and execute recovery efforts that minimize social disruption and reduce the consequences of future earthquakes (Bruneau et al. 2003). Figure 1(a) illustrates the concept of seismic resilience. As shown in the figure, once an earthquake occurs, the functionality (Q) rapidly declines to Q_r. During the recovery period (δ_r), various recovery measures are implemented to restore the functionality of the system to the desired level (Q_t). Here, δ_r = τ_r - τ_i, where τ_i is the start time of recovery activities, calculated as τ₀ + δ_i. In this equation, τ₀ signifies the onset of the earthquake, and δ_i is the time interval between the occurrence of the earthquake and the start of recovery activities, while τ_r is the completion time of the recovery process.

Fig. 1

Illustration of seismic resilience quantification: a concept of seismic resilience, b recovery functions, and c fragility curves used to determine resilience

Based on this definition, a commonly used expression for quantitatively calculating the resilience (R) (Frangopol and Bocchini 2011) is adopted in this study, as follows:

$$R=\frac{1}{{{\tau _h} - {\tau _0}}}\int_{{{\tau _0}}}^{{{\tau _h}}} {Q\left( \tau \right)} d\tau$$

(1)

where $\:{\tau}_{\text{h}}$ represents the time horizon under investigation, and Q($\:\tau$) denotes the functionality level of the considered infrastructure. The rapidity of the post-earthquake recovery process of a structure is related to its damage state (DS). According to Biondini et al. (2015), the functionality recovery profile $\:{{Q}}^{\text{}({k})}\left(\tau\right)$ for a given damage state k can be effectively described by the following mathematical expression:    

$${Q^{\left( j \right)}}\left( \tau \right)={Q_r}+H\left( {\tau - {\tau _0} - {\delta _i}} \right) \cdot {R_f}\left( {\frac{{\tau - {\tau _0} - {\delta _i}}}{{{\delta _r}}}} \right) \cdot \left( {{Q_t} - {Q_r}} \right)$$

(2)

where j is the indicator of the damage state, and j = 1, 2, 3, and 4 correspond to slight, moderate, extensive, and complete damage state, respectively. H(·) denotes the Heaviside step function, and R_f(·) represents the recovery functions. Different R_f(·) functions can reflect varying speeds of recovery, corresponding to different damage states of the structure after the earthquake. In this study, three types of R_f(·) functions are adopted as shown in Fig. 1(b), and can be expressed as follows:

$$\begin{gathered} R_{f}^{n} (\eta ) = 1 - e^{{ - \omega \eta }} \hfill \\ R_{f}^{s} (\eta ) = \frac{{1 - \cos (\pi \eta )}}{2} \hfill \\ R_{f}^{p} (\eta ) = e^{{ - \omega (1 - \eta )}} \hfill \\ \end{gathered}$$

(3)

where ω is a shape parameter with a value of 10 (Biondini et al. 2015; Pang and Wang 2021a; Zhou et al. 2024), and η = $\:\left({\tau}-{\tau}_{i}\right)/{\delta }_{r}\in\left[\text{0,1}\right]\:$ represents the normalized time variable. The negative-exponential-type function $\:{{R}}_{{f}}^{{n}}$(η) models recovery processes where the majority of functionality is rapidly restored following the seismic event. It is applied to the recovery process of a slightly damaged state. The function $\:{{R}}_{{f}}^{{s}}$(η), a sinusoidal-type, depicts the gradual restoration of functionality over time and applies to the moderate damage state. Lastly, the positive-exponential-type function $\:{{R}}_{{f}}^{{p}}$(η) effectively models recovery processes where functionality is primarily restored towards the end of the recovery period, applicable to both extensive and complete damage states.

Due to the complex and stochastic nature of structural post-earthquake recovery in engineering practice, uncertainties in parameters are considered during the development of the recovery profile. Table 1 gives the probabilistic distributions of key parameters in the functionality recovery profile across different damage states. Using the provided probability distributions, key parameter samples of the functionality recovery curve corresponding to the specified damage state are obtained through 10⁴ Monte Carlo simulations (MCS). In each simulation, every set of sampled key parameters, combined with the recovery function R_f(·) from Eq. (3), is applied in Eq. (2) to calculate the bridge functionality $\:{{Q}}^{({k})}\left({\tau}\right)$ at specified time instants for a specified damage state. Ultimately, the expected functionality curve for each damage state is generated by averaging the bridge functionality from 10⁴ MCS at every time instant.

Table 1 Distribution and values of random variables associated with the recovery evolution of functionality across different damage states (Decò et al. 2013; Pang and Wang 2021a)

The functionality of the structure under a given seismic intensity measure (IM) level, Q($\:{\tau}\text{|}{IM}$), can be estimated using Eq. (4) (Decò et al. 2013; Pang and Wang 2021a):

$$Q\left( {\tau \left| {IM} \right.} \right)=P_{f}^{0}\left( {IM} \right) \cdot 1+\sum\limits_{{k=1}}^{4} {P_{f}^{{\left( j \right)}}} \left( {IM} \right) \cdot {Q^{\left( j \right)}}\left( \tau \right)$$

(4)

where Q^(j)($\:{\tau}$) is the expected functionality considering uncertainties in the recovery process (see Eq. (2), $\:{{P}}_{{f}}^{\text{0}}\left({IM}\right)$ represents the probability of remaining in the undamaged state following the earthquake, while$\:\:{{P}}_{{f}}^{\text{(}{j}\text{)}}\left({IM}\right)$ indicates the probability of being in the j^th damage state (where j = 1,2,3, and 4). Both probabilities can be derived from the fragility curves depicted in Fig. 1(c). The mathematical expression for obtaining the probability of damage (i.e., fragility) typically follows the following general form (Nielson and DesRoches 2007):      

$${\text{Fragility = }}P\left[ {D \geqslant C\left| {IM} \right.} \right]=P\left[ {C - D \leqslant 0\left| {IM} \right.} \right]$$

(5)

where D represents the seismic demand on the structure, and C denotes the structural capacity. In this study, system fragility is employed to derive functionality and resilience. It is assumed that when any component reaches a specific damage state (Fail_comp−i), the system is considered to have reached that damage state (Fail_sys), such that:

$$P\left[ {Fai{l_{{\text{sys}}}}} \right]=P\left[ {\bigcup\limits_{{i=1}}^{n} {Fai{l_{{\text{comp}} - i}}} } \right]$$

(6)

Detailed information on developing system fragility curves can be found in previous studies (Nielson 2005; Nielson and DesRoches 2007). Finally, substituting Eq. (4) into Eq. (1) gives the value of R for a given IM. Repeating this process for all relevant IM levels generates the R versus IM curve for the structure. It should be noted that the functionality recovery model used in this procedure does not account for the specific influence of repair or retrofit methods on the recovery process. If such information is available, recent approaches (e.g. Sharma et al. 2018) can be used to integrate these factors into the seismic resilience quantification.

2.2 Long-term seismic resilience model considering corrosion

Based on the resilience assessment method described in Sect. 2.1, and considering the randomness of earthquake occurrence and structural degradation due to corrosion over the life-cycle, this section introduces a methodology for estimating the long-term seismic resilience. Figure 2 illustrates the framework for the long-term seismic resilience assessment. This study assumes the entire lifespan of the structure (0, t_int] as the timeframe. For ordinary bridges, the t_int is basically set to 100 years (MOT 2015). Seismic events occurring within this lifespan are denoted by the index k. The non-negative random variables {T₁, T₂,…, T_k} represent the times when seismic events occur. The time intervals between any two consecutive seismic events are denoted by {W₁, W₂,…, W_k}. The relationship among the aforementioned time variables gives $\:{{T}}_{{{k}}} = {{W}}_{1} + {{W}}_{2} + \cdots + {{W}}_{{{k}}}$. This study utilizes the stochastic Poisson process to integrate the uncertainties associated with the occurrence and intensity of seismic activities for long-term forecasting (Matthews 2002; Li et al. 2020d). In the Poisson process, the total number of seismic events N(t_int) follows a Poisson distribution. The time-interval time W follows an exponential distribution, and the probability density function of the inter-arrival time W is given by:

$${f_W}\left( x \right)=\lambda \exp \left( { - \lambda x} \right)$$

(7)

where $\:\lambda$ is the constant occurrence rate, and $\:\lambda$ herein denotes the annual rate of exceedance conditioned on an IM level.

Fig. 2

Schematic representation of the long-term seismic resilience assessment framework

Based on the established stochastic model, the seismic resilience of the bridge, both without and with the consideration of corrosion associated with each seismic occurrence, can be described by {R₁, R₂,…, R_k} and {R_c1, R_c2,…, R_ck}, respectively. The reduction in resilience due to corrosion can be defined as a series of random variables, denoted by {$\Delta$R₁, $\Delta$R₂,…, $\Delta$R_k}. Consequently, for the k-th seismic event, the associated resilience considering corrosion can be described as.

$${R_{ck}}={R_k} - \Delta {R_k}$$

(8)

Assuming that the necessary interventions are appropriately implemented after each seismic event to restore the functionality of the structure to its undamaged state (Decò et al. 2013; Pang and Wang 2021a), the resilience without accounting for corrosion can be determined using the method described in Sect. 2.1. In this context, the objective is to assess the long-term resilience considering corrosion by cumulating the values of {R_ck}. The challenge lies in quantifying the decrease in seismic resilience {$\Delta$R₁, $\Delta$R₂,…, $\Delta$R_k} due to corrosion. According to Fig. 2, the occurrence of seismic events is a stochastic process and the associated $\Delta$R_k is random, but the values of $\Delta$R increase over time as corrosion progresses (e.g., $\Delta$R_k > $\Delta$R_k−1). To incorporate such time-dependent characteristics, the random $\Delta$R_k can be modeled by the expected reduced resilience $\:\text{E}{[ \Delta {R_k}]}$, and it can be obtained by multiplying the time interval W_k and the average annual resilience reduction $\:\text{E}{[ \Delta {R_a}]}$ due to corrosion. Thus Eq. (8) can be rearranged as.

$${R_{ck}}={R_k} - {W_k}E[\Delta {R_a}]$$

(9)

Since the life-cycle for bridge t_int is usually taken as 100 years, if $\:\text{E}{[ \Delta {R_a}]}$ is determined by taking the average of the yearly reductions in resilience, the computation can be extremely expensive and time-consuming. To simplify the computational process, instead of calculating the annual resilience decrease due to corrosion, a time interval $\Delta$t_c (taken as 10 years) is selected to assess the reduced resilience due to corrosion. It should be noted that this simplification is for numerical simplification. If sufficient computational capability is available, $\Delta$t_c can be set to 1 year to compute the annual reduction. Based on the defined interval, the number of intervals N_c selected to assess resilience reduction due to corrosion can be defined as

$${N_c}=\frac{{{t_{\operatorname{int} }}}}{{\Delta {t_c}}}$$

(10)

  

For instance, with $\Delta$t_c =10 years, the number of intervals N_c within the 100-year lifespan is 10. For the newly defined time interval, the average annual resilience reduction $\:\text{E}[ \Delta$R_a] can be calculated through:

$$E\left[ {\Delta {R_a}} \right]=\frac{{\Delta {R_{{t_{\operatorname{int} }}}}}}{{{N_c}\Delta {t_c}}}$$

(11)

where $\Delta {R_{{t_{\operatorname{int} }}}}$is the resilience reduction due to corrosion in the final year of the service life of the bridge. In this study, $\Delta {R_{{t_{\operatorname{int} }}}}=\Delta {R_{100}}$. Additionally, the average annual resilience reduction is denoted as $\Delta {R_{a\_avg}}$, and it is assumed to follow a probability distribution with a mean of $\:\text{E}{[ \Delta {{R}}_{{a}}]}$. Meanwhile, the uncertainty of R_k is also considered, and an appropriate probability distribution is assigned to it. In this context, the seismic resilience considering corrosion $\:{{R}_{ck}}$ can be solved according to Eq. (9) by using MCS. With the m-th MCS, the reduced seismic resilience $\:{{R}}_{{ck}}^{\text{(}{m}\text{)}}$ can be computed as

$$R_{{ck}}^{{(m)}}=R_{k}^{{(m)}} - W_{k}^{{(m)}}\Delta R_{{a\_avg}}^{{(m)}}$$

(12)

Consequently, the long-term seismic resilience (R_LT) can be calculated through the m-th MCS as.

$$R_{{LT}}^{{(m)}} = \sum\limits_{{k = 1}}^{{N(t_{{\text{int} }} )^{{(m)}} }} {R_{{ck}}^{{(m)}} }$$

(13)

In this study, the total number of MCS is 10⁵. The final long-term seismic resilience is calculated as the mean of the results from these 10⁵ MCS. Figure 3 illustrates the flowchart of the long-term seismic resilience assessment framework.

Fig. 3

Flowchart of the framework for assessing the long-term seismic resilience

3 Corrosion-induced deterioration modelling

Chloride-induced corrosion is a major cause of degradation in RC bridges (Cui et al. 2018; Zhou et al. 2024). When chloride ions penetrate RC members, they reduce the cross-sectional area and mechanical properties of the reinforcing steel, as well as the strength of the concrete. This section describes the quantitative methods adopted for assessing the deterioration of reinforcing steel and concrete due to chloride-induced corrosion.

3.1 Corrosion initiation

The infiltration of chloride ions into RC structures primarily occurs through the diffusion process (Li et al. 2024). After being exposed for a certain duration, chloride ions penetrate the concrete cover, compromising the passivation layer of the embedded steel, and triggering the onset of corrosion (Pinto et al. 2024). This process subsequently results in the deterioration of reinforced concrete structures. One of the most critical parameters influencing the corrosion process in reinforced concrete (RC) is determining the initial time of corrosion (T_init). The commonly used probabilistic model for estimating T_init can be expressed as (Engelund et al. 2000):

$${T_{init}}={\left[ {\frac{{{x^2}}}{{4{k_e}{k_c}{D_0}\left( {t_{0}^{n}} \right)}}{{\left[ {er{f^{ - 1}}\left( {1 - \frac{{{C_{cr}}}}{{{C_s}}}} \right)} \right]}^{ - 2}}} \right]^{\frac{1}{{1 - n}}}}$$

(14)

where x is the thickness of the concrete cover; k_e and k_c are the environmental factor and curing factor, respectively; D₀ denotes the reference diffusion coefficient, determined via compliance testing; t₀ denotes the age of concrete, measured in years, at the time the compliance assessment is carried out; n represents the age exponent that accounts for densification of the cement paste resulting from ongoing hydration; erf (·) is the Gaussian error function; C_cr refers to the chloride concentration threshold when corrosion initiates, and C_s signifies the equilibrium concentration of chloride present on the exposed surface of the concrete, and can be calculated through:

$${C_s}={A_{cs}}\left( {w/b} \right)+{\varepsilon _{cs}}$$

(15)

where A_cs is the regression parameter for chloride content at the surface, ε_cs is the error term, and w/b denotes the water-binder ratio. The probability distributions of parameters estimating the initial time of corrosion (T_init) are listed in Table 2.

Table 2 Probabilistic distribution of parameters for calculating the corrosion initiation time

3.2 Deterioration model for reinforcing steel

Once corrosion starts in the reinforcement, the corrosion rate becomes the primary factor determining the speed of its propagation. The corrosion rate (r_corr(t_p), in mm/year, after t_p years of corrosion can be empirically expressed as (Vu and Stewart 2000):

$${r_{corr}}\left( {{t_p}} \right)=0.0116 \times {i_{corr}}\left( {{t_p}} \right)$$

(16)

where i_corr (t_p) refers to the corrosion rate in µA/cm² and can be calculated through:

$${i_{corr}}\left( {{t_p}} \right)=0.85{i_{corr,0}}{t_p}^{{ - 0.29}}$$

(17)

where i_corr,0 represents the corrosion rate once the corrosion starts propagating, and is given by:

$${i_{corr,0}}=\frac{{37.8{{\left( {1 - w/c} \right)}^{ - 1.64}}}}{x}$$

(18)

The deterioration process in reinforcing steel can be categorized into two forms: uniform corrosion and pitting corrosion (Du et al. 2005a; Kashani et al. 2014; Ghosh and Sood 2016; Cui et al. 2018; Shekhar et al. 2018; Li et al. 2020a; Pinto et al. 2024). A schematic diagram illustrating different types of corrosion mechanisms of reinforcement is shown in Fig. 4. More detailed descriptions of the two corrosion models are provided below.

Fig. 4

Schematic diagram of combined corrosion mechanisms for reinforcing steel: a uncorroded rebar, b uniform corrosion, c pitting corrosion, and d combined uniform and pitting corrosion

In the case of uniform corrosion, it is assumed that the reduction in cross-sectional area of the reinforcing steel is consistent along its circumference and length following the onset of corrosion (see Fig. 4b). The residual cross-sectional area of the steel rebar under uniform corrosion, A_u(t), can be calculated using the following relationship:

$${A_u}\left( t \right)=\frac{{\pi {{\left[ {{d_0} - 2\int_{{{T_i}}}^{t} {{r_{corr}}\left( {{t_p}} \right)d{t_p}} } \right]}^2}}}{4}$$

(19)

where t = T_i + t_p is the service time, and d₀ is the diameter of the pristine rebar.

In addition to the uniform area loss observed in steel rebars, the issue is further worsened by localized corrosion at various points along the rebars. This leads to the creation of deep pits, a phenomenon also referred to as pitting sectional area loss. This study employs the commonly-used hemispherical pit formation method to determine the residual steel area (denoted as A_p(t), see Fig. 4d) (Ghosh and Sood 2016; Broomfield 2023). The mathematical expression for A_p(t) is given in Eq. (20).

$${A_p}\left( t \right)=\left[ {1 - \frac{{a\left( t \right)}}{{2{d_0}}}} \right]\left[ {{A_u}\left( t \right) - {A_0}} \right]+{A_{DP}}\left( t \right)$$

(20)

where A₀ is the area of the pristine rebar, a(t) and A_DP(t) represent the width of the pit and residual area, respectively, without considering uniform section loss, as depicted in Fig. 4c, and can be calculated through:

$$a\left( t \right)=2p\left( t \right)\sqrt {1 - {{\left( {\frac{{p\left( t \right)}}{{{d_0}}}} \right)}^2}}$$

(21)

$${A_{DP}}\left( t \right)=\left\{ {\begin{array}{*{20}{c}} {{A_0} - {A_1} - {A_2}{\text{ for }}p\left( t \right) \leqslant \frac{{{d_0}}}{{\sqrt 2 }}{\text{ }}} \\ {{A_1} - {A_2}{\text{ for }}{d_0}{\text{ }}} \end{array}} \right.{\text{ }}$$

(22)

where the pit depth p(t), and areas A₁ and A₂ can be represented as:

$$p(t) = R\int_{{T_{i} }}^{t} {r_{{corr}} (t_{p} )dt_{p} }$$

(23)

$$A_{1} = 0.5\left[ {\theta _{1} (0.5d_{0} )^{2} - a(t)\left| {0.5d_{0} - \frac{{p(t)^{2} }}{{d_{0} }}} \right.} \right];\theta _{1} = 2\arcsin \left[ {\frac{{a(t)}}{{d_{0} }}} \right]$$

(24)

$$A_{2} = 0.5\left[ {\theta _{2} (p(t))^{2} - a(t)\frac{{p(t)^{2} }}{{d_{0} }}} \right];\theta _{2} = 2\arcsin \left[ {\frac{{a(t)}}{{2p(t)}}} \right]$$

(25)

where R is the pitting factor, representing the ratio of maximum pit depth to average depth under uniform corrosion. In this study, R is assumed to follow the Extreme Value Type I Gumbel distribution with parameters µ₀ = 5.56 and α₀ = 1.16, respectively (Stewart and Al-Harthy 2008; Pinto et al. 2024).

Besides diminishing the cross-sectional area of the steel, pitting has been shown in previous studies to also reduce the yield strength and ultimate strength of the steel (Du et al. 2005a, b). The reduced yield strength and ultimate strength, represented as f_y(t) and f_u(t) respectively, can be determined using Eqs. (26) and (27).

$${f_y}\left( t \right)=\left[ {1 - 0.005{Q_{corr}}\left( t \right)} \right]{f_y}\left( 0 \right)$$

(26)

$${f_u}\left( t \right)=\left[ {1 - 0.005{Q_{corr}}\left( t \right)} \right]{f_u}\left( 0 \right)$$

(27)

where f_y (0) is the yield strength of uncorroded rebars, f_u (0) is the ultimate strength of the uncorroded rebars, Q_corr denotes the percentage area loss of reinforcing steel induced by corrosion at time t (Panchireddi and Ghosh 2019), and can be calculated as:

$${Q_{corr}}\left( t \right)=\frac{{{A_0} - {A_r}\left( t \right)}}{{{A_0}}} \cdot 100$$

(28)

where A_r (t) represents the time-varying residual cross-sectional area of reinforcement resulting from corrosion. Similarly, the corrosion also leads to a reduction in the elastic modulus of the steel reinforcement (Pinto et al. 2024). The reduced elastic modulus at time t is given by Eq. (29):

$${E_s}\left( t \right)=\left[ {1 - 0.01 \cdot {Q_{corr}}\left( t \right)} \right]{E_{s0}}$$

(29)

where E_s0 refers to the elastic modulus of the pristine rebars. Additionally, due to the influence of corrosion, the volumetric stirrup ratio (ρ_s (t)) of RC members will also change over time. The consideration of corrosion in stirrups is related to the deterioration model of longitudinal reinforcement, and the relationship between ρ_s (t) and Q_corr is as follows (Firouzi et al. 2020):

$${\rho _s}\left( t \right)=\left[ {1 - 0.3113 \cdot {Q_{corr}}\left( t \right)} \right]{\rho _{s0}}$$

(30)

where ρ_s0 is the pristine volumetric stirrup ratio, and is given by (Mander et al. 1988) in Eq. (31)

$${\rho _{s0}}=\frac{{4{A_{sp}}}}{{{d_s}s}}$$

(31)

where A_sp denotes the cross-sectional area of the stirrups, and d_s represents the diameter of spiral between bar centers, and s is the center-to-center spacing of the stirrups.

3.3 Deterioration model for concrete

The expansive force caused by the rust products of the rebars leads to the formation of micro-cracks in the cover concrete. This subsequently results in secondary effects such as a reduction in strength for both the cover and core concrete, as well as the decreased compressive strain and ultimate strain of the core concrete (Otárola et al. 2022; Pinto et al. 2024; Zhou et al. 2024). The loss of cover concrete strength can be considered using models proposed in previous studies (Coronelli and Gambarova 2004; Rao et al. 2017), as follows:

$$f_{c}^{\prime }\left( t \right)=\frac{{f_{c}^{\prime }\left( 0 \right)}}{{1+\kappa \frac{{{\varepsilon ^*}\left( t \right)}}{{{\varepsilon _{cv}}}}}}$$

(32)

where $\:{{f}}_{\text{c}}^{{\prime\:}}\left(\text{0}\right)$ denotes the compressive strength of the cover concrete in uncorroded RC bridge piers, κ is a coefficient influenced by the diameter and surface texture of the reinforcing bar and is assumed to be 0.1 for ribbed bars of medium diameter, ε_cv is the concrete compressive strain at the stress $\:{{f}}_{\text{c}}^{{\prime\:}}\left(\text{0}\right)$, and $\:{ \varepsilon }^{\text{*}}\left(\text{t}\right)$ represents the mean tensile strain across the transverse direction that initiates microcracking as a consequence of corrosion, as defined in Eq. (33):

$${\varepsilon ^*}\left( t \right)=\frac{{{n_{bars}}{w_{cr}}\left( t \right)}}{{{D_c}}}$$

(33)

where n_bar denotes the number of reinforcing bars in the cross-section of the bridge pier, D_c represents the diameter of the cross-section, and w_cr(t) signifies the cumulative crack width associated with a specific level of corrosion, which can be determined through Eq. (34) as suggested by Molina et al. (1993):

$${w_{cr}}=4\pi \left( {{\nu _{rs}} - 1} \right){r_{rs}}$$

(34)

where v_rs represents the volumetric expansion ratio of rust products relative to the volume of the uncorroded steel rebar, which is assumed to be 2 according to Molina et al. (1993), r_rs indicates the extent of corrosion penetration, which is equal to the reduction in the radius of the bar.

Due to the corrosion of stirrups, the compressive strength of the core concrete will decrease. The reduced compressive strength ($\:{{f}}_{\text{c}\text{c}}^{{\prime\:}}\left({t}\right)$) of the degraded confined concrete is considered in this study based on the theoretical stress-strain model proposed by Mander et al. (1988), as follows:

$$f_{{cc}}^{\prime }\left( t \right)=f_{c}^{\prime }\left( 0 \right)\left( {2.254\sqrt {1+\frac{{7.94{f_l}\left( t \right)}}{{f_{c}^{\prime }\left( 0 \right)}}} - 2\frac{{{f_l}\left( t \right)}}{{f_{c}^{\prime }\left( 0 \right)}} - 1.254} \right)$$

(35)

where f_l (t) represents the time-varying lateral confining pressure exerted on concrete and can be obtained from:

$${f_l}\left( t \right)=\frac{1}{2}{k_e}{\rho _s}\left( t \right){f_y}\left( t \right)$$

(36)

where k_e is the confinement effectiveness coefficient and can be determined for circular hoops using the following equation:

$$k_{e} = \frac{{\left( {1 - \frac{{s^{\prime}}}{{2d_{s} }}} \right)^{2} }}{{1 - \rho _{{cc}} \left( t \right)}}$$

(37)

where s’ signifies the clear vertical spacing between the hoops, and ρ_cc(t) is the ratio of the area of longitudinal reinforcement to the area of the core.

Additionally, the time-dependent strain at maximum strength ($\:{ \varepsilon }_{{cc}}\left({t}\right)$) and the crushing strength ($\:{ \varepsilon }_{{cu}}\left({t}\right)$) of the core concrete can be obtained as per Mander et al. (1988) through Eqs. (38) and (39), respectively.

$${\varepsilon _{cc}}\left( t \right)={\varepsilon _{cv}}\left( {1+5\left[ {\frac{{f_{{cc}}^{\prime }\left( t \right)}}{{f_{c}^{\prime }\left( 0 \right)}} - 1} \right]} \right)$$

(38)

$${\varepsilon _{cu}}\left( t \right)=0.004+1.4\left[ {\frac{{{\rho _s}\left( t \right){f_y}\left( t \right){\varepsilon _{su}}}}{{f_{{cc}}^{\prime }\left( t \right)}}} \right]$$

(39)

where $\:{ \varepsilon }_{\text{su}}$ represents the strain corresponding to the peak strength of the stirrups. The deterioration model used in this study has been widely employed in previous research (e.g. Ghosh and Padgett 2010; Ghosh and Sood 2016; Shekhar et al. 2018; Panchireddi and Ghosh 2019; Pinto et al. 2024) for life-cycle vulnerability assessment of corroded bridges. Future research also needs to account for other specific effects, such as the influence on ductility, low-cycle fatigue, and compression buckling, into the modeling of corroded reinforcing bars (Kashani 2024; Kashani et al. 2024), in order to capture more accurately the cumulative damage in corroded bridges.

4 Structural configurations and numerical models

4.1 Details of benchmark bridges

Three four-equal-span RC continuous highway girder bridges, with different span lengths and pier heights, are selected as benchmark bridges in this study. Figure 5 depicts the geometric configuration of the three bridges, with span lengths of 20 m, 30 m, and 40 m. The deck consists of a single-box single-cell section, with a width of 8.5 m and a depth of 1.9 m. The deck is supported at both external far ends by two piers with a diameter of 1.2 m (i.e. P1, P2 and P6, P7 in Fig. 5), while the interior spans rest on single-column piers with a diameter of 1.5 m (P3, P4, P5). The deck is connected to each pier using spherical steel bearings.

Fig. 5

Geometric configuration of the selected benchmark bridges

Unidirectional sliding bearings (B1 and B6 in Fig. 5) are installed on Piers P1 and P6, bi-directional sliding bearings (B2, B3, B5, and B7) are placed on Piers P2, P3, P5, and P7, and fixed bearings are used to connect the deck to the middle pier. The piers are built on a 5.6 m × 4.4 m pile caps. The cap for double-column piers is supported by 3 × 3 pile groups, while the cap for single-column piers rests on 4 × 4 pile groups. Each pile has a diameter of 0.6 m. The deck is constructed from concrete with an axial compressive strength of 32.4 MPa (MOT 2018). The piers are built using concrete that has an axial compressive strength of 26.8 MPa, reinforced with 25 mm-diameter longitudinal bars (yield strength of 400 MPa) and 16 mm-diameter stirrups placed at 0.1 m vertical intervals. The foundations are composed of concrete with a compressive strength of 23.4 MPa.

4.2 Finite element modeling

The finite element models of the benchmark bridges are constructed using the nonlinear analysis platform OpenSees (McKenna 2011), as illustrated by the overview in Fig. 6. The deck is modeled using elastic beam-column elements, with the sectional properties of the box-girder assigned to these elements. The total mass of the deck is distributed at each node, with the translational masses at each node calculated based on the sectional properties and element length. The spherical steel bearings are modeled using zero-length elements. For the sliding bearings, the horizontal force-displacement relationship is simulated using an elastic-perfectly plastic material object. The initial stiffness (k_e) of the bearings is calculated using the relationship k_e = µ·F_v/d_y, where µ is the coefficient of friction, taken as 0.02 (MOT 2020), F_v is the vertical reaction force, and d_y is the sliding-initiated displacement, taken as 0.003 m. The vertical constitutive relationship of the bearings is modeled using an elastic uniaxial material object. Given that the vertical displacement of the bearings is minimal, a large value is assigned to the vertical stiffness (k_v). For the fixed bearings, an elastic uniaxial material object is used for modeling. A stiffness of 1 × 10⁸ kN/m (Feng et al. 2024) is assigned to each translational direction (i.e., k_h and k_v), ensuring minimal displacement under seismic loads. The modeling approaches used for the bridge piers are similar to those adopted in previous studies (Pinto et al. 2024; Zhou et al. 2024), which have been validated against experimental tests or in-situ data. The piers are simulated using nonlinear beam-column elements with fiber sections. Within the fiber sections, the stress-strain relationship of the reinforcement is modeled using the Steel02 material (Filippou et al. 1983), while the constitutive relationships for the cover and core concrete is defined using the Concrete04 material (Mander et al. 1988). The pile cap is simulated using elastic beam-column elements. The translational nodal masses are determined based on the sectional properties and element length and lumped at their respective nodes. Soil-structure interaction is not considered herein, hence the nodes at the bottom of the pile cap have fixed constraints. The influence of deterioration caused by corrosion is considered by updating the model with reduced parameters (e.g., diameter and strength of the rebar, and compressive strength of concrete) as determined from the corrosion model.

The implicit Newmark integration scheme (Newmark 1959) is employed for the dynamic analysis of the bridges, with coefficients β and γ set to 0.25 and 0.5, respectively. The bridge models incorporate Rayleigh damping, with the mass and stiffness proportional coefficients (α_M and β_K) determined using the circular frequencies of the first mode and the sum of modal effective masses reaching 90%, along with a modal damping ratio of 0.05. Uncertainties in the modeling parameters are considered in the simulation process. Based on the probabilistic distributions in Table 3, the Latin hypercube technique (Olsson et al. 2003) is used for random sampling to obtain modeling parameter samples.

Fig. 6

Finite element model of the benchmark bridges

Table 3 Consideration of uncertainties in bridge modeling parameters

5 Ground motion sequences

For the seismic resilience assessment, a sufficient number of representative ground motion sequences should be selected based on the earthquake scenario and the seismic hazard specific to the bridge site (Feng et al.2021). In this study, the benchmark bridges are assumed to be located in California in a site of Soil Class C according to the NEHRP site classification (NEHRP 2001), characterized by an average shear wave velocity (V_S30) in the top 30 m of 480 m/s. The Joyner-Boore distance (R_jb) is 5 km, and the fault type of the rupture is strike-slip. Based on the aforementioned target earthquake scenario, the earthquake ground motion records matching an unconditional spectrum are selected from the PEER NGA-West2 Database (PEER 2013) using the method proposed by Baker and Lee (Baker and Lee 2018). Although many ground motion records satisfy the target earthquake rupture scenario, those that include aftershocks are very limited. Ground motion from events with aftershocks are selected and the mainshock and aftershock are matched accordingly. Considering the inherent correlation and differences in ground motion characteristics between mainshocks and aftershocks (Chiou and Youngs 2008; Ruiz-García and Negrete-Manriquez 2011; Abrahamson et al. 2013), this study defines the following mainshock-aftershock matching criteria: (i) the mainshock and the aftershock are from the same seismic event, (ii) the mainshock and the aftershock are recorded at the same station, and (iii) the magnitude of the mainshock is greater than or equal to that of the aftershock.

Based on the above criteria, 40 natural mainshock and aftershock ground motion records are selected. Table 4 summarizes the main information regarding the selected MS-AS ground motion sequences. From Table 4, it is evident that the overall selected ground motion records have relatively low intensity. To produce larger inelastic responses, in addition to the natural unscaled records, these ground motion acceleration time histories are also uniformly amplified by a factor of 2, resulting in an additional 40 ground motion records. It is worth noting that this scale factor is generally acceptable from an engineering practice perspective (Bommer and Acevedo 2004; Watson-Lamprey and Abrahamson 2006). Accordingly, a total number of 80 MS-AS record pairs are used for the subsequent seismic fragility and resilience assessment. The range of magnitudes (M_w) for the mainshocks varies from 5.77 to 7.62, with seismic intensity, represented in terms of PGA in this study, ranges from 0.1 g to 1.79 g. In contrast, the magnitudes of the aftershocks range from 4.80 to 6.30, with PGAs between 0.04 g and 1.25 g. Figure 7 illustrates the response spectra of the 80 mainshock and aftershock ground motions. Within the analysis, a time gap of 30 s is set between the application of the mainshock and aftershock ground motion (Wen et al. 2018) in order to ensure that the bridge stops oscillating before the aftershock phase.

Table 4 List and information on selected ground motion sequences

Fig. 7

Acceleration response spectra of the selected ground motion sequences: a mainshock ground motions and b aftershock ground motions

6 Seismic performance and fragility assessment of deteriorating bridges under MS-AS sequences

6.1 Assessment procedures

In the seismic performance assessment of bridges, the piers are primary components (Ramanathan et al. 2015) which typically govern the vertical stability and load-carrying capacity of the structure. Due to the continuous loading process of possible mainshock and aftershock sequences, the structure may undergo multiple loading cycles, resulting in cumulative damage. This study primarily focuses on the cumulative damage of bridge piers. A well-established damage quantification metric, the Park-Ang damage index (Park and Ang 1985) (denoted as DI_PA hereafter), is utilized and can be expressed as follows:

$$D{I_{PA}}=\frac{{{\delta _m}}}{{{\delta _u}}}+\frac{\beta }{{{Q_y}{\delta _u}}}\int {dE}$$

(40)

where δ_m represents the maximum deformation observed during nonlinear dynamic analysis, while δ_u denotes the ultimate deformation under monotonic loading; Q_y stands for the yield strength, dE signifies the incremental absorbed hysteretic energy, and β is a non-negative parameter that indicates the structural damage induced by cyclic loading. The expression for β is given by (Park and Ang 1985):

$$\beta =\left( { - 0.447+0.073\frac{l}{d}+0.24{n_0}+0.314{\rho _t}} \right) \times {0.7^{{\rho _w}}}$$

(41)

where l/d represents the shear span ratio, n₀ is the normalized axial load, ρ_t is the longitudinal reinforcement ratio as a percentage, and ρ_w stands for the confinement ratio. Based on previous research (Panchireddi and Ghosh 2019; Zhou et al. 2023), this study assumes that the Park-Ang damage index is applicable to corroded structures. It is worth noting that other studies (Afsar Dizaj and Kashani 2020) suggest that future research should consider further assessment of the cumulative damage index to account for the effects of corrosion.

Due to the consideration of deterioration induced by corrosion over the entire service life (i.e., 100 years), this study divides the lifespan of the bridge into intervals of 10 years, starting from Year 0 (uncorroded) to Year 100. Additionally, the study examines two seismic excitation scenarios: mainshock-aftershock (MS-AS) and mainshock only (MS), in order to investigate the effects of aftershocks on the cumulative damage. Cloud Analysis (Shome et al. 1998) is used in the nonlinear time history analysis (NTHA) to characterize the relationship between the cumulative damage and seismic intensity. In total, 22 analysis scenarios (2 excitation scenarios × 11 corrosion scenarios) are considered, resulting in 5, 280 NTHAs (80 ground motions × 3 benchmark bridges × 22 analysis scenarios). Based on the results of nonlinear dynamic and static analyses, the cumulative damage and fragility of benchmark bridges are assessed, as discussed in detail in Sect. 6.2 and 6.3 below.

6.2 Cumulative seismic damage evaluation

This section separately examines the impact of aftershocks and corrosion-induced deterioration on the cumulative damage to the bridge. To study the influence of aftershocks alone, the structural response in the uncorroded scenario is primarily analyzed. Using Pier P4 which has the highest cumulative damage among all the piers as a representative, Fig. 8 illustrates the comparative results of the average DI_PA over different service times for MS-AS and MS only scenarios. It can be observed from the figure that the aftershock can exacerbate the cumulative damage to the bridge. This can be explained through Eq. (40), in which the maximum deformation ratio δ_M/δ_u is represented by r_δ, and the energy dissipation $\beta /(Q_{y} \delta _{u} )\int {dE}$ is denoted by r_E. Figure 9 compares r_δ and r_E of benchmark bridges under MS only and MS-AS conditions. It is shown that, overall, the value of r_δ is greater than that of r_E, indicating that the peak displacement of the pier during earthquakes predominantly influences the DI_PA. Since the peak is associated with the spectral response, Fig. 10 compares the average response spectra of the mainshock and aftershock for strong motion records where r_d under MS-AS is significantly greater than those under MS only (i.e., Record nos. 2, 35, 37, 38, 42, 43, 75, 76, 77, 78, 79). It can be observed that the average spectral acceleration of the aftershock corresponding to the fundamental period of the bridge is significantly greater than that of the mainshock. This shows that the influence of aftershocks on the cumulative damage to bridge piers depends on the difference in spectral seismic intensity between the aftershock and the mainshock. The higher intensity of the aftershock compared to the mainshock is the primary cause of the increased cumulative damage.

Fig. 8

Comparison of the average cumulative damage index of the middle pier of different benchmark bridges under the MS and MS-AS over their lifetime: a Bridge #1, b Bridge #2, and c Bridge #3

Fig. 9

Comparison of maximum deformation ratio (r_d) and energy dissipation (r_E) of benchmark bridges under MS only and MS-AS: ar_d and br_E

Fig. 10

Comparison of the average response spectra of the mainshock and aftershock for strong motion records where r_d under MS-AS is significantly greater than those under MS only

To further quantify the influence of aftershocks on cumulative damage, an additional indicator, denoted the cumulative damage ratio (r_c), is defined as follows:

$${r_c}=\frac{{D{I_{PA(MS - AS)}}}}{{D{I_{PA(MS)}}}}$$

(42)

in which DI_PA(MS−AS) and DI_PA(MS) are the average Park-Ang damage index under MS-AS and MS only, respectively. Figure 11 depicts the average cumulative damage ratio (r_c) of the benchmark bridges over their entire lifespan. It is observed that the aftershock has a greater influence on the cumulative damage for bridges with relatively tall piers compared to those with comparatively short piers, as a consequence of their dynamic characteristics. To further illustrate this, Fig. 12 presents the average response spectra of all the mainshock and aftershock ground motion records. The fundamental periods of the three benchmark bridges vary within specific ranges. Specifically, Bridge #1 exhibits a fundamental period between 0.53 and 0.60 s, Bridge #2 falls between 1.26 and 1.44 s, and Bridge #3 ranges between 2.14 and 2.52 s. As shown in Fig. 12, for relatively long-period bridges, the spectral accelerations under MS and AS corresponding to the fundamental period tend to become comparable. Therefore, the increase in the relative cumulative damage due to aftershocks is greater for Bridge #2 and Bridge #3 compared to Bridge #1. Although Bridge #3 has the longest period, the response spectrum indicates that the intensity of the seismic excitation is very low. Moreover, as shown in Fig. 9c, the bridge remains largely in an elastic state under most MS and AS excitation scenarios. However, for Bridge #2 (see Fig. 9b), aftershocks still cause more energy dissipation than the mainshocks under some strong motion records. Therefore, the influence of aftershocks on the relative cumulative damage to Bridge #2 is greater than that for Bridge #3. It should be noted that the possible influence of low-cycle fatigue of longitudinal rebars is not considered herein in the numerical modeling of bridge piers. This factor may exacerbate the cumulative damage of the structure under repeated cyclic loading, as suggested in recent studies (Salami et al. 2019; Afsar Dizaj et al. 2022), potentially amplifying the influence of aftershocks on cumulative damage.

Fig. 11

Average cumulative damage ratio (r_c) of the benchmark bridges over their entire lifespan

Fig. 12

Average response spectra of all the mainshocks and aftershocks

Figure 11 also illustrates the impact of corrosion-induced deterioration on the cumulative damage under MS-AS sequences. It is indicated that, overall, as the level of corrosion in bridges intensifies over the years, the influence of aftershocks on the cumulative damage ratio increases. For Bridges #1, #2, and #3, the r_c values increase from 1.53 to 1.92, 2.59 to 3.07, and 2.06 to 2.12, respectively. This is because as the corrosion level increases, the fundamental period elongates, and the spectral response values for MS and AS corresponding to the fundamental period tend to converge. From the comparative average response spectra for MS and AS (see Fig. 12), it is evident that for Bridge #3, the change in the bridge period due to corrosion has a minimal impact on spectral acceleration. Therefore, the influence of corrosion on the seismic damage in Bridge #3 is negligible.

6.3 Time-dependent fragility curves

To establish the time varying fragility curves under different seismic excitation scenarios, the capacity limit states of the bridge piers are also required. Table 5 shows the damage states (DS) considered in this study along with their corresponding ranges of DI_PA. It is worth noting that this study draws on the damage state definitions used in previous relevant studies (Ghosh et al. 2015; Panchireddi and Ghosh 2019), based on criteria proposed by Park et al. (1985). The effect of time-dependent corrosion (Dizaj et al. 2018) on DS definitions, may require further consideration in future work to assess the possible influence on the cumulative damage index.

Table 5 Definition of different damage states of bridge piers based on DI_PA (Park and Ang 1985; Ghosh et al. 2015; Panchireddi and Ghosh 2019)

According to Eqs. (5) and (6), the system fragility curves for different years under MS and MS-AS conditions are developed based on the cumulative damage of each pier and the defined DS. Using the uncorroded scenario, the 50th year, and the 100th year as representative cases, Fig. 13 shows the fragility curves of the benchmark bridges. Overall, as the bridge period increases, the probability of damage (denoted as P_f [DS|IM]) for the bridge decreases, hence the likelihood of complete damage for Bridge #2 and Bridge #3 is extremely low.

In general, as indicated in the figures, the aftershocks significantly increase the fragility of the bridge system. As the damage state worsens, the influence of the aftershock gradually decreases. Using Bridge #1 in the uncorroded scenario as an example, the maximum differences in system fragility between MS and MS-AS for various DS are 0.19, 0.16, 0.14, and 0.10, respectively. This result indicates that bridges in lower DSs are severely affected by earthquake sequences. Also, the effect of aftershocks on fragility increases with the passage of service time. For instance, for Bridge #1 in the slight DS, the maximum difference in fragility between MS and MS-AS is 0.19 for the uncorroded scenario. When the earthquake occurs in the 50th and 100th years, this difference increases to 0.25 and 0.28, respectively. Furthermore, corrosion-induced structural deterioration increases the probability of damage to the bridge system. However, as the service time increases, the influence of corrosion on the fragility diminishes. For example, for Bridge #1 in a moderate DS, the maximum difference in fragility under MS is 0.24 between the uncorroded scenario and the 50th year, while the difference between the 50th year and the 100th year is 0.05. Similarly, for MS-AS, the maximum difference in fragility is 0.23 between the uncorroded scenario and the 50th year, and 0.07 between the 50th year and the 100th year. Additionally, as the fundamental period of the bridge increases, the influence of corrosion on fragility decreases. For instance, for the ‘extensive DS’, the maximum differences in fragility between the uncorroded scenario and the 50th year under MS and MS-AS are 0.28 and 0.28, respectively, for Bridge #1; while the corresponding values are 0.08 and 0.10 for Bridge #2; decreasing to 0.01 and 0.02 for Bridge #3.

Fig. 13

Time-dependent system fragility curves for different bridges under MS and MS-AS: a Bridge #1, b Bridge #2, and c Bridge #3

7 Long-term seismic resilience evaluation

Substituting the system fragility P_f into Eq. (4) provides the functionality recovery profiles of the benchmark bridges under different seismic intensities for both the MS and MS-AS cases. Similarly, considering the uncorroded condition, the 50th and the 100th year as representative scenarios, Fig. 14 shows the time-dependent functionality recovery curves. Overall, as the fundamental period of the bridge increases, the expected loss of functionality decreases. For the same bridge, seismic intensity has a significant influence on the functionality. As the PGA increases, the residual functionality (Q_r in Fig. 1a) of the bridge decreases, and the duration required to restore it to its target functionality lengthens. In general, the aftershocks can have a significant effect on the functionality of the bridge. For Bridge #1, when the PGA is 0.9 g, the maximum difference in expected functionality between the MS and MS-AS cases reaches 0.17. Corrosion-induced deterioration affects bridges with different fundamental periods differently. For Bridge #1, the maximum differences in expected functionality between the uncorroded condition and the 100th year are 0.23 and 0.25 under MS-only and MS-AS conditions, respectively. As the service time increases, the influence of corrosion on functionality gradually decreases. For Bridge #2 and Bridge #3, the effect of corrosion on functionality also decreases over time, with the maximum differences in expected functionality between the uncorroded condition and the 100th year being 0.08 and 0.02, respectively.

Fig. 14

Time-dependent functionality recovery profiles for different bridges under MS and MS-AS: a Bridge #1, b Bridge #2, and c Bridge #3

Fig. 15

Seismic resilience of different bridges over their lifespan under MS and MS-AS: a Bridge #1, b Bridge #2, and c Bridge #3

Fig. 16

Average corrosion rate (r_corr(t_p)) over the lifetime of the benchmark bridges

According to Eq. (1), integrating the functionality recovery curves over time for different seismic intensities of MS and MS-AS conditions at various years provides the seismic resilience (R) of the benchmark bridges over their lifespan, as shown in Fig. 15. It is found that, overall, the seismic resilience increases with the elongation of the fundamental period of the bridge. For the same bridge, higher seismic intensity leads to reduced resilience. The aftershocks reduce seismic resilience, but the extent of this effect varies among different bridges. For instance, in the uncorroded scenario with PGA = 0.3 g, the maximum difference in the R value between MS and MS-AS is 0.056 for Bridge #1, while for Bridge #2 and Bridge #3, the maximum differences are reduced to 0.046 and 0.019, respectively. Furthermore, corrosion-induced deterioration reduces the resilience of the bridge, with the magnitude of this reduction decreasing as the fundamental period of the bridge increases. In particular, for Bridge #3, the effect of deterioration can be considered negligible. It is worth noting that for Bridge #1 and Bridge #2, seismic resilience decreases most significantly during the initial 10 years of their service life, primarily due to the higher corrosion rate (r_corr(t_p)) during this time period, as shown in Fig. 16.

After determining the seismic resilience values for the bridges under various excitation scenarios at different years, and considering the risk of earthquake occurrence as well as the randomness of seismic resilience and resilience loss (ΔR) due to deterioration, the method outlined in Sect. 2.2 is used to calculate the long-term resilience (R_LT) of the benchmark bridges under MS and MS-AS conditions. In this study, it is assumed that R and ΔR follow a normal distribution, with coefficients of variation of 0.1. Moreover, three different earthquake intensity scenarios with return periods of 475 years (2% in 50 years), 975 years (5% in 50 years), and 2475 years (10% in 50 years) are considered for the investigated bridges. The corresponding occurrence rates (λ) in Eq. (7) are 0.0021, 0.001, and 0.0004, respectively. Given these return periods, the PGA values for the three earthquake intensity levels are 0.46 g, 0.6 g, and 0.8 g, respectively.

Figure 17 illustrates the results of the life-cycle resilience for various bridges. Overall, as the fundamental period of the bridge increases, the R_LT of the bridge also increases. Conversely, R_LT decreases with the increase in seismic intensity. Moreover, it is found that aftershocks have a non-negligible impact on the long-term resilience of Bridge #1 and Bridge #2, while the influence on Bridge #3 is minimal. The maximum differences in R_LT between MS-AS and MS for Bridge #1 and Bridge #2 are 9.6% and 8.7%, respectively, whereas for Bridge #3, the maximum difference in R_LT is only 4.4%. This is mainly because aftershocks induce substantial additional cumulative damage to Bridge #1 and Bridge #2, significantly reducing their seismic resilience (see Fig. 15a and b). In contrast, aftershocks have a minimal influence on the cumulative damage of Bridge #3, resulting in only a slight decrease in seismic resilience and, consequently, minimal variation in its life-cycle resilience. Additionally, it is also observed that the effect of aftershocks on the long-term resilience is marginally influenced by seismic intensity. For the investigated bridges, the differences in the reduction of resilience caused by aftershocks under different PGAs are within 3%.

Fig. 17

Life-cycle resilience (R_LT) for benchmark bridges under MS and MS-AS: a Bridge #1, b Bridge #2, and c Bridge #3

8 Conclusions

This study investigated the cumulative damage of deteriorating RC highway bridges exposed to mainshock and aftershock (MS-AS) sequences. An integrated framework for assessing the life-cycle seismic resilience was proposed and applied to representative bridge structures. The investigation included the following underlying components: (i) development of life-cycle seismic resilience assessment methodology for deteriorating structures; (ii) identification of the influence of aftershocks and corrosion-induced deterioration on cumulative damage in bridges; (iii) generation of time-dependent system fragility curves under MS-AS sequences; and (iv) implementation of life-cycle seismic resilience evaluation for deteriorating RC highway bridges under MS-AS sequences. The role of repair and maintenance in improving bridge functionality throughout its entire life cycle can be considered in future developments of the framework. The main findings of the study are summarized below.

• (1) The influence of aftershocks on the cumulative damage to bridges depends on the difference in seismic intensity between the aftershock and the mainshock. In some cases, depending on the ground motion sequence and the dynamic characteristics of the structure, it was shown that the aftershocks can significantly exacerbate the cumulative damage to the benchmark bridges.

• (2) Aftershocks were shown to have a higher influence on the relative cumulative damage in bridges with relatively taller piers (i.e. Bridges #2 and #3 in this study), which typically exhibit relatively long fundamental vibration periods. For comparatively long fundamental periods, the spectral accelerations for the MS and AS scenarios tend to become comparable, resulting in more relative cumulative damage. However, the increase in relative cumulative damage with the fundamental period is dependent on whether the structure experiences notable inelastic demands. It was shown that the bridge with the longest period (i.e. Bridge #3) remained largely elastic, hence exhibiting a lower cumulative damage ratio compared to the medium-period structure (i.e. Bridge #2).

• (3) As the level of corrosion in bridges escalates, with increasing service years, the influence of AS on the cumulative damage ratio becomes more pronounced. For the considered benchmark structures, Bridges #1, #2, and #3, the cumulative damage ratio (r_c) increased from 1.53 to 1.92, 2.59 to 3.07, and 2.06 to 2.12, respectively. This is because as corrosion progresses, the fundamental period of the bridge elongates, causing the spectral response values for the MS and AS cases corresponding to the fundamental period to converge.

• (4) Aftershocks were shown to significantly increase the fragility of the bridge system, especially for lower damage states (DSs). As the DS increases, the influence of AS gradually decreases. For Bridge #1, which exhibited the highest fragility, the maximum difference in system fragility between the MS and MS-AS scenarios for various DS were 0.19, 0.16, 0.14, and 0.10, respectively. The influence of aftershocks on the fragility also increased with the service life of the bridge. For instance, in the ‘slight DS’ of Bridge #1, the maximum difference in fragility between MS and MS-AS was 0.19 for the uncorroded scenario. When seismic events occur in the 50th and 100th years, this difference increased to 0.25 and 0.28, respectively.

• (5) Aftershocks were also shown to diminish the seismic resilience, albeit with varying effect across different bridges. For instance, under the uncorroded scenario with PGA = 0.3 g, Bridge #1 exhibited the largest difference in R value between MS and MS-AS at 0.056, while Bridge #2 and Bridge #3 showed reduced differences of 0.046 and 0.019, respectively. Additionally, corrosion-induced deterioration decreased the bridge resilience, with the magnitude of reduction diminishing as the fundamental period of the bridge increased. While the influence of deterioration on Bridge #3 was minimal, due to its largely elastic response, Bridges #1 and #2 experienced the most pronounced decline in seismic resilience within the initial 10 years of service life, primarily due to higher corrosion rates during this period.

• (6) Aftershocks were found to have a non-negligible effect on the life-cycle resilience (R_LT) of Bridges #1 and #2, while their influence on Bridge #3 was minimal. The maximum differences in R_LT between MS-AS and MS for Bridge #1 and Bridge #2 were 9.6% and 8.7%, respectively, whereas the maximum difference in R_LT for Bridge #3 was only 4.4%. The influence of aftershocks on the long-term resilience was also marginally influenced by seismic intensity. For the considered bridges, the differences in resilience reduction caused by aftershocks under different PGAs were within 3%.

Overall, the findings of this investigation highlight the merits of the proposed framework in evaluating the life-cycle seismic resilience of bridges for different hazard scenarios and deterioration conditions.