Bridge Digital Twinning Using an Output-Only Bayesian Model Updating Method and Recorded Seismic Measurements

Rapid post-earthquake damage diagnosis of bridges can guide decision-making for emergency response management and recovery. This can be facilitated using digital technologies to remove the barriers of manual post-event inspections. Prior mechanics-based Finite Element (FE) models can be used for post-event response simulation using the measured ground motions at nearby stations; however, the damage assessment outcomes would suffer from uncertainties in structural and soil material properties, input excitations, etc. For instrumented bridges, these uncertainties can be reduced by integrating sensory data with prior models through a model updating approach. This study presents a sequential Bayesian model updating technique, through which a linear/nonlinear FE model, including soil-structure interaction effects, and the foundation input motions are jointly identified from measured acceleration responses. The efficacy of the presented model updating technique is first examined through a numerical verification study. Then, seismic data recorded from the San Rogue Canyon Bridge in California are used for a real-world case study. Comparison between the free-field and the foundation input motions reveals valuable information regarding the soil-structure interaction effects at the bridge site. Moreover, the reasonable agreement between the recorded and estimated bridge responses shows the potentials of the presented model updating technique for real-world applications. The described process is a practice of digital twinning and the updated FE model is considered as the digital twin of the bridge and can be used to analyze the bridge and monitor the structural response at element, section, and fiber levels to diagnose the location and severity of any potential damage mechanism.


Introduction
The field of computational structural mechanics has advanced to a mature level to facilitate high-fidelity and computationally-efficient seismic response simulation of bridges [1][2][3]. Practitioners and researchers use mechanics-based Finite Element (FE) models for response prediction of complex bridge structures. Nevertheless, these models include inherent uncertainties when it comes to mirroring real-world response behavior. These can include uncertainties in soil and structural material models and parameters, dynamic input excitations, uncertainties in damping energy dissipation mechanisms, etc. The uncertainties can be quantified and/or reduced by integrating mechanics-based FE models with measured responses of the structure through a data assimilation approach. This approach consists of training/updating models with measurements to estimate uncertain/unknown model parameters [4] (i.e., input-output model updating) or uncertain/unknown model and input parameters [5,6] (i.e., output-only model updating). The trained/updated model In this paper, to transcend the aforementioned technical challenges in the application of digital twinning and virtual sensing, an output-only model updating technique in the time domain is presented. Through this presented technique, the uncertain parameters of a linear or nonlinear mechanics-based FE model along with the FIMs can be estimated using sparsely measured acceleration responses recorded during an earthquake. The model updating technique presented herein is mainly based on a sequential Bayesian inference algorithm originally developed in previous works [5,21]. The efficacy of the presented model updating technique is verified through numerically simulated data using the San Roque Canyon (SRC) bridge model as a testbed. Then, the seismic data collected from the SRC bridge are used for real-world case studies. Although, the amplitudes of the available recorded motions are low, and thus the bridge mainly behaves in the linear-elastic regime, the process can be readily applied to strong earthquakes and nonlinear structural behavior. The main novelty of this paper is the application of the time-domain output-only model updating technique in real-world settings to investigate its efficacy and limitations.
The outline of the paper is as follows. First, the sequential Bayesian inference for output-only model updating technique and the identifiability analysis approach to determine the identifiable model parameters are briefly presented in Section 2. Section 3 is focused on the verification of the presented model updating technique using numerically simulated data. The FE model of SRC bridge is presented in detail in Section 3.1, and the identifiability analysis is performed in Section 3.2. Following that, the results of the verification study in a numerically simulated environment are presented in Section 3.3. Five case studies are carried out in Section 4 using real-world earthquake datasets and the results are discussed. In Section 5, the application of DT and virtual sensing is presented. Finally, the conclusion and future steps are discussed. study in a numerically simulated environment are presented in Section 3.3. Five case studies are carried out in Section 4 using real-world earthquake datasets and the results are discussed. In Section 5, the application of DT and virtual sensing is presented. Finally, the conclusion and future steps are discussed.

Sequential Bayesian Inference Method and Identifiability Analysis
The model updating technique for joint system and input estimation is based on a sequential Bayesian inference method using the unscented transformation approach [22] for uncertainty propagation. This technique will be reviewed briefly in Section 2.1. The Identifiability analysis process to quantify the information content of the measurement data and potential identifiability of model parameters is presented in Section 2.2.

Sequential Bayesian Inference Method for Output-Only FE Model Updating
The model updating technique for joint system and input estimation is schematically shown in Figure 1. This technique is closely similar to the works presented in [5,21], with some tweaks and improvements as will be outlined here. As shown in Figure 1, the unknown model parameters and inputs (here FIMs) are augmented into the unknown parameter vector , the uncertainties of which are expressed with a Gaussian Probability Density Function (PDF). These uncertainties are propagated into the FE model = ( ), in which (… ) is the nonlinear response function of the FE model. Next, a simulation (or prediction) error model ( ) is defined to correlate the FE-predicted response ( ) with the measured response ( ) collected by the sensors. Finally, the Bayes' theorem is used to find the posterior PDF of the unknown parameters, which is then used as the prior PDF for the next sequence of measured responses. A sequential estimation window approach is used in this study to improve the efficacy of the estimation process. In this approach, the time domain is divided into ≥ 2 successive overlapping estimation windows. The estimation window spans from time step to time step , and = − is the length of the estimation window. Also, the overlap between the and ( + 1) estimation windows is defined as , = − , ∀ ≤ − 1. The estimation problem is solved at each estimation window iteratively to estimate the mean vector and covariance matrix of the unknown parameter vector and then moves to the next estimation window until completion. The estimates of the unknown model parameters at the end of each estimation window are transferred to the next estimation window and used as initial estimates. However, to A sequential estimation window approach is used in this study to improve the efficacy of the estimation process. In this approach, the time domain is divided into n w ≥ 2 successive overlapping estimation windows. The w th estimation window spans from time step t w 1 to time step t w 2 , and t w l = t w 2 − t w 1 is the length of the w th estimation window. Also, the overlap between the w th and (w + 1) th estimation windows is defined as t w,w+1 The estimation problem is solved at each estimation window iteratively to estimate the mean vector and covariance matrix of the unknown parameter vector and then moves to the next estimation window until completion. The estimates of the unknown model parameters at the end of each estimation window are transferred to the next estimation window and used as initial estimates. However, to transfer the estimates of FIMs, each estimation window is divided into two parts. The estimates of FIMs in the first part, spans from t w 1 to t w+1 1 in the w th estimation window, which does not overlap with the next estimation window, are considered as final estimates. The second part overlaps with the next estimation window, spans from t w+1 1 to t w 2 in the w th estimation window, and the estimates of FIMs in this part are transferred to the next estimation window to be considered as initial estimates. The sequential estimation window approach is schematically shown in Figure 2 and further discussed in the following section.
transfer the estimates of FIMs, each estimation window is divided into two parts. The estimates of FIMs in the first part, spans from to in the estimation window, which does not overlap with the next estimation window, are considered as final estimates. The second part overlaps with the next estimation window, spans from to in the estimation window, and the estimates of FIMs in this part are transferred to the next estimation window to be considered as initial estimates. The sequential estimation window approach is schematically shown in Figure 2 and further discussed in the following section.
In Equation (1) in which : ∈ ℝ ( × )× is the simulation error vector at the estimation window and accounts for the misfit between the measured and FE-predicted responses of the bridge [23]. By neglecting the effects of modeling error, the simulation error at each time step is ideally modeled as a spatially and temporally independent zero-mean Gaussian white noise process (i.e., ~( , ) ). Hence, is a block diagonal matrix whose block diagonals are the simulation error covariance matrix . The objective of the model updating is to find the estimates of The FE-predicted response of a bridge at the w th estimation window, , to a general case of multiple-support earthquake excitation can be expressed as a nonlinear function of the model parameter vector, θ ∈ R n θ ×1 , and the time history of the multiple-support FIMs, . This is shown in Equation (1). The terms n θ , n s and n y are, respectively, the number of model parameters, the number of supports with different translational input motions, and the number of measurement channels.
, h t ∈ R n y ×1 is the nonlinear response function of the FE model at time step t, and .. u g n,1:t ∈ R (t×3)×1 is the time history of the three translational components of the FIMs at the n th support from time step 1 to time step t. The measured response vector of the bridge at the w th estimation window, y t w , is related to the FE-predicted response through a simulation error model as v t w (2) in which v t w 1 :t w 2 ∈ R (n y ×t w l )×1 is the simulation error vector at the w th estimation window and accounts for the misfit between the measured and FE-predicted responses of the bridge [23]. By neglecting the effects of modeling error, the simulation error at each time step t is ideally modeled as a spatially and temporally independent zero-mean Gaussian white noise process (i.e., v t ∼ N(0, R)). Hence, v t w is a block diagonal matrix whose block diagonals are the simulation error covariance matrix R. The objective of the model updating is to find the estimates of the unknown parameters for which the discrepancies between the measured and FE-predicted responses are minimized. The unknown parameter vector at the w th estimation window is defined as unknown parameter vector is modeled as a random vector, the evolution of which is characterized by a Gaussian Markov process-also known as a random walk. A statespace model is set up, in which the state equation governs the evolution of the random unknown parameter vector and the measurement equation corresponds to the simulation error model [24], i.e., ϕ t w in which γ w,k ∼ N(0, Q w ) is the process noise vector at the k th iteration of the w th estimation window, and Q w ∈ R n ϕ,w ×n ϕ,w is the process noise covariance matrix at the w th estimation window. Equations (3) and (4) represent a state-space model with unknown states (i.e., model parameters and FIMs herein). An Unscented Kalman Filtering (UKF) [25] method is used to estimate the unknown states. The estimation process is iterative at each estimation window. Therefore, the subscript k is added in Equations (3) and (4) to denote the iteration number.
To propagate the parameter uncertainties into the model, a scaled Unscented Transformation (UT) method is employed [26], in which the model is evaluated separately at a set of deterministically selected realizations of the parameter vector, which are referred to as the Sigma Points (SPs). The SPs at the w th estimation window are selected based on the prior mean vector (φ − t w 1 :t w 2 ) and prior covariance matrix ( P − ϕϕ t w 1 :t w 2 ) of the unknown parameters. The vector of SPs at the w th estimation window is defined as The mean vector ( y t w ) and covariance matrix of the FE-predicted responses ( P yy t w ) at the w th estimation window, as well as the cross-covariance matrix of vectors ϕ t w 1 :t w 2 and y t w , shown as P ϕy t w , are computed using a weighted sampling method as follows.
P ϕy t w In Equations (5)-(7), W j m and W j e are the mean and covariance weighting coefficients [26], respectively, andŷ t w 1 :t w 2 ϑ j w is the FE-predicted response at the w th estimation window evaluated at ϑ j w . Now, the UKF prediction-correction procedure can be employed to estimate the posterior mean vectorφ + t w 1 :t w 2 ,k+1 and posterior covariance matrix P + ϕϕ t w 1 :t w 2 ,k+1 of parameter vector at the (k + 1) th iteration. Moreover, to move from the w th estimation window to the (w + 1) th estimation window, convergence criteria, including the maximum number of iterations at each estimation window and allowable convergence tolerance in the posterior mean vector, are checked. To avoid unphysical estimates of posterior parameters, a constrain correction approach based on [27] is also implemented. The model updating algorithm is summarized in Figure 3.

Formulation for the Identifiability Analysis
Successful estimation of the unknown model parameters depends on the information that the measurement data may contain about those model parameters, as well as the sensitivity of the FE model responses with respect to those parameters. To quantify the information content of the measurement data and, therefore, to assess the identifiability of the model parameters, an approach similar to one presented in [28] is used. In this approach, the information gain of the model parameter ( ) from measurement data, ∆ ( ), is expressed as the difference between the a priori and a posteriori information entropy, which can be calculated as where is the a priori variance of and is the diagonal element of the Fisher Information Matrix defined as

Formulation for the Identifiability Analysis
Successful estimation of the unknown model parameters depends on the information that the measurement data may contain about those model parameters, as well as the sensitivity of the FE model responses with respect to those parameters. To quantify the information content of the measurement data and, therefore, to assess the identifiability of the model parameters, an approach similar to one presented in [28] is used. In this approach, the information gain of the i th model parameter θ i from measurement data, ∆H θ i , is expressed as the difference between the a priori and a posteriori information entropy, which can be calculated as where p i is the a priori variance of θ i and [I] ii is the i th diagonal element of the Fisher Information Matrix defined as in which n is the total number of time steps. The term R is the covariance matrix of the simulation error vector as defined earlier andθ is a maximum a posteriori (MAP) estimate, which is approximated with the initial estimates based on the recommendations provided in [28]. To calculate the sensitivity terms (i.e., ∂ŷ t /∂θ), a Finite Difference Method is employed. Then, the information gain of different model parameters is compared to sort out the model parameters with the highest information gain, which are relatively more likely to be identifiable unless they have strong dependence on other model parameters.
The mutual entropy gain between θ i and θ j , ∆M θ i , θ j , can be quantified through a mutual gain metric defined as where P ij is the a priori covariance matrix of θ i and θ j . In case of a strong dependency between two parameters, the model parameter with smaller information gain (based on Equation (8)) is fixed and the other model parameter will be estimated.

Verification Case Study Using the San Roque Canyon Bridge
A precast reinforced concrete bridge, referred to as the San Roque Canyon (SRC) bridge, is used as a testbed to examine the efficacy of the model updating technique. The SRC bridge, located in Santa Barbara County, CA, crosses the San Roque Creek river and is a 149-m long continuous concrete box girder bridge with a 14-m wide deck and two lanes of traffic. The SRC bridge has two concrete piers with octagonal cross-section and seat-type abutments. The bridge deck and pier cross-sections are shown in Figure 4a,b. This bridge was instrumented in 1996 with six uniaxial accelerometers on the bridge and three uniaxial accelerometers at a nearby free-field station, see Figure 4c,d. Three out of six accelerometers on the bridge measure the response of the deck in the transverse direction (channels 6, 8, and 9), and one accelerometer records the vertical motion of the deck at its center (channel 7). The channels 4 and 5 record the longitudinal response of the abutment and deck, respectively. Also, channels 1, 2, and 3 are at the free-field to collect FFMs in the horizontal (channels 1 and 3) and vertical directions (channel 2).
in which is the total number of time steps. The term R is the covariance matrix of the simulation error vector as defined earlier and is a maximum a posteriori (MAP) estimate, which is approximated with the initial estimates based on the recommendations provided in [28]. To calculate the sensitivity terms (i.e., ∂ ∂ ⁄ ), a Finite Difference Method is employed. Then, the information gain of different model parameters is compared to sort out the model parameters with the highest information gain, which are relatively more likely to be identifiable unless they have strong dependence on other model parameters. The mutual entropy gain between and , ∆ ( , ), can be quantified through a mutual gain metric defined as where is the a priori covariance matrix of and . In case of a strong dependency between two parameters, the model parameter with smaller information gain (based on Equation (8)) is fixed and the other model parameter will be estimated.

Verification Case Study Using the San Roque Canyon Bridge
A precast reinforced concrete bridge, referred to as the San Roque Canyon (SRC) bridge, is used as a testbed to examine the efficacy of the model updating technique. The SRC bridge, located in Santa Barbara County, CA, crosses the San Roque Creek river and is a 149-m long continuous concrete box girder bridge with a 14-m wide deck and two lanes of traffic. The SRC bridge has two concrete piers with octagonal cross-section and seat-type abutments. The bridge deck and pier cross-sections are shown in Figure 4a,b. This bridge was instrumented in 1996 with six uniaxial accelerometers on the bridge and three uniaxial accelerometers at a nearby free-field station, see Figure 4c,d. Three out of six accelerometers on the bridge measure the response of the deck in the transverse direction (channels 6, 8, and 9), and one accelerometer records the vertical motion of the deck at its center (channel 7). The channels 4 and 5 record the longitudinal response of the abutment and deck, respectively. Also, channels 1, 2, and 3 are at the free-field to collect FFMs in the horizontal (channels 1 and 3) and vertical directions (channel 2).  Since its instrumentation until 2021, SRC bridge has recorded seven earthquakes whose data is publicly available through the Center for Engineering Strong Motion Data (CESMD) [30]. However, only five earthquakes with Peak Ground Acceleration (PGA) greater than 0.01 g are available and all five earthquakes are considered in this study. These five earthquakes are listed in Table 1 along with their date, distance between the epicenter and the bridge, PGA, and Peak Structural Accelerations (PSAs) in different directions. The PGA for all the earthquakes, except for the 2004 IslaVista, are in horizontal direction. The recorded acceleration data for these five earthquakes are shown in Figure 5a, in which the measurement channels used later for the identifiability analysis and model updating are specified by red dashed lines. Moreover, the pseudo spectral acceleration for FFMs in transverse and longitudinal directions considering 5% damping ratio are shown in Figure 5b,c. The period of the first longitudinal and transverse modes of SRC are also shown in this figure. Since its instrumentation until 2021, SRC bridge has recorded seven earthquakes whose data is publicly available through the Center for Engineering Strong Motion Data (CESMD) [30]. However, only five earthquakes with Peak Ground Acceleration (PGA) greater than 0.01 g are available and all five earthquakes are considered in this study. These five earthquakes are listed in Table 1 along with their date, distance between the epicenter and the bridge, PGA, and Peak Structural Accelerations (PSAs) in different directions. The PGA for all the earthquakes, except for the 2004 IslaVista, are in horizontal direction. The recorded acceleration data for these five earthquakes are shown in Figure  5a, in which the measurement channels used later for the identifiability analysis and model updating are specified by red dashed lines. Moreover, the pseudo spectral acceleration for FFMs in transverse and longitudinal directions considering 5% damping ratio are shown in Figure 5b

FE Model of the SRC Bridge
A detailed FE model of the SRC bridge is developed in OpenSees [31] following the guidelines provided in [32]. All the nominal material properties are taken from the as-built structural drawings and Caltrans Seismic Design Criteria [33]. Figure 6 is a schematic representation of the FE model, in which the model components are numbered from 1 to 13 and explained in the following text. Moreover, the model parameters are numbered from 1 to 34 and their nominal values are listed in Table 2. The far-field soil-embankment effect is modeled through three parallel zeroLength spring and dashpot elements (Component #12), with the stiffness and damping properties calculated based on [42]. The springs have elastic-no-tension (ENT) uniaxial material with a modulus of elasticity equal to 3 GPa. The dashpots have Viscous uniaxial material with a damping coefficient defined as the summation of linear radiation and material damping , coefficients. The parameters and , are taken as 38 and 140 MN. sec m ⁄ , respectively.
At the foundations, the near-field SSI effects are neglected and the far-field effects under the piers and abutments foundations (Component #13) in all six degrees of freedom are modeled using zeroLength elements (springs and dashpots) with stiffness and damping coefficients calculated based on [43] In addition to the explicit sources of damping energy dissipation in the model (dashpots and nonlinear materials), the inherent damping property of the structure is modeled using Rayleigh damping with the mass ( ) and stiffness ( ) proportional coefficients equal to 0.6 and 0.003, respectively, which results in 5% damping for modes 1 and 6. The Rayleigh damping is modeled using the initial stiffness matrix.  Table 2.  Table 2.  The deck (Component #1) is modeled using elasticBeamColumn elements, as it is a capacity protected element [33]. Each span is meshed using ten elements; the cross-sectional properties are calculated based on the section geometry shown in Figure 4a and the material modulus of elasticity (E d ) is 27.8 GPa. To account for the torsional vibrations, the rotational mass moment of inertia is calculated based on the section geometry and is added to deck nodes (Component #2). The pier foundations are modeled using four 4 m × 4 m linearelastic shellMITC4 elements (Component #3) with 1.7 m thick ElasticMembranePlateSection section and modulus of elasticity of 27.8 GPa. To connect columns to the deck, rigidLink elements (Component #4) are used to account for the rigidity provided by the cap beams. Fiber-section forceBeamColumn elements with five integration points [34] are utilized to model the piers (Component #5). The elastic shear and torsional stiffnesses are aggregated to the pier fiber sections. The confinement effects in the pier sections are considered using the Mander's model [35]. The confined core is modeled using Concrete04 with compressive strength f c,c of 40.4 MPa achieved at 0.37% strain and initial modulus of elasticity (E c ) of 27.8 GPa. Moreover, reinforcement is modeled using Steel02 material with a yield strength of 46.9 GPa and a modulus of elasticity of 200 GPa. To connect the deck to the abutment system, rigid frame elements (Component #6) are added to the two ends of the deck with the number of nodes equal to the number of bearing pads. These nodes are connected to the rigid abutment body through bearing pads and shear keys. Bearing pads in the vertical direction are modeled using zeroLength elements with ElasticPPGap uniaxial material with compression-only vertical compliance up to a 0.01 m deformation, known as the initial gap (Component #7). The stiffness in bearing pads (Component #8) is modeled using zeroLength elements with uniaxial bilinear steel01 material assuming 10 8 N m and 10 7 N m shear stiffness in the longitudinal k b L and transverse k b T directions, respectively, and 1% strain-hardening ratio. Once the gap between the deck and the abutment backwall is closed, the backwall provides resistant in the longitudinal direction. The backwall stiffness (Component #9), which is calculated using the theory of plate and shell [36], is modeled using compression-only elastic perfectly-plastic gap behavior. For this purpose, an ElasticPPGap uniaxial material with a 0.05 m gap is used. The transverse response of the deck is controlled by the zeroLength elements representing the shear keys (Component #10) with a shear stiffness of 100 MN/m. According to the structural drawings, the shear keys of SRC bridge are ductile. So, the model proposed in [37] can be used to model the nonlinear compression-only behavior of the shear keys after the gap between the deck and shear keys is closed. However, since the level of excitation in this study is low, the shear keys are not modeled and it is assumed that their stiffness is lumped to the stiffness of bearing pads.

Identifiability Analysis
SSI effects are modeled to consider near-field and far-field effects [19,20,38]. At the abutments, the near-filed effects include the passive soil pressure behind the backwall, which provides resistance against the abutment movement. This passive soil pressure is modeled using zeroLength soil springs (Component #11) with a HyperbolicGapMaterial based on the Generalized Hyperbolic Force-Displacement (GHFD) backbone curve [39]. The backbone curve is developed using the abutment wall height (2.44 m), wall-soil interface friction angle (40 • ), and Mohr-Coulomb shear strength parameters with a common silty-sand soil type [40], whose properties are taken from [41]. In addition to the explicit sources of damping energy dissipation in the model (dashpots and nonlinear materials), the inherent damping property of the structure is modeled using Rayleigh damping with the mass (a 0 ) and stiffness (a 1 ) proportional coefficients equal to 0.6 and 0.003, respectively, which results in 5% damping for modes 1 and 6. The Rayleigh damping is modeled using the initial stiffness matrix.

Identifiability Analysis
An identifiability analysis is performed to find the most identifiable model parameters based on the approach described in Section 2.2. The identifiability analysis is an inputoutput approach measuring the sensitivity of model responses with respect to different model parameters. The analysis is completed here using the significant portion of the 2004 Isla Vista earthquake, and the FFMs are used in lieu of FIMs. Based on the instrumentation plan in Figure 4, channel 7 records the vertical component of the bridge response and is likely polluted by the traffic-induced vibrations and thus not included for the identifiability analysis. Consequently, the vertical component of the input motion (channel 2), is also removed from the analysis. The identifiability analysis process is discussed in the following part.
A list of potential model parameters and their nominal values are shown in Table 2. The longitudinal, transverse, and vertical directions mentioned in this table are based on the directions shown in Figure 6. Model parameters are numbered from 1 to 34. These parameter numbers are used in Figure 6 to relate the model parameters to the corresponding model components. To find the most identifiable model parameters, the relative information gain for each model parameter is calculated based on Equation (8) and the results are shown in Figure 7. As can be seen in this figure, there are 10 model parameters that are potentially identifiable due to their high relative information gain in comparison to other model parameters. These model parameters, sorted in descending order with respect to their relative information gain, are: a 0 , a 1 , c a L , k p R,L , and k p R,T . As can be seen in Figure 7, the model parameters related to the nonlinear response of bridge, including f c,c , are likely unidentifiable. This is because the considered earthquake in this study has small amplitude. Figure 8a shows the relative mutual information gain among the model parameter pairs obtained from Equation (10). Moreover, to better observe the dependency between model parameters, a scaled version of Figure 8a is presented in Figure 8b. In this figure, all the mutual information gains in each row are scaled to the maximum value in the corresponding row. Also, the diagonal values are nullified (actual values are replaced by zero). As can be seen, there is no significant dependence between the first four model parameters with high relative information gain (k b L , E d , k b T , and E c ). So, these parameters are chosen as unknown model parameters to be estimated through the model updating process. Parameters k p R,T and k p R,L are removed from the unknowns because they are dependent on parameters k b L and E c , respectively. Furthermore, parameters a 0 and a 1 are moderately dependent; however, both are kept among the unknowns to be estimated. Parameters k a L and c a L are dependent on E d and k b L . This is likely due to the competing effects that parameters k a L , E d , and k b L have on the stiffness of the superstructure in the longitudinal direction. Since E d is already selected as an unknown model parameter to be estimated, k a L and c a L are excluded from the unknown parameter vector for the verification studies. In summary, parameters E d , E c , k b T , k b L , a 0 , and a 1 , highlighted in Table 2, are selected as unknown model parameters to be estimated.     Table 2.

Verification Study Using Numerically Simulated Data
To simulate the measurements, the FE response of the SRC bridge to the 2004 Isla Vista earthquake is simulated and polluted with an artificial zero-mean Gaussian noise with 5% Root Mean Square (RMS) noise-to-signal ratio. For the simulation, the measured free field motions are used as the FIMs. The time history analysis is performed using time step size of 0.01 sec and the Newmark-beta average acceleration method for time integration, with the nominal model parameter values presented in Table 2. The noisy simulated responses are used for model updating to estimate the six unknown model parameters and the FIM time histories in the longitudinal and transverse directions. It should be noted that due to the short length of the SRC bridge, the FIMs are assumed to be uniform. As

Verification Study Using Numerically Simulated Data
To simulate the measurements, the FE response of the SRC bridge to the 2004 Isla Vista earthquake is simulated and polluted with an artificial zero-mean Gaussian noise with 5% Root Mean Square (RMS) noise-to-signal ratio. For the simulation, the measured free field motions are used as the FIMs. The time history analysis is performed using time step size of 0.01 sec and the Newmark-beta average acceleration method for time integration, with the nominal model parameter values presented in Table 2. The noisy simulated responses are used for model updating to estimate the six unknown model parameters and the FIM time histories in the longitudinal and transverse directions. It should be noted that due to the short length of the SRC bridge, the FIMs are assumed to be uniform. As discussed before, measurement channels 7 and 2 are not considered and only channels 4, 5, 6, 8, and 9 are used in the model updating.
The model updating is carried out using 31 number of overlapping estimation windows (n w = 31). The length of each estimation window is 27 time steps (t w l = 27), and the overlap between each two estimation windows is equal to 10 time steps (t w,w+1 o = 10). In order to initialize the model parameter vector (θ 0 ), −20% estimation error is assumed in the initial estimates of all unknown model parameters where estimation error is defined as Estimation error(%) = Estimated model parameter value Nominal model parameter value − 1 × 100 (11) and the nominal values are shown in Table 2. The assumed −20% deviation of model parameters from their nominal values is chosen as a reasonable initial error based on engineering judgement. The final estimates of the model parameters are sensitive to the initial estimates. However, as long as the initial estimates are not unreasonably far from the nominal values, no significant change in the final estimates is expected [44]. The term P θθ t 0 1 :t 0 2 is set as a diagonal matrix with i th diagonal entry equal to 0.1θ 0,i 2 , whileθ 0,i is the i th entry ofθ 0 . Moreover, P.. is also initialized as a t 0 l × t 0 l diagonal matrix with diagonal entries equal to (0.001g) 2 . The first n θ diagonal entries in Q w are equal to 0.001θ w,i 2 , whereθ w,i is the i th entry of theθ w , and the rest are time-invariant and equal to 10 −8 -note that the termθ w refers to the estimate of model parameter vector at the w th estimation window. Finally, the diagonal entries of the simulation error covariance matrix, R, representing the measurement noise variances, are set as (0.003%g) 2 . Moreover, the maximum number of iterations and the allowable convergence tolerance are selected equal to 9 and 0.02. These thresholds are based on experience to achieve a compromise between estimation accuracy and computational demand. The time histories of the estimation error for unknown model parameters are shown in Figure 9. In this figure, the error in the estimation of model parameters are calculated at the end of each estimation window. As can be seen, the estimation errors for most unknown model parameters are significant in the first five seconds. However, as time passes and more information is collected, the estimation errors decrease and converge to zero with less than 4% error. Also, as seen in this figure, elastic modulus of deck and initial elastic modulus of columns converge to their nominal values faster than the other unknown model parameters. This is because they are associated with the initial linear-elastic stiffness of the bridge structural components.
ant and equal to 10 -note that the term refers to the estimate of mode vector at the estimation window. Finally, the diagonal entries of the simu covariance matrix, , representing the measurement noise variances, (0.003% ) . Moreover, the maximum number of iterations and the allowa gence tolerance are selected equal to 9 and 0.02. These thresholds are based on to achieve a compromise between estimation accuracy and computational de The time histories of the estimation error for unknown model parameter in Figure 9. In this figure, the error in the estimation of model parameters ar at the end of each estimation window. As can be seen, the estimation errors f known model parameters are significant in the first five seconds. However, as and more information is collected, the estimation errors decrease and conv with less than 4% error. Also, as seen in this figure, elastic modulus of deck elastic modulus of columns converge to their nominal values faster than th known model parameters. This is because they are associated with the initial l stiffness of the bridge structural components. The estimated FIMs in the longitudinal and transverse directions are c their true counterparts in Figure 10a,b, respectively. Moreover, to evaluate t The estimated FIMs in the longitudinal and transverse directions are compared to their true counterparts in Figure 10a,b, respectively. Moreover, to evaluate the accuracy of the estimated FIMs, the Relative Root Mean Square Error (RRMSE) of the estimated FIM time histories are shown in Figure 10c. The RRMSE at time t n is calculated as where the termsŝ i and s i are the estimated and true time history values at the i th time step, respectively. As seen in Figure 10, the estimated FIMs are less accurate at the beginning of estimation, when model parameters are not accurately estimated yet, and as time passes the RRMSEs decrease.
RRMSE. (%) = ∑ ( ) × 100 (12) where the terms and are the estimated and true time history values at the time step, respectively. As seen in Figure 10, the estimated FIMs are less accurate at the beginning of estimation, when model parameters are not accurately estimated yet, and as time passes the RRMSEs decrease.  Figure 11 shows a comparison between measured and estimated responses at the measurement channels. As can be seen, the match between the estimated and measured response time histories is reasonable, which verifies the model updating technique. The RRMSE values for the estimated responses decrease in time due to improved accuracy in the estimated FIMs and unknown model parameters.

Case Study Using Real-World Earthquake Data
In this section, the SRC bridge model is integrated with the real data collected during earthquakes shown in Table 1. Using the presented model updating technique, the 6 unknown model parameters, already specified in Section 3.2, are estimated jointly with the two FIMs in the longitudinal and transverse directions. Details of the model updating process is similar to those explained in the previous section, and the nominal parameter  × 100 (12) where the terms ̂ and are the estimated and true time history values at the ℎ time step, respectively. As seen in Figure 10, the estimated FIMs are less accurate at the beginning of estimation, when model parameters are not accurately estimated yet, and as time passes the RRMSEs decrease.  Figure 11 shows a comparison between measured and estimated responses at the measurement channels. As can be seen, the match between the estimated and measured response time histories is reasonable, which verifies the model updating technique. The RRMSE values for the estimated responses decrease in time due to improved accuracy in the estimated FIMs and unknown model parameters.

Case Study Using Real-World Earthquake Data
In this section, the SRC bridge model is integrated with the real data collected during earthquakes shown in Table 1. Using the presented model updating technique, the 6 unknown model parameters, already specified in Section 3.2, are estimated jointly with the two FIMs in the longitudinal and transverse directions. Details of the model updating process is similar to those explained in the previous section, and the nominal parameter

Case Study Using Real-World Earthquake Data
In this section, the SRC bridge model is integrated with the real data collected during earthquakes shown in Table 1. Using the presented model updating technique, the 6 unknown model parameters, already specified in Section 3.2, are estimated jointly with the two FIMs in the longitudinal and transverse directions. Details of the model updating process is similar to those explained in the previous section, and the nominal parameter values (shown in Table 2) are used as the initial estimates for the unknown model parameters. Figure 12 shows the estimation time history of the unknown model parameters. In this figure, the parameters' estimates are normalized to their nominal (initial) values, and the term E ave is the average of estimates for the concrete modulus of elasticity in deck (E d ) and the initial modulus of elasticity in columns (E c ). The event-to-event variabilities in the final estimates of the unknown model parameters are shown in Figure 13. In this figure, the nominal values are shown as dashed red lines. As seen in Figure 13a, while the final estimates for the concrete modulus of elasticity are close to the nominal values, there is minor event-to-event variability. These minor variations could be due to the concrete aging, modeling uncertainties, and other sources of modeling error. The variation in the final estimates of elastomeric bearing pad's shear stiffness in longitudinal and transverse directions are shown in Figures 13b and 13c, respectively. As can be seen, the level of the final estimates in the transverse direction is higher than the ones in the longitudinal direction. This is likely because the model parameter k b T represents the combined effects of the stiffness of the elastomeric pad and shear key, while the model parameter k b L only accounts for stiffness of the elastomeric pad. The level of the final estimates in Figure 13b,c are close to each other for different earthquakes with the exception of the 2018 Santa Cruz earthquake (earthquake No. 5). This variation can be due to the fact that the 2018 Santa Cruz earthquake is one of the weakest events used in this study with small PGA and PSAs (see Table 1). In such a weak event, the bearing pad stiffness is likely high enough to prevent any movement. In this case, the measurements may not be sensitive to the stiffness of the bearing pads, which could result in large estimates for k b L and k b T . To better investigate the estimated Rayleigh damping, the final estimates of the damping coefficients (a 0 and a 1 ) are used to calculate the damping ratio as a function of frequency as shown in Figure 14. As can be seen, the model updating results show almost no damping during the first earthquake, i.e., 2003 San Simeon. This is because that the bridge behaves almost quasi-statically in this excitation due to the superior presence of low-frequency components in 2003 San Simeon in comparison to the other earthquakes (see Figure 5). On the other hand, the highest damping ratios are estimated from the 2013 Isla Vista earthquake with the highest PGA, which suggests a correlation between the inherent damping and earthquake intensity. Based on these results, considering a 3-4% Rayleigh damping ratio for the first and eighth modes (with frequencies of 1.12 and 5.82 Hz, respectively) can be a rational assumption for the SRC bridge in the case of low to moderate earthquakes, i.e., 2013 Isla Vista. However, under weak excitations, e.g., the 2004 Isla Vista earthquake, lower damping ratio, around 1-2%, can be considered for these modes.  Table 2) are used as the initial estimates for the unknown model parameters. Figure 12 shows the estimation time history of the unknown model parameters. In this figure, the parameters' estimates are normalized to their nominal (initial) values, and the term is the average of estimates for the concrete modulus of elasticity in deck ( and the initial modulus of elasticity in columns . The event-to-event variabilities in the final estimates of the unknown model parameters are shown in Figure 13. In this figure, the nominal values are shown as dashed red lines. As seen in Figure 13a, while the final estimates for the concrete modulus of elasticity are close to the nominal values, there is minor event-to-event variability. These minor variations could be due to the concrete aging, modeling uncertainties, and other sources of modeling error. The variation in the final estimates of elastomeric bearing pad's shear stiffness in longitudinal and transverse directions are shown in Figure 13b and Figure 13c, respectively. As can be seen, the level of the final estimates in the transverse direction is higher than the ones in the longitudinal direction. This is likely because the model parameter represents the combined effects of the stiffness of the elastomeric pad and shear key, while the model parameter only accounts for stiffness of the elastomeric pad. The level of the final estimates in Figure 13b,c are close to each other for different earthquakes with the exception of the 2018 Santa Cruz earthquake (earthquake No. 5). This variation can be due to the fact that the 2018 Santa Cruz earthquake is one of the weakest events used in this study with small PGA and PSAs (see Table 1). In such a weak event, the bearing pad stiffness is likely high enough to prevent any movement. In this case, the measurements may not be sensitive to the stiffness of the bearing pads, which could result in large estimates for and .   To better investigate the estimated Rayleigh damping, the final estimates of the damping coefficients ( and ) are used to calculate the damping ratio as a function of frequency as shown in Figure 14. As can be seen, the model updating results show almost no damping during the first earthquake, i.e., 2003 San Simeon. This is because that the bridge behaves almost quasi-statically in this excitation due to the superior presence of low-frequency components in 2003 San Simeon in comparison to the other earthquakes (see Figure 5). On the other hand, the highest damping ratios are estimated from the 2013 Isla Vista earthquake with the highest PGA, which suggests a correlation between the inherent damping and earthquake intensity. Based on these results, considering a 3-4% Rayleigh damping ratio for the first and eighth modes (with frequencies of 1.12 and 5.82 Hz, respectively) can be a rational assumption for the SRC bridge in the case of low to moderate earthquakes, i.e., 2013 Isla Vista. However, under weak excitations, e.g., the 2004 Isla Vista earthquake, lower damping ratio, around 1-2%, can be considered for these modes.   To better investigate the estimated Rayleigh damping, the final estimat damping coefficients ( and ) are used to calculate the damping ratio as a fu frequency as shown in Figure 14. As can be seen, the model updating results sho no damping during the first earthquake, i.e., 2003 San Simeon. This is because bridge behaves almost quasi-statically in this excitation due to the superior pr low-frequency components in 2003 San Simeon in comparison to the other ear (see Figure 5). On the other hand, the highest damping ratios are estimated from Isla Vista earthquake with the highest PGA, which suggests a correlation betwee herent damping and earthquake intensity. Based on these results, considering a 3 leigh damping ratio for the first and eighth modes (with frequencies of 1.12 and respectively) can be a rational assumption for the SRC bridge in the case of low t ate earthquakes, i.e., 2013 Isla Vista. However, under weak excitations, e.g., the Vista earthquake, lower damping ratio, around 1-2%, can be considered for thes The comparisons between the FFMs and the estimated FIMs and the corresponding RRMSEs in the longitudinal and transverse directions are shown in Figures 15 and 16, respectively. Here, RRMSE presents the misfit between the predicted FIMs and collected FFMs. So, greater values for RRMSEs can imply more significant SSI effects. Note that as channel 3 malfunctioned during the 2017 Montecito and 2018 Santa Cruz earthquakes, measurements of channels 4 and 1 are considered to report the recorded FFMs in the longitudinal and transverse directions, respectively, for these two earthquakes. In Figure 15, the estimated FIMs in the longitudinal direction are compared with the recorded measurements in channel 4, and as can be seen, these two signals match well for these two earthquakes. In other words, what is recorded by measurement channel 4 on the abutment overlays the estimated FIM. This means that there is likely no significant bridge-abutment interaction. The best fit between the FIMs and FFMs is seen for the 2003 San Simeon earthquake, which is not surprising due to the low-frequency content of this motion as discussed earlier.
Overall, the estimated longitudinal FIMs in all earthquakes show close similarity with the recorded FFMs. Therefore, neither kinematic interaction nor inertial interaction is likely significant for this bridge under the studied ground motion intensities in the longitudinal direction; however, a higher level of differences is observed between recorded FFMs and estimated FIMs in the transverse direction.
quakes, measurements of channels 4 and 1 are considered to report the recorded FFMs in the longitudinal and transverse directions, respectively, for these two earthquakes. In Figure 15, the estimated FIMs in the longitudinal direction are compared with the recorded measurements in channel 4, and as can be seen, these two signals match well for these two earthquakes. In other words, what is recorded by measurement channel 4 on the abutment overlays the estimated FIM. This means that there is likely no significant bridge-abutment interaction. The best fit between the FIMs and FFMs is seen for the 2003 San Simeon earthquake, which is not surprising due to the low-frequency content of this motion as discussed earlier. Overall, the estimated longitudinal FIMs in all earthquakes show close similarity with the recorded FFMs. Therefore, neither kinematic interaction nor inertial interaction is likely significant for this bridge under the studied ground motion intensities in the longitudinal direction; however, a higher level of differences is observed between recorded FFMs and estimated FIMs in the transverse direction.   Figures 17-21 show the comparison between the measured and estimated (from the updated models) acceleration responses. As seen, the predicted responses match well the measured responses in all events. However, as previously discussed in Section 3.3 and can be seen in the following figures, the accuracy in estimation of the responses is poor at the beginning of the earthquake event and improves in time.   show the comparison between the measured and estimated (from the updated models) acceleration responses. As seen, the predicted responses match well the measured responses in all events. However, as previously discussed in Section 3.3 and can be seen in the following figures, the accuracy in estimation of the responses is poor at the beginning of the earthquake event and improves in time.  show the comparison between the measured and estimated (from the updated models) acceleration responses. As seen, the predicted responses match well the measured responses in all events. However, as previously discussed in Section 3.3 and can be seen in the following figures, the accuracy in estimation of the responses is poor at the beginning of the earthquake event and improves in time.

Digital Twin and Virtual Sensing Application
The updated FE model is a replicate of the SRC bridge in the digital world given the training data, and is referred to as its digital twin (DT). The DT can be used to reconstruct local element-and material-level responses to diagnose damage at local levels of the bridge. For this purpose, the DT is used in a forward simulation for a given input motion and the local responses (e.g., section moment curvature, material stress-strain) are collected and used for damage localization and quantification. The process is often referred to as virtual sensing [45]. Although the level of available earthquake motions was small in this study and thus, no damage is expected in the structure, the application of the DT can still be demonstrated. It should be noted that the phrase "Digital Twin" may have different meanings and applications in different industries. Here, the digital twinning method is used for post-earthquake assessment of bridges. The data recorded in seismic events are considered as discontinuous real-time data assuming that no significant changes will happen between the events. Therefore, the DT will be updated after each seismic event and used for post-earthquake damage diagnosis. In this section, the SRC DT is established using the final estimates of the unknown model parameters from the 2013 IslaVista earthquake data (earthquake No. 3). Then, the estimated FIMs are used with the DT to monitor/estimate the structural responses. Figure 22 shows two examples of such monitored/estimated responses at local levels. Figure 22a presents the estimated momentcurvature response at the lowest section of the east pier about the transverse direction.

Digital Twin and Virtual Sensing Application
The updated FE model is a replicate of the SRC bridge in the digital world given the training data, and is referred to as its digital twin (DT). The DT can be used to reconstruct local element-and material-level responses to diagnose damage at local levels of the bridge. For this purpose, the DT is used in a forward simulation for a given input motion and the local responses (e.g., section moment curvature, material stress-strain) are collected and used for damage localization and quantification. The process is often referred to as virtual sensing [45]. Although the level of available earthquake motions was small in this study and thus, no damage is expected in the structure, the application of the DT can still be demonstrated. It should be noted that the phrase "Digital Twin" may have different meanings and applications in different industries. Here, the digital twinning method is used for post-earthquake assessment of bridges. The data recorded in seismic events are considered as discontinuous real-time data assuming that no significant changes will happen between the events. Therefore, the DT will be updated after each seismic event and used for post-earthquake damage diagnosis. In this section, the SRC DT is established using the final estimates of the unknown model parameters from the 2013 IslaVista earthquake data (earthquake No. 3). Then, the estimated FIMs are used with the DT to monitor/estimate the structural responses. Figure 22 shows two examples of such monitored/estimated responses at local levels. Figure 22a presents the estimated moment-curvature response at the lowest section of the east pier about the transverse direction. Also, the predicted stress-strain response of the extreme concrete fiber at the lowest section of the west pier is shown in Figure 22b. As can be seen, the structure behaves in its linear-elastic regime during the 2013 IslaVista earthquake. However, a similar concept is applicable for stronger earthquakes to monitor the response nonlinearity at the local levels to infer potential damage location and extent.
Also, the predicted stress-strain response of the extreme concrete fiber at the lowest section of the west pier is shown in Figure 22b. As can be seen, the structure behaves in its linear-elastic regime during the 2013 IslaVista earthquake. However, a similar concept is applicable for stronger earthquakes to monitor the response nonlinearity at the local levels to infer potential damage location and extent.

Conclusions
This study presented the procedure for bridge digital twinning and virtual sensing through the application of an output-only time-domain model updating technique for post-earthquake damage diagnosis. In this technique the output-only seismic responses are integrated with the mechanics-based Finite Element (FE) model of a bridge through a Bayesian inference approach to jointly estimate the Foundation Input Motions (FIMs) and the unknown model parameters, and so the development of the bridge's digital twin. The model updating is implemented using a sequential estimation window approach, in which the time domain is divided into sequential overlapping estimation windows. At each estimation window, the Probability Density Functions (PDF) of input and parameters are updated iteratively using an Unscented Kalman filtering method and transferred to the next estimation window. The technology solution for post-earthquake damage diagnosis will necessitate seamless implementation using parallel processing and high-performance computing (HPC) to enable a near real-time (e.g., within hours) processing capability.
The output-only time-domain model updating technique was first verified in a numerically simulated environment. For this purpose, a detailed FE model of the San Roque Canyon (SRC) bridge, located in CA, was developed in OpenSees. An identifiability analysis was performed to select the likely identifiable model parameters among 34 parameter candidates. Due to instrumentation sparsity and the low-amplitude seismic inputs, only a few (six) linear-related model parameters (concrete modulus of elasticity of deck, initial concrete modulus of elasticity of piers, stiffness of bearing pads, and Rayleigh damping coefficients) were found to be likely identifiable. The verification study presented good performance of the model updating technique in estimating the unknown model parameters jointly with the time history of FIMs.
In the next step, the model updating technique was repeated for five real-world earthquake data recorded at the SRC bridge station from 2003 to 2018. The estimated model parameters were compared between different earthquakes and the likely reasons for discrepancies were discussed. Estimates of the Rayleigh damping coefficients for different earthquakes revealed dependency of damping ratio on the intensity of the seismic excitation, i.e., a stronger earthquake with greater PGA resulted in a larger estimate of damping ratio. The estimates of bearing pad stiffness had large event-to-event variations, Figure 22. Local level response monitored in SRC bridge using its digital twin during 2013 IslaVista earthquake event: (a) moment-curvature response (about the transverse direction) at the lowest section of the east pier, and (b) estimated stress-strain response in the extreme fiber at the lowest section of the west pier.

Conclusions
This study presented the procedure for bridge digital twinning and virtual sensing through the application of an output-only time-domain model updating technique for post-earthquake damage diagnosis. In this technique the output-only seismic responses are integrated with the mechanics-based Finite Element (FE) model of a bridge through a Bayesian inference approach to jointly estimate the Foundation Input Motions (FIMs) and the unknown model parameters, and so the development of the bridge's digital twin. The model updating is implemented using a sequential estimation window approach, in which the time domain is divided into sequential overlapping estimation windows. At each estimation window, the Probability Density Functions (PDF) of input and parameters are updated iteratively using an Unscented Kalman filtering method and transferred to the next estimation window. The technology solution for post-earthquake damage diagnosis will necessitate seamless implementation using parallel processing and high-performance computing (HPC) to enable a near real-time (e.g., within hours) processing capability.
The output-only time-domain model updating technique was first verified in a numerically simulated environment. For this purpose, a detailed FE model of the San Roque Canyon (SRC) bridge, located in CA, was developed in OpenSees. An identifiability analysis was performed to select the likely identifiable model parameters among 34 parameter candidates. Due to instrumentation sparsity and the low-amplitude seismic inputs, only a few (six) linear-related model parameters (concrete modulus of elasticity of deck, initial concrete modulus of elasticity of piers, stiffness of bearing pads, and Rayleigh damping coefficients) were found to be likely identifiable. The verification study presented good performance of the model updating technique in estimating the unknown model parameters jointly with the time history of FIMs.
In the next step, the model updating technique was repeated for five real-world earthquake data recorded at the SRC bridge station from 2003 to 2018. The estimated model parameters were compared between different earthquakes and the likely reasons for discrepancies were discussed. Estimates of the Rayleigh damping coefficients for different earthquakes revealed dependency of damping ratio on the intensity of the seismic excitation, i.e., a stronger earthquake with greater PGA resulted in a larger estimate of damping ratio. The estimates of bearing pad stiffness had large event-to-event variations, which were related to the level of intensity and potential modeling errors. Bearing pads were modeled as linear springs, while their stiffness often depends on the level and frequency content of the excitation. Furthermore, the estimated FIMs were compared with the Free-Field Motions (FFMs), and the discrepancies were discussed based on the expected soil-structure interaction effects. As expected, smaller discrepancies were observed between the FFMs and the estimated FIMs for the events with low frequency excitation. Furthermore, the measured responses of the bridge were compared with those estimated from the updated model (or digital twin). While the agreement between the measured and digital twin predicted responses was acceptable, a comparison of the relative root mean square errors between the numerically simulated case study and the real-world case studies clearly showed the inherent effects of modeling error in the FE model updating technique. Modeling errors refer to the mathematical idealization and simplifications in the model, which can result in biased and incorrect model updating outcomes. Finally, the digital twin was developed using the strongest studied earthquake. To demonstrate the virtual sensing application, the bridge's digital twin was used to predict the local level response of the bridge given the strongest studied earthquake, an important capability that can be used for post-event damage diagnosis (damage detection, localization, and quantification). Although the levels of considered earthquakes were low in this study, the approach can be applied to large-scale and nonlinear models to diagnose damage as the result of material nonlinearity.
The objective of this study was to examine the application of the time domain outputonly Bayesian FE model updating for damage identification and bridge digital twinning through a real-world case study. Further studies of this type with various structural systems (simple and complex, experimental, and real-world) and under various loading conditions are needed to better understand the limitation and capabilities of the Bayesian model updating approach as a technology solution for structural health monitoring and damage diagnosis.