Seismic damage assessment for concrete bridge columns using maximum and residual displacement data synthesized via generative AI

Abstract

Seismic damage index for concrete bridges often quantifies only peak deformation and neglects residual response, which reflects cumulative damage and governs postearthquake recovery. Meanwhile, their use in risk-based fragility assessments renders the results highly sensitive to data scarcity, particularly for rare high-damage states under extreme earthquakes. This study develops an integrated probabilistic framework that couples damage quantification with data augmentation. A new damage index is proposed to combine maximum and residual displacements into a single interpretable measure. Optimized limit states are calibrated accordingly. To address data scarcity, a fragility-informed conditional generative adversarial network (cGAN) is introduced. The model incorporates adversarial, physics-based, and fragility-based losses to generate realistic combinations of intensity measures and engineering demand parameters for calibration and uncertainty quantification. By enriching the dataset with cGAN-generated samples, the method produces fragility curves with reduced dispersion and improved reliability in severe-damage regimes for the bridge demonstration. The integration of enhanced damage index and statistically consistent data augmentation offers a robust basis for seismic performance evaluation and supports informed decisions on inspection, repair, retrofit, and resilience planning.

Keywords

Seismic damage assessment concrete bridge columns generative AI conditional generative adversarial networks (cGANs)data augmentation damage index limit state fragility curves

Introduction

Bridges are fundamental elements of transportation infrastructure, and maintaining their safety and postevent functionality is a central objective of seismic design and asset management (ASCE, 2022; Caltrans, 2025; Housner and Thiel Jr, 1995). In highway systems, reinforced-concrete (RC) columns are the principal load-bearing elements that concentrate inelastic demand during strong ground motions and thus largely determine both collapse risk and postearthquake usability (Caltrans, 2025; Kim et al., 2005; Todorov and Billah, 2021). Effective seismic assessment and decision-making therefore require damage metrics that meaningfully connect structural responses to practical outcomes such as inspection, repair and continued service. Performance-based seismic design (PBSD) provides the procedural foundation for linking seismic intensity measures (IMs) to engineering demand parameters (EDPs) and ultimately to performance decisions (Ghobarah, 2001; Zameeruddin and Sangle, 2016). Specifically, probabilistic seismic demand models (PSDMs) characterize the statistical relationship between IMs (e.g., peak ground acceleration) and EDPs (e.g., maximum displacement and residual offsets) and can be used to derive conditional exceedance probabilities. In many studies these probabilities are expressed as fragility curves to summarize the likelihood of exceeding predefined limit thresholds (Cornell et al., 2002; Du et al., 2019, 2023; Yazdani and Yazdannejad, 2019; Yin and Sun, 2005). Fragility curves serve as a fundamental component of risk-informed seismic assessment, yet their practical reliability depends on accurate and unbiased damage quantification, particularly the definition of the scalar damage index (DI) and limit states (LS) that link PSDM outputs to engineering decisions.

Most conventional bridge RC column assessments adopt a single metric, commonly used maximum displacement (MD) or a cumulative ductility metric such as residual displacement (RD) index, to define LS (Choi et al., 2004; Kim et al., 2005; Liu and Sun, 2005; Todorov and Billah, 2021; Yoon et al., 2019; Zhong et al., 2023). A MD-based DI was proposed to quantify column damage as an explicit function of peak displacement demand and the ultimate and yield displacements from the monotonic pushover backbone (Yoon et al., 2019, 2022). The DI is linked to six principal LS, each defined based on visually observable plastic-region damage from laboratory tests, which establishes a probabilistic connection between seismic peak demands and limit states of RC columns, serving as a fundamental component of PBSD for highway bridges. Peak metrics capture transient inelastic ductility demand and are straightforward to obtain from time-history analysis, but they do not directly quantify postearthquake functionality. Residual deformations govern postearthquake functionality and repairability, even when collapse does not occur (Gan et al., 2019; Mackie and Stojadinovic, 2004; Seyed Ardakani and Saiidi, 2018; Wei et al., 2023; Xiang et al., 2020). Accordingly, engineering practice and resilience assessment increasingly require damage assessment including both severity (peak response) and usability (residual response). Existing treatments for the damage index incorporating the joint behavior of MD and RD typically either define independent thresholds for each variable and combine them into discrete damage state matrix, or assume a bivariate lognormal model with a fixed correlation coefficient (Wei et al., 2023). Such approaches rely on strong parametric assumptions and coarse discretization, limiting their ability to capture conditional dependence, skewness, and tail dependence that are critical for representing severe and functional damage behaviors (Li and Sun, 2020; Saini and Saiidi, 2014; Wu et al., 2023).

Moreover, the scarcity of experimental and high-fidelity simulation data, particularly in the severe-damage regime, constrains reliable estimation of exceedance behavior and increases uncertainty in low-probability damage states. Recent studies have applied machine learning and surrogate modeling to improve the efficiency of seismic fragility assessment for bridge systems and reinforced concrete bridge columns, including neural-network-based fragility models (Nguyen et al., 2026), Gaussian-process surrogates with active learning (Ning et al., 2024), and machine-learning-aided frameworks for response estimation and fragility construction (Liu et al., 2026). These approaches have demonstrated clear advantages in reducing the computational burden associated with nonlinear analyses. However, these approaches are primarily oriented toward learning response or fragility mappings from available data, without explicitly constraining the resulting models to preserve consistency in joint probabilistic relationship between IMs and EDP or exceedance behavior across all the damage states and are often not validated against independent datasets. This limitation becomes particularly relevant when IM-EDP samples are unevenly distributed and damage assessment depends jointly on maximum and residual responses, for which variability is inherently higher and the severe-damage regime is often underrepresented.

To address data scarcity in characterizing the joint probabilistic relationship between IMs and EDP, a conditional generative adversarial network (cGAN) is employed within a probabilistically constrained framework (Kang et al., 2025; Marano et al., 2024; Wu et al., 2023). The model is trained using a loss formulation that integrates adversarial, physics-based, and fragility-consistency terms, enabling the generated samples to preserve conditional dependence and exceedance behavior across damage states. The discriminator module in the cGAN model is further guided by an independent set of IM-EDP results obtained from simulations under recorded ground motions, which serves as a probabilistic constraint for distributional consistency during training. Within those settings, the cGAN is used to reconstruct the joint IM-EDP distribution, particularly in underrepresented severe-damage regions.

The proposed method is demonstrated through a case study of an ordinary highway bridge, where the fragility-informed cGAN, trained from the OpenSees nonlinear time-history results, generates synthetic IM-EDP datasets to calibrate and validate the proposed DI and LS framework. Combined with the statistically optimized θ_J, the proposed D_J establishes a balanced, probabilistically consistent tool for evaluating both peak and residual demands. The cGAN-augmented data enhances fragility reliability, particularly for high-damage states, and enables robust characterization of joint response statistics under data-scarce conditions. Consequently, the proposed DI and LS framework provides a parsimonious, physically interpretable link between PSDM-derived responses and decision-making for inspection, repair prioritization, retrofit planning, and resilience assessment.

Methodology

The proposed method integrates two complementary threads required for seismic damage fragility assessment: (i) a univariate fragility-equivalent damage index that summarizes both peak inelastic deformation and permanent deformation in a single, interpretable scalar metric, and (ii) a conditional generative model that augments scarce IM-EDP (intensity-response) data to stabilize estimation and calibration of exceedance of probability functions. The following subsections describe the damage-index formulation and derivation (see “Damage index and optimal limit states”) and the components of the cGAN model and the training procedures (see “cGAN augmentation”). The practical engineering implementation and tuning used in the case study are reported in “Case study of ordinary concrete highway bridge”.

Damage index and optimal limit states

MD-based and RD-based metrics are adopted as the primary DI in this study. The MD-based index can be calculated by $\frac{Δ_{D} - Δ_{Y}}{Δ_{U} - Δ_{Y}}$ ( $Δ_{D}$ : maximum displacement demand, $Δ_{Y}$ and $Δ_{U}$ : yield and ultimate displacement from pushover analysis) in Table 1. Following the probabilistic damage control application (PDCA) proposed by Saini and Saiidi (2014), the MD-based index quantifies the peak transient inelastic ductility and is directly associated with four failure modes in plastic hinge, including Major Spalling (Limit state 1), Exposed Reinforcement (Limit state 2), Core Shedding (Limit state 3), and Reinforcement Rupture (Limit state 4) in Figure 1 (Saini and Saiidi, 2014). Its primary advantage lies in its derivation from a large dataset of columns with diverse configurations and failure modes, enabling statistically robust and engineering-oriented representation of seismic damage across a wide range of scenarios (Yoon et al., 2019).

Figure 1.

Damage limit states (from Saini and Saiidi, 2014).

Table 1.

Damage limit states via performance for MD-base and RD-based DI.

DI	θ₁	θ₂	θ₃	θ₄
MD-based	0.375	0.6	0.822	1.0
RD-based	0.2%	0.5%	1.0%	1.75%

The residual drift metric, expressed as a percentage (%) of drift and adopted based on Wei et al. (2023) as summarized in Table 1, measures the permanent lateral offset remaining after ground shaking and governs postevent functionality and repairability. Four discrete RD-based limit states are defined, including Quick Repair (Limit State 1), Repairable (Limit State 2), Difficult to Repair (Limit State 3), and Reconstruction (Limit State 4), corresponding to increasing levels of residual damage.

All the LS with the DI are specified by set of thresholds θ_1,k for RD-based metric (D₁) and θ_2,k for MD-based metric (D₂) (index k enumerates damage levels). Because either peak demand or permanent deformation can independently trigger practical consequences, the joint exceedance rule used to define the target fragility surface must reflect a balanced engineering judgement between safety and economic considerations. A composite joint criterion is therefore adopted. For a chosen IM level, the joint exceedance event for limit state k is defined by the weighted normalized sum:

a_{k} \frac{D_{1}}{θ_{1, k}} + (1 - a_{k}) \frac{D_{2}}{θ_{2, k}} \geq 1

(1)

where

a_{k} \in [0, 1]

represents the relative importance of residual demand at state k. The joint criterion occupies the engineering middle ground between the extreme logical rules “AND” and “OR”: it avoids the under-conservatism of “AND” and the excessive conservatism of “OR” by requiring a balanced contribution from D₁ and D₂. The conditional joint exceedance probability conditioned on IM is defined as

P_{joint, k} (D_{join, i} \geq θ_{k} | I_{i}) = P (a_{k} \frac{D_{1}}{θ_{1, k}} + (1 - a_{k}) \frac{D_{2}}{θ_{2, k}} \geq 1 | I_{i})

(2)

where i is the number of each IM-EDP pair obtained from experiment, history record, or simulation result. In practice

P_{joint, k} (D_{join, i} \geq θ_{k} | I_{i})

is estimated either by direct frequency counting within IM bins or by fitting a smoothed parametric fragility function to the binary exceedance indicator defined by equation (1). A lognormal functional form is convenient and interpretable for positive, multiplicative observables. Thus, a lognormal fragility model is fitted in log (I) space:

P_{joint, k} (D_{join, i} \geq θ_{k} | I_{i}) = Φ (\frac{\ln (I_{i}) - μ_{k}}{β_{k}})

(3)

where Φ(⋅) is the standard normal cumulative density function (CDF), and

μ_{k}

β_{k}

are the mean and standard deviation in log-space, which is estimated by maximum likelihood estimation from the binary exceedance IM-EDP data. Accordingly, a univariate fragility-equivalent damage index D_J is constructed so that its exceedance probabilities reproduce the composite probabilities of equation (2). The interpretable candidate is a normalized linear combination:

D_{J, k} = a_{k} {\tilde{D}}_{1} + (1 - a_{k}) {\tilde{D}}_{2}, {\tilde{D}}_{1} = \frac{D_{1}}{θ_{1, m a x}}, {\tilde{D}}_{2} = \frac{D_{2}}{θ_{2, m a x}}

(4)

where

θ_{1, m a x}

and

θ_{2, m a x}

are the most severe limit state threshold for

D_{1}

and

D_{2}

, respectively. The scalar threshold associated with limit state k is then defined

θ_{J, k}

so that exceedance of approximates the joint event in equation (1) in a probabilistic form.

As illustrated in Figure 2, the parameters ${a_{k}, θ_{J, k}}$ are calibrated by minimizing the discrepancy between the joint exceedance probabilities and those derived from D_J across all IM levels and LS. The IM range is first partitioned into representative bins. Within each bin, empirical joint exceedance probabilities are estimated from equation (2), and compared with the exceedance probabilities obtained from the scalar damage index using equation (4). The scalar exceedance probabilities are obtained by fitting a parametric fragility model to the binary exceedance indicator defined by the proposed D_J. The calibration loss is defined as a weighted sum of squared differences between the joint and scalar exceedance probabilities by equation (5).

L ({a_{k}, θ_{J, k}}) = \sum_{k = 1}^{N} w_{k} {\sum_{i = 1}^{M} {[P_{joint, k} (D_{joint, i} \geq θ_{k} | I_{i}) - P_{J, k} (D_{J, i} \geq θ_{J, k} | I_{i})]}^{2}}^{1 / 2}

(5)

where

P_{joint}

and

P_{J}

denote empirical estimates of the joint and scalar exceedance probabilities, respectively. N is the number of limit state, M is the number of IM-EDP pairs, and

w_{k}

are user-specified importance weights which enable acceleration of convergence. The constrained optimization is solved with robust solvers (multi-start global search in MATLAB, bounds

0 \leq a_{k} \leq 1

Figure 2.

Development of a new damage index D_J that integrates both MD-based and RD-based damage index with principle of exceedance probability equivalence for each limit state.

In this manner, the resulting damage index D_J provides a compact surrogate that preserves the probabilistic characteristics of the original joint formulation while enabling a univariate fragility representation.

cGAN augmentation

Robust calibration of ${a_{k}, θ_{J, k}}$ and reliable characterization of $P_{J} (D_{J, i} \geq θ_{J, k} | I_{i})$ require adequate and unbiased sampling of the conditional distribution $P (D_{1}, D_{2} | I_{i})$ . When empirical or simulated datasets under-represent high IM or rare joint outcomes, parameter estimates, and fragility tails become unstable. A conditional generative adversarial network (cGAN) is employed to improve sampling of the joint IM-EDP space, particularly in regions associated with rare or extreme response levels, thereby enabling realistic generation of enough IM-EDP pairs for fragility calibration. Figure 3 illustrates the components of cGAN model, consisting of two competing neural networks: a generator (G) and a discriminator (D). The generator (G) produces synthetic EDP pairs conditioned on the input IMs (matrix Z and K), while the discriminator (D) evaluates whether the samples are consistent with the reference data r which is separate set of IM-EDP results obtained from simulations under recorded ground motion inputs, including variations in magnitude, source-to-site distance, fault mechanism, and near-fault pulse characteristics.

Figure 3.

Schematic of cGAN model.

The discriminator is trained by minimizing the minibatch loss:

L_{D} = - \frac{1}{M} \sum_{i = 1}^{M} [logD (r_{i}, k_{i}) + \log (1 - D (s_{i}, k_{i}))]

(6)

while the adversarial component of the generator loss is

L_{G} = - \frac{1}{M} \sum_{i = 1}^{M} logD (s_{i}, k_{i})

(7)

Adversarial training is driven by the minimax objective function:

\underset{G}{m i n} \underset{D}{m a x} E_{x_{real}, c} [\log D (x_{real}, c)] + E_{x_{real}, c} [\log (1 - D (G (z, c), c))]

(8)

where x_real denotes the true EDP data, c the conditioning IMs, and z ∼ N(0, I) the latent noise. Through adversarial training, the generator learns to reproduce the conditional distribution of response data given the input IMs. To ensure admissibility of the generated samples for fragility analysis, two additional constraints are incorporated into the generator objective. A physics penalty enforces simple, differentiable constraints C (⋅) (for the penalty leverages C (D₁) = −D₁ and C(D₂) = −D₂ to enforce the non-negativity of maximum and residual displacement, respectively, by generating a C > 0 loss only when the predictions D < 0):

L_{phy} = \frac{1}{M} \sum_{i = 1}^{M} ‖ C (s_{i}, k_{i}) ‖_{2}^{2}

(9)

A fragility-consistency penalty is introduced to align the exceedance frequencies computed from generated samples with the empirical joint exceedance probability ${\hat{P}}_{data, k} (D_{J, i} \geq θ_{J, k} | I_{i})$ within each IM bin. For a minibatch, synthetic exceedance probabilities are computed per IM bin IM_j and per limit state θ_J,k as:

{\hat{P}}_{gen, k} (D_{J, i} \geq θ_{J, k} | I_{i}) = \frac{1}{| B_{i} |} \sum_{i \in B} 1 {joint (s_{i}, θ_{J, k})}

(10)

where B_i indexes the number of the minibatch samples in bin I_i and 1{⋅} is the indicator of the joint

θ_{joint, k}

is exceeded by

D_{J, i}

. The fragility loss is therefore

L_{frag} = \sum_{k = 1}^{N} w_{k} {{\sum_{j = 1}^{M} ({\hat{P}}_{data, k} (D_{J, i} \geq θ_{J, k} | I_{i}) - {\hat{P}}_{gen, k} (D_{J, i} \geq θ_{J, k} | I_{i}))}^{2}}^{1 / 2}

(11)

The total generator objective loss minimized in practice is

L_{total} = L_{G} + λ_{phy} L_{phy} + λ_{frag} L_{frag}

(12)

which

L_{total}

matches the “G Total Loss” shown in Figure 3. This loss formulation enables the model to reconstruct the joint IM-EDP distribution in a manner consistent with both the underlying data and the target fragility behavior systematically. In addition, an independent set of IM-EDP results r (Figure 3) obtained from simulations under recorded ground motions is incorporated through the discriminator as a probabilistic reference, providing an external constraint on the learned distribution.

To ensure robust and stable model training, four major hyperparameters are needed to be systematically tuned (Gao et al., 2024; Wu et al., 2023): (1) Learning rate of the discriminator (lrD), which controls how aggressively the discriminator updates its parameters in response to classification errors. An excessively high lrD leads to discriminator dominance and generator collapse. (2) Learning rate of the generator (lrG), which determines the generator ability to adapt to the discriminator and loss feedback. (3) Noise dimension (noise dim), representing the length of the random latent vector z input to the generator, which governs the diversity of generated samples. (4) Batch size, defining the number of samples per training iteration. Smaller batches introduce high gradient variance. During training, the cGAN was optimized via alternating gradient updates using the Adam optimizer (Goodfellow et al., 2014; Marano et al., 2024). Spectral normalization and gradient penalties were imposed on D to prevent exploding gradients and to ensure Lipschitz continuity (Goodfellow et al., 2014; Marano et al., 2024). This training configuration ensured the convergence of both networks toward a Nash equilibrium (Marano et al., 2024; Ratliff et al., 2013), yielding physically consistent and statistically robust synthetic datasets that preserved the conditional dependence between IM and damage response.

Figure 4 depicts the iterative procedure used to train the cGAN model of Figure 3 and to convert the learned generator into a validated sampler for fragility calibration. Observational IM-EDP pairs are encoded in the condition matrix k and empirical joint exceedance frequencies ${\hat{P}}_{data, k} (D_{J, i} \geq θ_{J, k} | I_{i})$ are precomputed. Training proceeds in minibatches: for each minibatch a noise pair (Z _i , K _i ) is drawn, and the generator produces s_i. The discriminator is updated first using equation (6) and then the generator is updated by backpropagation of gradients of $L_{t o t a l}$ in equation (11).

Figure 4.

The flowchart of training with cGAN model.

Because minibatch estimates can be noisy for sparsely populated IM bins, stabilization is applied by maintaining a running-average estimator

{\bar{P}}_{gen}^{t + 1} (D_{J, i} \geq θ_{k} | I_{i}) = a {\bar{P}}_{gen}^{t} (D_{J, i} \geq θ_{k} | I_{i}) + (1 - a) {\hat{P}}_{gen}^{t} (D_{J, i} \geq θ_{k} | I_{i})

(13)

chosen to balance recency and variance reduction. This running average is used in L_frag to reduce gradient variance and avoid destabilizing updates from extreme minibatches. Model validation is performed periodically on the matrix V shown in Figure 4, and checkpoints are saved and ranked using a composite validation score that balances marginal realism, quantified by Wasserstein distance and KS statistics, against fragility fidelity, measured by the root-mean-square difference between empirical and synthetic exceedance frequencies over IM bins and LS:

{RSME}_{frag} = \frac{1}{N \cdot M} \sum_{k = 1}^{N} w_{k} {{\sum_{j = 1}^{M} ({\hat{P}}_{data, k} (D_{J, i} \geq θ_{J, k} | I_{i}) - {\hat{P}}_{gen, k} (D_{J, i} \geq θ_{J, k} | I_{i}))}^{2}}^{1 / 2}

(14)

The combined adversarial loss plus physics/fragility loss supervision therefore produces variance-reducing, engineering-consistent synthetic data suitable for robust calibration and validation. Note that the Wasserstein distance and KS statistics are used only for evaluation of marginal distributions and checkpoint selection. The cGAN is trained with the standard Binary Cross-Entropy (BCE)-based adversarial loss, rather than a Wasserstein GAN (WGAN) objective.

The above framework is then applied to a representative ordinary standard highway bridge to demonstrate practical implementation and to achieve the cGAN augmentation on calibration stability, fragility curves. Modeling details and numerical results are reported in the case study in “Case study of ordinary concrete highway bridge”. It should be noted that the proposed damage-index formulation, based on probabilistic equivalence to joint exceedance behavior, is general and can be applied to different structural systems and response measures. The calibrated quantities, including weighting coefficients and limit-state thresholds, are obtained from the underlying modeling dataset and therefore depend on the structural model, selected response measures, ground-motion ensemble, and site conditions considered in the analysis. Accordingly, the parameter values reported in the case study should be interpreted as representative of the specific configuration analyzed. For other bridge-column typologies or modeling assumptions, the same framework can be applied to obtain a corresponding parameter set for fragility analysis in risk assessment.

Case study of ordinary concrete highway bridge

Prototype bridge model

The prototype of a conventional two-span cast-in-place/post-tensioned (CIP/PT) box-girder highway bridge is presented in Figure 5, which is sited in downtown Los Angeles. Each span measures 45.7 m and the superstructure is supported at midspan by a single column bent. The bent column is circular (1.2 m diameter), contains 18 No.29 (US #9) longitudinal bars and No.25 (US #8) transverse hoops at 0.165 m spacing, and extends 9.1 m from the bent-cap centroid to the column base. Accounting for the tributary superstructure mass, the first natural period of the system is approximately 2.0 s. Foundation soils are characterized by a site average shear-wave velocity over the upper 30 m V_s30 of 270 m/s.

Figure 5.

Illustration of the numerical box-girder highway bridge model for the case study.

The column was modeled in OpenSees using fiber elements to capture flexural inelasticity and energy dissipation. The section was discretized into confined and unconfined concrete fibers and longitudinal steel layers, with cyclic constitutive laws calibrated for large-strain RC behavior. Flexural yielding was concentrated in end regions modeled with “Force-Beam-Column” elements based on the iterative force-based formulation, which represents both distributed and concentrated plasticity. Numerical integration employed a modified two-point Gauss-Radau scheme, placing integration points at the element ends and at 8/3 of the hinge length inside the element, ensuring exact representation of linear curvature and consistent softening length (Scott, 2011). A plastic hinge length of 1.0 m, determined from the empirical expressions in (Caltrans, 2025), was assigned to the column ends. Rigid links enforced kinematic compatibility between top of the bent column and girder. Geometric nonlinearity (P-Δ effects) was included. Ancillary components such as shear keys and abutment backfill were represented by elastic elements connected to the girder and restrained by the ground, as illustrated in Figure 5. This simplified representation was adopted to capture their boundary restraint effects while isolating the dominant transverse column response and avoiding additional system-level uncertainties in demonstrating the proposed fragility-informed cGAN framework. Rayleigh damping was applied to represent modal energy dissipation consisting with RC bridge practice. Following the approach of Yoon et al. (2019), a monotonic pushover analysis of the fiber-section model yielded a yield displacement $Δ_{Y}$ of 0.14 m in Figure 6. The ultimate displacement $Δ_{U}$ is obtained from a calibrated confined-concrete model and represents an experimentally consistent upper-bound deformation (e.g., Limit state 4 in Figure 2) capacity (Yoon et al., 2019). These quantities are used as reference capacities for normalization in the subsequent fragility and damage-index calculations.

Figure 6.

Force-displacement curve from pushover analysis.

Ground-motion dataset for training

Seismic inputs comprised both recorded pulse-like and nonpulse-like ground motions, alongside site-specific synthetic realizations. A total of 488 recorded pulse-like motions, each with two orthogonal horizontal components and identified by Shahi and Baker (2014) via an efficient MATLAB algorithm from the NGA-West2 database, were adopted. Besides these, 1000 synthetic ground motions were generated following the MATLAB procedure of Yamamoto and Baker (2013), based on the fault and site parameters summarized in Figure 7, ensuring broad coverage of intensity and pulse characteristics relevant to near-fault and long-period bridge response.

Figure 7.

Information of the fault at the bridge site.

Additionally, a set of 100 recorded ground motions (two orthogonal horizontal components each) selected from the NGA database following Mackie and Stojadinovic (2005) is incorporated to provide an independent probabilistic reference. The dataset spans a wide range of shaking intensities (PGA from 0 to 0.80 g) and includes diverse seismic characteristics, such as variations in magnitude, source-to-site distance, fault mechanism, and near-fault pulse effects. Rather than being used directly in fragility construction, this dataset is employed during the adversarial training process, where the corresponding IM-EDP results obtained from OpenSees simulations serve as a reference distribution for the discriminator. This setup enables the cGAN to align the generated samples with an independently simulated probabilistic representation of structural response, particularly in sparsely sampled regions of the response space. As a result, the generated data support a more stable reconstruction of the joint IM-EDP distribution and associated fragility behavior, especially in the moderate-to-severe damage regime.

Although two orthogonal horizontal components are available for each recorded motion, the nonlinear time-history analyses were performed using unidirectional input applied in the transverse direction of the bridge, consistent with the focus on transverse column response. This simplification allows the dominant damage-related response to be isolated while avoiding additional system-level complexities associated with multi-directional excitation.

The peak ground velocity (PGV) was chosen as the primary intensity measure due to its strong correlation with combined peak and residual column displacements, while alternative IMs (e.g., PGA, spectral acceleration (Sa)) were assessed in following studies. Nonlinear time-history analyses of the ground-motion suite yielded maximum and residual transverse displacements at the column top, with post-processing removing nonphysical responses and imposing a numerical floor to stabilize log-space fragility fits. The resulting (IM, EDP) pairs, together with held-out validation subsets, form the empirical database for cGAN training, synthetic ensemble generation, and calibration of the proposed D_J and optimal θ_J.

Results and discussion

The selection of an appropriate IM is an important aspect of fragility analysis, as it directly influences the characterization of seismic demand and the resulting fragility estimates. Previous studies have shown that different IM choices can lead to notable variations in fragility relationships, particularly in terms of dispersion and predictive capability, and that single scalar IMs may be insufficient to capture complex response features (Giannini et al., 2022).

In this study, three commonly used scalar IMs, namely PGA, PGV, and Sa, are examined using the dataset obtained from nonlinear time-history analyses (NTHAs) of the prototype bridge in Figure 8. Their relationships with maximum and residual displacement are evaluated through linear regression after normalizing the IMs to the range [0,1].

Figure 8.

Correlation between IMs and maximum-residual displacement based on ground-motion dataset for training.

The results indicate that all three IMs exhibit limited predictive capability for residual displacement, as reflected by relatively low R² values. This is consistent with the inherently dispersed nature of residual deformation, which depends on cumulative and path-dependent structural behavior and is therefore not well captured by a single scalar IM. Among the candidates considered, PGV shows comparatively stronger and more consistent correlation with both maximum and residual displacement, while PGA and Sa exhibit weaker associations for this benchmark case.

Accordingly, PGV is adopted as the primary IM in this study as a comparative and case-specific choice, rather than as a universally optimal predictor. Within the proposed framework, the IM serves primarily as a conditioning variable for representing the joint IM-EDP distribution, rather than as a deterministic predictor of structural response. In this context, PGV provides a suitable balance for conditioning the joint MD-RD response space and is therefore used for subsequent fragility calibration and cGAN training.

To ensure robust and stable model training, the sensitivity of the cGAN to four major hyperparameters (lrD, lrG, noise dimension, and batch size) is evaluated through a bootstrap resampling analysis (Reiff et al., 2016) of ln (PGV) when the correspond D_J exceeding the θ₁ ∼ θ₄. During bootstrap analysis, the robustness and stability of the training process are quantified based on the mean $μ$ and standard deviation $σ$ of ln(PGV) obtained from repeated resampling of the training data, with their distributional trends illustrated by box plots to characterize variability across bootstrap realizations in Figures 9 to 12. As shown in Figures 9 and 10, the generator (G) exhibits higher sensitivity to the learning rate, with its bootstrap results of $μ$ gradually increasing as the learning rate decreases beyond lrG = $6 \times 10^{- 5}$ . In contrast, the discriminator (D) remains generally stable over the range of IrD = $8 \times 10^{- 5}$ ∼ $1.2 \times 10^{- 3}$ . The main reason is that G loss combines both adversary and fragility terms (Figure 3), which generate gradients that can point in different directions. As a result, changes in lrG strongly affect how G balances producing realistic outputs while satisfying the fragility constraints. In contrast, the discriminator (D) only optimizes a simple adversarial D loss (Figure 3), so its behavior remains relatively stable across different learning rates of D. Accordingly, lrG = $6 \times 10^{- 5}$ and IrD = $1 \times 10^{- 4}$ were selected for subsequent hyperparameter sensitivity analyses.

Figure 9.

Bootstrap analysis of the mean $μ$ and standard deviation $σ$ for the G learning rate (lrG).

Figure 10.

Bootstrap analysis of the mean $μ$ and standard deviation $σ$ for the D learning rate (lrD).

Figure 11.

Bootstrap analysis of the mean $μ$ and standard deviation $σ$ for dimension of the input noise vector (noise dim).

Figure 12.

Bootstrap analysis of the mean $μ$ and standard deviation $σ$ for the number of samples per batch during training (batch size).

As demonstrated in Figures 11 and 12, the $μ$ and $σ$ of ln(PGV) obtained from bootstrap resampling remain largely unchanged across different noise dimensions (noise dim) and batch sizes, indicating that these hyperparameters have minimal influence on the stability and robustness of cGAN training. The observed phenomenon above is consistent with cGAN characteristics, where the conditional seismic inputs primarily govern the generated structural responses, and the stochastic noise vector contributes only minor variability. The chosen batch size provides reasonably stable gradient estimates during training, and the network architecture together with the available training data allows the generator to capture the main features of the underlying response distribution. Under these settings, the model training exhibits stable convergence and produces generated data that aligns well with the empirical distribution. In the following, lrG = $6 \times 10^{- 5}$ , lrD = $1 \times 10^{- 4}$ , noise dim = 24 and batch size = 300, are used in the subsequent cGAN training.

Using these parameters, the cGAN models are trained with a composite loss combining adversarial, physics, and fragility components. Figure 13 illustrates the smooth evolution of G and D losses, physical loss, and fragility loss as training progresses. Early fluctuations gradually diminish, and both physics and fragility penalties stabilize at low values, confirming that the generator has learned to produce physically consistent conditional responses.

Figure 13.

Evolution of generator and discriminator losses during cGAN training. The plot shows adversarial loss (G Loss), discriminator loss (D Loss), physical loss, and fragility loss versus training epochs.

Figures 14 and 15 indicate that at epochs 100 and 1000, the generated data underrepresent extreme responses, and the joint probability density distributions are incomplete. With training extended to epoch 10,000, extreme values are well captured, and the joint probability density distributions closely match the empirical distribution, confirming stable convergence and accurate reproduction of the relationship between IM and EDP.

Figure 14.

Scatter plots of generated versus real EDPs (left: maximum displacement; right: residual displacement) with respect to PGV. (a) Results at initial training epoch = 100; (b) results at final training epoch = 1000; and (c) results at final training epoch = 10,000.

Figure 15.

Joint probability density distributions of generated versus real EDPs (left: maximum displacement; right: residual displacement) with respect to PGV. (a) Results at initial training epoch = 100; (b) results at final training epoch = 1000; and (c): results at final training epoch = 10,000.

As shown in Figure 16, the bootstrap-derived $μ$ and $σ$ of the generated ln(PGV) for D_J exceeding LS₁ ∼ LS₄ shrink steadily with increasing training epochs after epoch of 10k, reflecting enhanced stability and robustness of the cGAN outputs. Correspondingly, Figure 17 demonstrates that, as training progresses, the fragility curves for LS₁ ∼ LS₄ transition from initially biased forms (e.g., epoch 1000) to curves that closely approximate those obtained from the original input data (e.g., epoch 10,000).Among epoch of 1000 ∼ 15,000, the fragility curve at epoch of 10,000 shows the best agreement with the input data and is thus chosen for further analyses.

Figure 16.

Bootstrap analysis of the mean $μ$ and standard deviation $σ$ for the round of training (epoch).

Figure 17.

Fragility curves with 95% confidence intervals for different cGAN-generated epochs (LS₁ ∼ LS₄).

The evolution of the residual-demand weighting factor $a_{k}$ and the corresponding D_J limit state thresholds θ_J,k is summarized in Tables 2 and 3. With a limited number of synthetic samples (e.g., 1k), noticeable deviations in $a_{k}$ are observed relative to the OpenSees reference, particularly for higher damage states. As the amount of generated data increases, the weights for LS₁ ∼ LS₂ tend to stabilize near 0.1, while the weight for LS₄ shows greater variability, reflecting the sensitivity of high-damage calibration to the limited availability of extreme-response samples. A similar trend is observed for the calibrated thresholds θ_J,k. With increasing data volume, the estimates become more stable overall, although greater variability persists in LS₃ ∼ LS₄ due to the sparse sampling in the tail region. For larger datasets, the calibrated values approach those obtained from OpenSees (e.g., LS₃ and LS₄ are approximately 0.788 and 0.99, compared to 0.725 and 0.975), indicating consistency at the level of fragility behavior rather than exact parameter agreement. These results suggest that increasing the amount of augmented data improves the stability of joint calibration.

Table 2.

Residual demand weight $a_{k}$ at θ_J,k versus cGAN data volume.

Data source	$a_{1}$	$a_{2}$	$a_{3}$	$a_{4}$
OpenSees	0.100	0.131	0.725	0.600
cGAN(1k)	0.163	0.225	0.350	0.100
cGAN(5k)	0.108	0.120	0.303	0.600
cGAN(10k)	0.102	0.100	0.303	0.600
cGAN(20k)	0.102	0.100	0.319	0.975
cGAN(50k)	0.102	0.100	0.317	0.975

Table 3.

Optimal θ_J,k thresholds of D_J versus cGAN data volume.

Data source	θ _J,1	θ _J,2	θ _J,3	θ _J,4
OpenSees	0.342	0.584	0.725	0.975
cGAN(1k)	0.350	0.498	0.788	0.975
cGAN(5k)	0.342	0.569	0.788	1.006
cGAN(10k)	0.346	0.567	0.789	0.991
cGAN(20k)	0.345	0.564	0.788	0.998
cGAN(50k)	0.345	0.579	0.788	0.998

The influence of cGAN-generated sample volume on fragility estimation is shown in Figures 18 and 19. As the generated IM-EDP sample size increases from 1k to 50k, the fragility curves become close to the results from input data and their 95% confidence intervals shrink, indicating reduced statistical uncertainty and more stable estimates. Although the narrower confidence intervals observed with larger synthetic datasets are primarily a result of increased sample size, this does not undermine the value of the cGAN-based data expansion. The key advantage lies in the model's ability to reproduce the underlying joint distribution of seismic IM-EDP relationships, thereby enabling more statistically stable and physically consistent fragility estimation under data-scarce conditions. Figures 18 and 19 show that with a small generated sample (1k), fragility curves deviate substantially from the input-data-based results. Increasing the sample size to 50k leads to convergence of $μ$ and $σ$ , with fragility curves increasingly matching the empirical results. The stabilization effect is particularly notable for high-limit states (LS₃ and LS₄), demonstrating that large-scale data augmentation is essential for reliable fragility estimation when real data are sparse.

Figure 18.

Bootstrap analysis of the mean $μ$ and standard deviation $σ$ for generated data volume of the cGAN.

Figure 19.

Fragility curves with 95% confidence intervals for different cGAN-generated data volumes (LS₁ ∼ LS₄).

Figure 20 compares the fragility curves derived from the OpenSees simulation dataset with those obtained from the cGAN-generated dataset. For LS₁ ∼ LS₃, the median exceedance probabilities are in close agreement, indicating that the cGAN preserves the central fragility behavior implied by the simulation data.

Figure 20.

Fragility curves with 95% confidence intervals based on the proposed D_J for LS₁ ∼ LS₄: comparison between OpenSees modeling and cGAN training. (with lrG = $6 \times 10^{- 5}$ , lrD = $1 \times 10^{- 4}$ , noise dim = 24, batch size = 300, epoch =10,000, and data volume = 20k).

In the highest damage state (LS₄), the cGAN-based fragility curve is observed to be lower than that obtained from the OpenSees dataset. This difference reflects the nature of the underlying estimates. The OpenSees-based fragility curve is derived from a finite set of simulations and, in the severe-damage range, is governed by a limited number of extreme-response samples, making the tail estimate sensitive to sampling variability.

In contrast, the cGAN is trained to reproduce the joint IM-EDP distribution under fragility constraints. During training, the discriminator evaluates the generated samples against an independent set of IM-EDP probabilistic results obtained from ground-motion-based simulations, which serve as a statistical reference for distributional consistency. This reference is not used in constructing the OpenSees fragility curves in Figure 20, but is incorporated through the adversarial training process.

The resulting LS₄ behavior therefore reflects a distribution-level reconstruction of the sparse tail informed by both the original simulation data and the independent statistical reference, rather than direct inheritance of fluctuations associated with limited samples. From an engineering perspective, this leads to a more stable characterization of extreme-damage probabilities for fragility-based assessment and decision-making. The role of the cGAN is therefore not to match simulation results pointwise, but to provide a more robust representation of exceedance behavior under limited-data conditions.

Figure 21 further compares D_J-based fragility curves with those obtained from the conventional “OR” and “AND” rules applied separately to MD and RD. The D_J curves consistently fall between the two extremes, reducing the over-conservatism of the “OR” rule and the under-conservatism of the “AND” rule. The weighted-sum formulation successfully integrates peak and residual damage effects, resulting in a single, physically interpretable index that is statistically reliable for seismic damage assessment.

Figure 21.

Comparing fragility curves based on the proposed D_J with exceedance probabilities calculated by “OR” and “and” rules for LS₁ ∼ LS₄.

Conclusions

This study develops a unified damage index D_J assisted with a fragility-informed cGAN model for seismic damage assessment of reinforced-concrete columns of bridges. The proposed D_J simultaneously captures both peak damage and post-earthquake responses by MD-based and RD-based damage indexes into a single metric. To address data scarcity and improve fragility reliability, a cGAN model incorporating a fragility-loss penalty is introduced for statistic-consistent data augmentation. The major conclusions are drawn as follows:

A probabilistic formulation is established to unify MD-based and RD-based damage index through a weighted joint exceedance rule that balances peak damage and postearthquake functionality considerations. A fragility-equivalent D_J with the probabilistic formulation is then achieved via a MATLAB-based optimization procedure, which determines the residual-demand weights $a_{k}$ and limit states thresholds θ_J,k by enforcing joint probability equivalence.

A fragility-informed cGAN is developed to augment IM-EDP data under limited-data conditions within a probabilistically constrained framework. The model integrates adversarial, physics-based, and fragility-consistency loss terms, enabling the generated samples to preserve exceedance behavior across limit states. Compared with fragility estimation based solely on a finite set of OpenSees simulations, this formulation reconstructs the joint IM-EDP behavior at the distribution level, reducing sensitivity to individual extreme-response observations and providing a more stable basis for fragility estimation, particularly in the sparsely sampled severe-damage regime.

Sensitivity analysis reveals that the generator learning rate (lrG) predominantly governs training stability, whereas discriminator rate, noise dimension, and batch size have minor effects. With the optimized lrG = $6 \times 10^{- 5}$ , the composite adversarial-physics-fragility loss converges smoothly, yielding synthetic IM-EDP data that closely match empirical distributions.

The inclusion of the independent IM-EDP dataset obtained from simulations under recorded ground motions introduces an external probabilistic reference into the training process, enabling distribution-level constraint beyond the original simulation data. This design mitigates the closed-loop nature of single-dataset training and supports a more stable reconstruction of exceedance behavior, particularly in sparsely sampled high-damage regimes.

The cGAN-based augmentation improves the conditioning of fragility estimation by increasing the effective sample size in underrepresented regions of the response space. With 10k or more generated samples, a reduction in dispersion of ln (PGV) and a narrowing of fragility confidence intervals are observed, particularly for higher damage states (LS₃ ∼ LS₄). These effects reflect a more stable representation of exceedance behavior at the distribution level, supporting more reliable fragility-based decision-making under limited-data conditions, where estimates from sparse simulation samples may be highly variable.

As the volume of augmented data increases, the residual-demand weights $a_{k}$ and calibrated thresholds θ_J,k of the proposed damage index D_J exhibit convergence and reduced variability. This behavior indicates improved stability of the joint probability calibration across multiple damage states, while remaining conditioned on the underlying simulation framework.

The proposed framework provides a general approach for constructing a unified damage index based on probabilistic equivalence to joint response behavior. The calibrated parameters, including weighting factors and limit-state thresholds, depend on the structural model, response measures, and ground-motion dataset, and should be re-identified for other bridge-column configurations. Extending the framework to a wider range of structural systems and hazard conditions remains an important direction for future work.

Footnotes

ORCID iD

Lizhi Sun

CRediT authorship contribution statement

JL: writing—original draft, validation, software, methodology, formal analysis, and data curation. JF: writing—review and editing, validation, and resources. LS: writing—review and editing, supervision, project administration, methodology, funding acquisition, and conceptualization.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Engineering Partners International, LLC, (Grant Number PC13803).

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Data availability statement

All data generated or used during the study are available from the corresponding author by reasonable request.

References

ASCE (2022) Minimum design loads and associated criteria for buildings and other structures. American Society of Civil Engineers, Reston, VA 20191, USA.

Caltrans (2025) Seismic Design Criteria (SDC) – Version 2.1. California Department of Transportation, Sacramento, CA, USA.

Choi

DesRoches

Nielson

(2004) Seismic fragility of typical bridges in moderate seismic zones. Engineering Structures 26(2): 187–199.

Cornell

Jalayer

Hamburger

, et al. (2002) Probabilistic basis for 2000 SAC federal emergency management agency steel moment frame guidelines. Journal of Structural Engineering 128(4). 4. American Society of Civil Engineers: 526–533.

Padgett

Shafieezadeh

(2019) A posteriori optimal intensity measures for probabilistic seismic demand modeling. Bulletin of Earthquake Engineering 17(2): 681–706.

Wang

Xie

, et al. (2023) Regional seismic risk and resilience assessment: Methodological development, applicability, and future research needs – An earthquake engineering perspective. Reliability Engineering & System Safety 233: 109104.

Gan

, et al. (2019) Identification of ground intrusion in underground structures based on distributed structural vibration detected by ultra-weak FBG sensing technology. Sensors 19(9): 2160.

Gao

Zhang

, et al. (2024) Data augmentation using improved conditional GAN under extremely limited fault samples and its application in fault diagnosis of electric submersible pump. Journal of the Franklin Institute 361(4). Elsevier: 106629.

Ghobarah

(2001) Performance-based design in earthquake engineering: State of development. Engineering Structures 23(8): 878–884.

10.

Giannini

Paolacci

Phan

, et al. (2022) A novel framework for seismic risk assessment of structures. Earthquake Engineering & Structural Dynamics 51(14). John Wiley & Sons, Ltd: 3416–3433.

11.

Goodfellow

Pouget-Abadie

Mirza

, et al. (2014) Generative adversarial nets. In: Proceedings of the 28th International Conference on Neural Information Processing Systems 2: 2672–2680.

12.

Housner

Thiel

Jr (1995) The continuing challenge: Report on the performance of state bridges in the Northridge earthquake. Earthquake Spectra 11(4): 607–636.

13.

Kang

Zhou

Kumar

, et al. (2025) Deep learning method for post-earthquake damage assessment of frame structures based on time–frequency analysis and CGAN. Earth Systems and Environment 9(1): 403–420.

14.

Kim

T-H

Lee

K-M

Chung

Y-S

, et al. (2005) Seismic damage assessment of reinforced concrete bridge columns. Engineering Structures 27(4): 576–592.

15.

Sun

(2020) Detectability of bridge-structural damage based on fiber optic sensing through deep-convolutional neural networks. Journal of Bridge Engineering 25(4): 04020012.

16.

Liu

Sun

(2005) Multiscale modeling of elastoplastic deformation and strengthening mechanisms in aluminum-based amorphous nanocomposites. Acta Materialia 53(9): 2693–2701.

17.

Liu

Zhou

Zhao

, et al. (2026) Physics-enhanced machine learning surrogate modelling for real-time seismic fragility assessment of reinforced concrete bridges. Engineering Structures 348: 121772.

18.

Mackie

Stojadinovic

(2004) Residual displacement and post-earthquake capacity of highway bridges. In: Proceedings of the 13th World Conference on Earthquake Engineering, Vancouver, British Columbia, Canada, 1–6 August, 13 WCEE Secretariat, Paper No. 1550, pp. 1–16.

19.

Mackie

Stojadinovic

(2005) Fragility basis for California highway overpass bridge seismic decision making. Report No. 2005/02, Pacific Earthquake Engineering Research Center, Berkeley, CA, USA.

20.

Marano

Rosso

Aloisio

, et al. (2024) Generative adversarial networks review in earthquake-related engineering fields. Bulletin of Earthquake Engineering 22(7): 3511–3562.

21.

Nguyen

Phan

Pham

, et al. (2026) Portfolio-scale seismic fragility of RC bridge columns with series-distributed neural networks. International Journal of Disaster Risk Reduction 132: 105955.

22.

Ning

Xie

Burton

, et al. (2024) Enabling efficient regional seismic fragility assessment of multi-component bridge portfolios through Gaussian process regression and active learning. Earthquake Engineering & Structural Dynamics 53(9): 2929–2949.

23.

Ratliff

Burden

Sastry

(2013) Characterization and computation of local nash equilibria in continuous games. In: 2013 51st Annual Allerton Conference on Communication, Control, and Computing (Allerton), pp.917–924. IEEE.

24.

Reiff

Sanayei

Vogel

(2016) Statistical bridge damage detection using girder distribution factors. Engineering Structures 109: 139–151.

25.

Saini

Saiidi

(2014) Probabilistic damage control approach for seismic design of bridge columns. California Department of Transportation Report No. CA14–2280, Sacramento, CA, USA.

26.

Scott

(2011) Numerical integration options for the force-based beam–column element in OpenSees. Technical report, Pacific Earthquake Engineering Research Center, Berkeley, CA 94720, USA.

27.

Seyed Ardakani

Saiidi

(2018) Simple method to estimate residual displacement in concrete bridge columns under near-fault earthquake motions. Engineering Structures 176: 208–219.

28.

Shahi

Baker

(2014) An efficient algorithm to identify strong-velocity pulses in multicomponent ground motions. Bulletin of the Seismological Society of America 104(5): 2456–2466.

29.

Todorov

Billah

AHMM

(2021) Seismic fragility and damage assessment of reinforced concrete bridge pier under long-duration, near-fault, and far-field ground motions. Structures 31: 671–685.

30.

Wei

Jia

Bai

, et al. (2023) Seismic resilience assessment of bridges considering both maximum and residual displacements. Engineering Structures 291: 116420.

31.

T-Y

R-T

Wang

P-H

, et al. (2023) Development of a high-fidelity failure prediction system for reinforced concrete bridge columns using generative adversarial networks. Engineering Structures 286: 116130.

32.

Xiang

Chen

Alam

(2020) Probabilistic seismic fragility and loss analysis of concrete bridge piers with superelastic shape memory alloy-steel coupled reinforcing bars. Engineering Structures 207: 110229.

33.

Yamamoto

Baker

(2013) Stochastic model for earthquake ground motion using wavelet packets. Bulletin of the Seismological Society of America 103(6). Seismological Society of America: 3044–3056.

34.

Yazdani

Yazdannejad

(2019) Estimation of the seismic demand model for different damage levels. Engineering Structures 194: 183–195.

35.

Yin

Sun

(2005) Elastic modeling of periodic composites with particle interactions. Philosophical Magazine Letters 85(4): 163–173.

36.

Yoon

Ataya

Mahan

, et al. (2019) Probabilistic damage control application: Implementation of performance-based earthquake engineering in seismic design of highway bridge columns. Journal of Bridge Engineering 24(7): 04019068.

37.

Yoon

Ataya

Mahan

, et al. (2022) Total risk model and loss analysis of probabilistic damage control application in seismic design of highway bridges. Journal of Bridge Engineering 27(5): 04022020.

38.

Zameeruddin

Sangle

(2016) Review on recent developments in the performance-based seismic design of reinforced concrete structures. Structures 6: 119–133.

39.

Zhong

Zhu

Mangalathu

, et al. (2023) Probabilistic curvature-and-drift limit states predictive models of high-strength bridge columns. Engineering Structures 295: 116859.