A Two-Step FE Model Updating Approach for System and Damage Identiﬁcation of Prestressed Bridge Girders

: This study presents a two-step FE model updating approach for health monitoring and damage identiﬁcation of prestressed concrete girder bridges. To reduce the effects of modeling error in the model updating process, in the ﬁrst step, modal-based model updating is used to estimate linear model parameters mainly related to the stiffness of boundary conditions and material properties. In the second step, a time-domain model updating is carried out using acceleration data to reﬁne parameters accounting for the nonlinear response behavior of the bridge. In this step, boundary conditions are ﬁxed at their ﬁnal estimates using modal-based model updating. To prevent the convergence of updating algorithm to local solutions, the initial estimates for nonlinear material properties are selected based on the ﬁrst-step model updating results. To validate the applicability of the two-step FE model updating approach, a series of forced-vibration experiments are designed and carried out on a pair of full-scale decommissioned and deteriorated prestressed bridge I-girders. In the ﬁrst step, parameters related to boundary conditions, including stiffness of supports and coupling beams, as well as material properties, including initial stiffness of concrete material, are estimated. In the second step, concrete compressive strength and damping properties are updated. The ﬁnal estimates of the concrete compressive strength are used to infer the extent of damage in the girders. The obtained results agree with the literature regarding the extent of reduction in concrete compressive strength in deteriorated concrete structures.


Introduction
Bridges are vital components of the transportation infrastructure. The average age of in-service bridges in the United States is increasing, which necessitates methods and tools to inform decision making related to the maintenance and/or replacement of these structures [1,2]. Finite Element (FE) model updating methods have emerged as a venerable procedure for operational health monitoring and post-event structural damage identification [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. In these methods, the initial/baseline FE model-developed using available as-built drawings-is updated using measured dynamic responses. During this process, uncertain model parameters-including material properties, damping parameters, boundary conditions, etc.-are calibrated/estimated. The deviation of final estimates of model parameters from their initial/baseline values reveals information regarding the location and extent of damage in the structure.
FE model updating approaches are mainly divided into two groups. The first group is modal-based model updating, wherein the initial FE model is updated to match the identified modal properties of the structure. In this method, the modal properties are first identified using modal identification methods (e.g., [18][19][20][21][22]). Then, the parameters characterizing the linear response behavior of the FE model are estimated to reduce the discrepancies between the identified and FE-predicted modal properties [3]. The accuracy of the identified modal signatures is controlled by the level of nonlinearity in the response behavior of structure, measurement noise, excitation frequency range, and sensor sparsity. The uncertainty in the identified modal properties propagates into the model updating process and is reflected in the final estimates of model parameters [23]. In addition, modal properties cannot be used to infer parameters characterizing the nonlinear material behavior. Consequently, the updated FE model may not be able to predict the dynamic response behavior of the structure correctly, especially when the structure is subjected to material nonlinearity. Although the estimated linear model parameters can be used for structural monitoring and damage identification [24], the application of modal-based model updating for damage identification of reinforced concrete bridges has been shown to be limited [25,26]. The second group of model updating approaches is referred to as time-domain model updating [17,27,28]. In this approach, the unknown model parameters characterizing linear and/or nonlinear response behavior of structure are updated to reduce the discrepancies between the measured and FE-predicted responses in time domain. In contrast to the modal-based model updating method, the measured dynamic responses of the structure are used directly for model updating in this approach. The direct application of dynamic responses in time-domain model updating eliminates the propagation of modal identification uncertainties into the model updating process.
Several studies in the literature (e.g., [10,23,[29][30][31][32][33][34][35]) have focused on the system and damage identification of bridges subjected to ambient or traffic excitation (i.e., operational conditions) using modal-based model updating. The performance of these methods is mostly evaluated by comparing the identified and posterior FE-predicted modal signatures of the bridge. Studies in [36,37] showed that the accuracy of the updated model is highly sensitive to the selection of unknown model parameters. Moreover, [38] indicates that damage detection of bridges would depend on the proper simulation of boundary conditions. A two-step FE model updating process is suggested in this study to resolve modeling errors due to boundary conditions. Using measurements other than modal properties for model updating and damage identification of bridges has attracted research interest recently. While static/pseudo-static responses (e.g., displacement measurements) have been used in previous studies for model updating [5,11,17,[39][40][41][42], acceleration measurements have not been used directly for the purpose of model updating of bridges under operational conditions. It is worth noting that static/pseudo-static responses contain limited information regarding the dynamic behavior of the bridge compared to acceleration measurements. Therefore, acceleration measurements can be more informative about the uncertain model parameters compared to static/pseudo-static responses. The target in this study is to use both the modal properties and acceleration responses for model updating and damage diagnosis through a two-step FE model updating process.
Moreover, previous studies have shown that weak identifiability and mutual dependency between model parameters, modeling errors, as well as convergence of parameters to local solution may challenge the model updating process [4,43,44]. These challenges are exacerbated in a real-world application, especially in cases with large number of uncertain model parameters and/or improper selection of initial values for the uncertain model parameters.
To resolve the above-mentioned issues, this study presents a sequential combination of modal-based and time-domain model updating for operational health monitoring and damage identification of aged bridges using acceleration responses. In this procedure, first, a deterministic modal-based model updating is carried out to estimate the linear model parameters of a bridge. These model parameters are related to boundary conditions and material properties. Then, in order to refine the parameter estimation and account for the nonlinear response behavior of the bridge, a time-domain model updating is carried out. In this step, nonlinear material properties, as well as the damping energy-dissipation-related In summary, the reasoning behind introducing the sequential combination of modalbased and time-domain model updating, which constitute the novelties of this work, can be summarized as follows.

•
Unknown boundary conditions often challenge the application of time-domain model updating for bridges since the model parameters are often dependent on the boundary conditions. Here, the modal-based model updating is used to identify the boundary conditions first.

•
The application of modal-based model updating for damage identification of bridges is limited. This is likely due to the uncertainties in the identified modal signatures that propagate through the parameter estimation process. Therefore, here, the estimation of model parameters for damage identification will be refined through a subsequent time-domain model updating.

•
The dynamic measurements can provide more information about the uncertain material parameters compared to the static/pseudo-static responses. Hence, here, the acceleration measurements are used directly in the time-domain model updating.

•
To improve the numerical stability and convergence of the model parameters the linear and nonlinear response behavior of the bridge are assimilated through the two-step model updating process.
To show the two-step FE model updating method and validate its applicability for damage identification in a real-world setting, a pair of full-scale precast prestressed bridge I-girders were used as testbed structures. These girders were in service from 1971 until 2009 before they were decommissioned and repurposed for research experiments [45]. A series of forced-vibration experiments were designed specifically for this study. The girders were subjected to sinusoidal force excitations, and their acceleration responses were measured at different locations. First, the collected acceleration responses are used to identify the modal signatures of the testbed structure. Then, the two-step FE model updating is carried out. In the first step, the initial FE model of girders is updated in the modal domain, and boundary conditions, including stiffness of supports and coupling beams, as well as material properties, including initial stiffness of concrete material, are estimated. The updated model is used as the prior model in the Bayesian model updating process to estimate concrete compressive strength and damping properties. Comparison between the posterior FE-predicted responses and field measurements shows a good agreement in the time domain. Moreover, the final estimates of concrete compressive strength result in a realistic damage identification/quantification for the girders. This process validates the applicability of the introduced two-step FE model updating approach for damage identification of bridge structures/components under operational conditions. While the input load used in this study varies from moving traffic load, this study proves the concept for future real-world application.
The paper is organized as follows. Section 2 is focused on test methodology and preliminary results including an introduction to the testbed structure and the experiments, modal identification, and development of the initial FE model. The two-step FE model updating method and the results are discussed in Section 3. Concluding remarks are provided in Section 4.

Material, Test Methodology, and Preliminary Results
In this section, first, a description of the field experiment including testbed structure, dynamic excitation system, wireless sensing network, and force-vibration tests is presented. Then, the modal identification process and corresponding results are shown in Section 2.2. Finally, the initial FE model is developed in Section 2.3. The testbed structure in this study includes two AASHTO precast prestressed bridge I-girders that were part of the Maryland 90 bridge. After decommissioning, the girders were salvaged and transferred to the Turner-Fairbank Highway Research Center (TFHRC) in McLean, VA, USA to be used as a research testbed [45,46]. The cross-section and elevation views of the girders are shown in Figure 1. presented. Then, the modal identification process and corresponding results are shown in Section 2.2. Finally, the initial FE model is developed in Section 2.3.

Testbed Structure
The testbed structure in this study includes two AASHTO precast prestressed bridge I-girders that were part of the Maryland 90 bridge. After decommissioning, the girders were salvaged and transferred to the Turner-Fairbank Highway Research Center (TFHRC) in McLean, VA, USA to be used as a research testbed [45,46]. The cross-section and elevation views of the girders are shown in Figure 1. The studied girders are 1.37 m deep and 26 m long. Reinforcements and prestressing steel strands are shown in Figure 1. A side view of the testbed structure (with a 20 cm slab on top of each girder) is shown in Figure 2a. After being transferred to the TFHRC, each girder was placed on two 1 × 1 m 2 bearing pads on top of two 2.3 × 1.6 × 1.6 m 3 geosynthetic reinforced soil (GRS) piers [47]. These can be seen in Figure 2b,c. The girders were placed parallel to each other with 2.9 m centerline spacing and were connected with four coupling beams with 0.3 × 1.4 m 2 cross-sectional area. Three out of four coupling beams are seeable in Figure 2d. The studied girders are 1.37 m deep and 26 m long. Reinforcements and prestressing steel strands are shown in Figure 1. A side view of the testbed structure (with a 20 cm slab on top of each girder) is shown in Figure 2a. After being transferred to the TFHRC, each girder was placed on two 1 × 1 m 2 bearing pads on top of two 2.3 × 1.6 × 1.6 m 3 geosynthetic reinforced soil (GRS) piers [47]. These can be seen in Figure 2b,c. The girders were placed parallel to each other with 2.9 m centerline spacing and were connected with four coupling beams with 0.3 × 1.4 m 2 cross-sectional area. Three out of four coupling beams are seeable in Figure 2d.
The studied girders were in service in a corrosive environment for almost 40 years. The environment simultaneously exposed the concrete matrix of the girders to physical and chemical deterioration processes. Concrete delamination and degradation, as well as steel corrosion, are the main damage mechanisms for concrete bridges in such environment [48]. Due to this, the girders experienced aging and deterioration in several locations, including cracking, steel reinforcement corrosion, spalling, etc. Figure 2e shows an example of the observed damage in girders.
Aside from the aging-related damage discussed above, in 2012, salt spray chambers were installed on each girder to accelerate deterioration in the girders. The chamber installed on the west girder sprayed a 15 weight percent (wt.%) NaCl solution and the chamber installed on the east girder sprayed a 3.5 wt.% NaCl solution. This was part of a study to develop protocols for non-destructive testing (NDT) methods for prestressed girder bridges [45]. The studied girders were in service in a corrosive environment for almost 40 years. The environment simultaneously exposed the concrete matrix of the girders to physical and chemical deterioration processes. Concrete delamination and degradation, as well as steel corrosion, are the main damage mechanisms for concrete bridges in such environment [48]. Due to this, the girders experienced aging and deterioration in several locations, including cracking, steel reinforcement corrosion, spalling, etc. Figure 2e shows an example of the observed damage in girders.
Aside from the aging-related damage discussed above, in 2012, salt spray chambers were installed on each girder to accelerate deterioration in the girders. The chamber installed on the west girder sprayed a 15 weight percent (wt.%) NaCl solution and the chamber installed on the east girder sprayed a 3.5 wt.% NaCl solution. This was part of a study to develop protocols for non-destructive testing (NDT) methods for prestressed girder bridges [45].

Dynamic Excitation System
The shaker used in this study was a small uniaxial hydraulic device capable of producing arbitrary displacement motions in the vertical direction. The device consisted of a hydraulic piston connected to a servo-hydraulic valve that controlled the motion of a stack of steel plates that combined to form 4450 N of moving weight. The motion was controlled by an MTS PID hydraulic controller using an LVDT sensor to provide feedback displacement. The shaker's hydraulic piston, moving weights, and steel frame of the shaker were supported on four 4450 N load cells, which were used to measure the total force generated by the shaker. The load cells were installed between the shaker plate and the clamping plate on the girder. The shaker plate with dimensions of 30 cm × 30 cm × 2 cm was located on the top center of the clamping plate with dimensions of 140 cm × 70 cm × 5 cm (see Figure 3). The shaker displacement and force time histories were collected through a LORD-Microstrain V-Link-200 wireless node [49]. The wireless sensing network is discussed further in the following section. Figure 3 shows a close view of the shaker, and Figure 4 Buildings 2023, 13, 420 6 of 27 shows the two locations (layout 1 and layout 2) at which the shaker was installed on the testbed structure. The hydraulic power for the shaker was provided by a diesel-powered mobile pump.
by the shaker. The load cells were installed between the shaker plate and the clam plate on the girder. The shaker plate with dimensions of 30 cm × 30 cm × 2 cm wa cated on the top center of the clamping plate with dimensions of 140 cm × 70 cm × (see Figure 3). The shaker displacement and force time histories were collected throu LORD-Microstrain V-Link-200 wireless node [49]. The wireless sensing network is cussed further in the following section. Figure 3 shows a close view of the shaker, Figure 4 shows the two locations (layout 1 and layout 2) at which the shaker was insta on the testbed structure. The hydraulic power for the shaker was provided by a di powered mobile pump.

Wireless Sensing Network
The wireless sensing network included ten battery-powered triaxial MEMS wir accelerometers (LORD-Microstrain G-Link-200) and a V-Link-200 node [49]. E ment. The shaker's hydraulic piston, moving weights, and steel frame of the shaker were supported on four 4450 N load cells, which were used to measure the total force generated by the shaker. The load cells were installed between the shaker plate and the clamping plate on the girder. The shaker plate with dimensions of 30 cm × 30 cm × 2 cm was located on the top center of the clamping plate with dimensions of 140 cm × 70 cm × 5 cm (see Figure 3). The shaker displacement and force time histories were collected through a LORD-Microstrain V-Link-200 wireless node [49]. The wireless sensing network is discussed further in the following section. Figure 3 shows a close view of the shaker, and Figure 4 shows the two locations (layout 1 and layout 2) at which the shaker was installed on the testbed structure. The hydraulic power for the shaker was provided by a dieselpowered mobile pump.

Wireless Sensing Network
The wireless sensing network included ten battery-powered triaxial MEMS wireless accelerometers (LORD-Microstrain G-Link-200) and a V-Link-200 node [49]. Each

Wireless Sensing Network
The wireless sensing network included ten battery-powered triaxial MEMS wireless accelerometers (LORD-Microstrain G-Link-200) and a V-Link-200 node [49]. Each accelerometer had an adjustable measurement range of ±8 g and could be configured for continuous, periodic, or event-triggered sampling modes to output acceleration, tilt, or derived vibration parameters (velocity, amplitude, etc.). The measured data could be transmitted in real-time and/or be saved to the onboard memory with storage capacity up to 8 × 10 6 data points. The accelerometers had a noise density of 25 µg √ Hz with a wireless range of up to 1 km and an adjustable sampling rate of up to 4 kHz. Each accelerometer had dimensions of 47 mm × 43 mm × 44 mm. To install the accelerometers, zinc-plated steel washers were glued on top surface of the girders. Then, each accelerometer was screwed to a magnetic base and attached to the washers. Figure 5 shows an installed wireless accelerometer. Data collection and coordination between the wireless nodes including the accelerometers and the V-Link-200 node, which was used to collect shaker data, were carried out through the wireless USB data acquisition gateway. The gateway used the lossless extended range synchronized (LXRS) data communication protocol and facilitated lossless data collection with node synchronization of ±50 µs. Synchronization was carried out by transmission of a continuous system-wide timing reference known as the beacon. The communication between the gateway and sensors was wireless through a license-free 2.405 GHz to 2.480 GHz radio frequency with 16 channels. The configuration of the network, data acquisition initialization, and sampling mode selection were managed through the SensorConnect software [50], which was installed on a host computer. The layout of the employed accelerometers is shown in Figure 4. In this study, the sampling rate was 128 Hz, and the acceleration data were collected in directions 1 (i.e., east-west) and 3 (i.e., up-down).
steel washers were glued on top surface of the girders. Then, each accelerome screwed to a magnetic base and attached to the washers. Figure 5 shows an install less accelerometer. Data collection and coordination between the wireless nodes in the accelerometers and the V-Link-200 node, which was used to collect shaker da carried out through the wireless USB data acquisition gateway. The gateway u lossless extended range synchronized (LXRS) data communication protocol an tated lossless data collection with node synchronization of ±50 μs. Synchroniza carried out by transmission of a continuous system-wide timing reference know beacon. The communication between the gateway and sensors was wireless th license-free 2.405 GHz to 2.480 GHz radio frequency with 16 channels. The config of the network, data acquisition initialization, and sampling mode selection we aged through the SensorConnect software [50], which was installed on a host co The layout of the employed accelerometers is shown in Figure 4. In this study, t pling rate was 128 Hz, and the acceleration data were collected in directions 1 (i west) and 3 (i.e., up-down).

Forced-Vibration Testing
The testbed structure was subjected to a series of designed sinusoidal forc tions through frequency sweeps with pre-defined frequency range, duration, and tude. Sweeps ranged between 2 Hz and 20 Hz, including 50 different frequencies ing logarithmically with a duration of 30 s for each frequency. Moreover, the sw cited the girders with three different target load amplitudes equal to 445 N, 2225 4450 N. Considering two layouts (see Figure 4) and three levels of load amplitud ers were tested under six frequency sweeps. A general view of the testbed struct ing the field experiments is shown in Figure 6. The excitation force time history stantaneous excitation frequency-calculated using a short-time Fourier transfor for each sweep are presented in Figure 7.

Forced-Vibration Testing
The testbed structure was subjected to a series of designed sinusoidal force excitations through frequency sweeps with pre-defined frequency range, duration, and amplitude. Sweeps ranged between 2 Hz and 20 Hz, including 50 different frequencies increasing logarithmically with a duration of 30 s for each frequency. Moreover, the sweeps excited the girders with three different target load amplitudes equal to 445 N, 2225 N, and 4450 N. Considering two layouts (see Figure 4) and three levels of load amplitudes, girders were tested under six frequency sweeps. A general view of the testbed structure during the field experiments is shown in Figure 6. The excitation force time history and instantaneous excitation frequency-calculated using a short-time Fourier transform [51]-for each sweep are presented in Figure 7.

Modal Identification
In this section, the modal properties of the testbed structure are identified from the forced-vibration experimental data. For the purpose of modal identification, the collected data from the sweep with layout 1 and 445 N target load amplitude are used (see Figure 7a). These data include the acceleration of 10 channels in directions 1 and 3 (results in 20 input signals-see Figure 4) and the shaker excitation force. A brief summary of the modal identification process and the identified modal properties are presented in this section.
Various modal identification techniques are available in the literature to identify modal properties-including natural frequencies, damping ratios, and mode shapes-from experimental data [19,21,22]. In this study, due to the nature of the excitation, the empirical frequency response functions (EFRFs) [22] are calculated using applied input (shaker excitation force) and measured outputs (acceleration responses). Then, the calculated multioutput EFRFs are used to estimate a state-space model. This two-stage frequency-domain approach helps to specify the frequency band of interest in which the modal identification needs to be carried out.
The EFRF between a measurement signal y(t) and input signal u(t)-depicted by G( f )-is defined as below: where Y( f ) and U( f ) are Fourier transforms of y(t) and u(t), respectively. In a realworld setting, data are always polluted with measurement noise. To reduce the effects of the measurement noise and obtain a smooth EFRF, Welch's averaging method [52] with 12,800 data points Hamming window is used for spectral estimation. This size of window is selected to ensure covering the full length of excitation and the ambient signal before and/or after it.
Having EFRFs for all 20 signals, an n-order state-space model is estimated to fit the estimated EFRFs in the frequency band of interest. In this study, to reduce uncertainties due to low-and high-frequency noises, the frequency band of interest is selected between 2 to 25 Hz. Blue curves in Figure 8 show the calculated EFRFs using measured data.
To estimate the state-space model, the subspace state-space identification method is used [53,54]. While this method is briefly explained here, the proof of theory and more details can be found in [53]. In the subspace state-space identification method, measurements are placed in a block Hankel matrix which is divided into a past and a future part. The identification algorithm proceeds with projecting the future measurements into the past measurements, while the projection matrix can be factorized as the product of an observability matrix and a state sequence. These two matrices are identified by applying the singular value decomposition (SVD) to the projection matrix, and the order of the system is calculated as the number of non-zero singular values. By applying one block shift in separation between past and future measurements in the Hankel matrix, another projection matrix, shifted observability, and state sequence matrices can be obtained. At this point, the system matrices can be calculated from the overdetermined set of linear equations.
The numerical algorithms for system identifications are available in the Matlab n4sid [55] function and are used in this study. As a classical remedy, the modal identification is carried out for a range of model orders, and a stability diagram is plotted on which true modes appear as stable modes [56]. For this purpose, the stability analysis is run considering model orders from 2 to 40 with 1% and 5% error tolerances for natural frequency and damping ratios, respectively. The stability analysis showed that a model order of n = 26 is the lowest model order to have all stable modes within the frequency band of interest. The fits between estimated and calculated EFRFs can be improved using the prediction error minimization algorithm and nonlinear least-squares objective functions. This approach is carried out using the Matlab ssest function [57], which initializes the model parameters based on the previously estimated state-space model, and then updates the parameters using an iterative search to minimize the prediction errors [21]. Red curves in Figure 8 display EFRFs of the estimated state-space model at measurement points. As can be seen, the estimated state-space model is able to approximate the calculated EFRFs acceptably. Identified natural frequencies ( f ID ) and damping ratios (ξ ID ) of the system are reported in Table 1. The identified mode shapes-those which will be used for modal-based model updating-are later shown in Section 3.1. To estimate the state-space model, the subspace state-space identification method is used [53,54]. While this method is briefly explained here, the proof of theory and more details can be found in [53]. In the subspace state-space identification method, measurements are placed in a block Hankel matrix which is divided into a past and a future part. The identification algorithm proceeds with projecting the future measurements into the past measurements, while the projection matrix can be factorized as the product of an observability matrix and a state sequence. These two matrices are identified by applying the singular value decomposition (SVD) to the projection matrix, and the order of the system is calculated as the number of non-zero singular values. By applying one block shift in separation between past and future measurements in the Hankel matrix, another projection matrix, shifted observability, and state sequence matrices can be obtained. At this point, the system matrices can be calculated from the overdetermined set of linear equations.

Finite Element Modeling of the Testbed Structure
The initial FE model of the testbed structure is developed in OpenSees [58] using the available as-built drawings. The girders are modeled using fiber-section force-based beamcolumn elements (forceBeamColumn) with approximately 30 cm length and 3 integration points. Moreover, the linear-elastic shear stiffness and torsional stiffness are aggregated to the fiber sections. The slab is also modeled using rectangular fiber sections as a part of the girders' sections. Steel reinforcement is modeled using single fibers located at the center of each bar. To model the profile of the draped strands-with a nominal cross-sectional area of 1 cm 2 -the strands' depth at the integration points of each element is calculated based on Figure 1. The strand is modeled using a single fiber at its corresponding depth. Concrete material for girders and slabs is modeled with nominal compressive strength ( f c,E and f c,W for the east and west girders) of 46 MPa, strain at maximum strength (ε c ) of 0.2%, strain at crushing strength (ε u ) of 0.5%, and zero crushing strength. As girders are already damaged, no tensile strength is assumed, and therefore, the Concrete01 material model is employed. ε c ) which is 46 GPa here. Steel is modeled using bilinear steel01 material with a modulus of elasticity of 200 GPa, and yielding strength of 455 MPa and 1720 MPa for reinforcing steels and strands, respectively. Moreover, the prestressing force in strands (128.60 kN per strand) is modeled by applying its resulting initial strain to the steel material. For this purpose, the normal strain in strands' steel resulting from the prestressing force is calculated by dividing the prestressing force by strand's axial rigidity (product of the strand's modulus of elasticity and gross section area). The calculated strain (0.65%) is assigned to the strands' steel01 material using InitStrainMaterial. The mass, 104 × 10 3 kg in total, and weight of the girders are assigned to the element nodes.
Coupling beams are modeled using elasticBeamColumn elements considering the gross cross-sectional area and an Elastic material with the modulus of elasticity equal to 30 GPa. The coupling beams are connected to the girders assuming rigid connections. The mass and weight of the coupling beams are assigned to the element nodes. To model each support, the girder nodes that are located along the bearing pad are constrained to a node at the center of the bearing pad using rigidLink constraint. To account for flexibility in the piers and bearing pads, supports are defined using a 6 degrees of freedom (DOFs) spring. This is modeled using ZeroLength elements with Elastic uniaxial material. Moreover, the energy dissipation in the piers and bearing pads is collectively modeled using a dashpot in the vertical direction. This is performed using ZeroLength elements with Viscous uniaxial material. The corresponding stiffness and damping parameters are initially selected based on engineering judgment [47,59] and are later updated using modal-based and timedomain model updating. In summary, the FE model consists of 214 beam-column elements, 16 rigidLink elements, 12 ZeroLength elements, and 230 joints. The list of model parameters that are later treated as unknown parameters in Section 3.1 and 3.2 and their initial values are summarized in Table 2. In this table, the equivalent values for concrete compressive strength and modulus of elasticities of girders are shown in the same row separated with '/'. The directions in this table are based on Figure 4. The nonlinear time history analysis is performed using the Newmark average acceleration method with a constant time step size of 0.0078 s equal to the measurement sampling rate. The Newton-Raphson method is used to iteratively solve the nonlinear equilibrium equations [60]. To define energy dissipation in the structural system aside from the material nonlinearity, modal damping is modeled for the first six modes. The damping ratios, i.e., ξ i , ∀ i ∈ {1, 2, 3, 4, 5, 6}, are set equal to the identified ones (see Table 1). The only exceptions are damping ratios for modes 1 and 2, which are set equal to 0.02. The reason is to have similar initial values-which is not very different from the identified values (0.0192 and 0.0178)-during time-domain model updating, which is discussed later.

First Step: Modal-Based Model Updating
The modal properties of the initial FE model are different from the identified ones ( Figure 9). Hence, the initial FE model needs to be updated using modal-based model updating to better fit the identified modal properties. As mentioned before, the modalbased model updating is limited to the linear response behavior of structures. Hence, only the linear model parameters are updated at this step. Modal-based model updating process is discussed next.
The modal-based model updating is a process to minimize the discrepancies between identified and FE-predicted modal frequencies and mode shapes by updating the linear parameters of the FE model [24,61]. In this process, an objective function, g(θ), is defined as shown in Equation (2). The discrepancies between modal frequencies, mode shapes, and a regularization term are respectively the first, second, and third terms in Equation (2).
In Equation (2), the parameter W f is the weighing scalar for frequency residuals and the term r f ∈ R N×1 is the vector including the square root of normalized differences between FE-predicted and identified modal frequencies. The parameter W M is the weighing scalar for modal assurance criteria (MAC) residuals. The term N is the total number of identified modes that are used for model updating, and the term MAC i indicates the MAC value for mode i. The parameter W θ is the weighing scalar for penalizing large deviation of the unknown FE model parameters from their initial values. The vector θ ∈ R n θM ×1 is the vector of unknown FE model parameters and n θM is the number of unknown model parameters for modal-based model updating. The vector θ 0 ∈ R n θ M ×1 is the initial estimates of the unknown model parameters. The i th entry of the vector r f , denoted as r f ,i , is defined as follows: The modal-based model updating is a process to minimize the discrepancies between identified and FE-predicted modal frequencies and mode shapes by updating the linear parameters of the FE model [24,61]. In this process, an objective function, ( ), is defined as shown in Equation (2). The discrepancies between modal frequencies, mode shapes, and a regularization term are respectively the first, second, and third terms in Equation (2).
In Equation (2), the parameter is the weighing scalar for frequency residuals and the term ∈ ℝ × is the vector including the square root of normalized differences between FE-predicted and identified modal frequencies. The parameter is the weighing scalar for modal assurance criteria (MAC) residuals. The term N is the total number of identified modes that are used for model updating, and the term indicates the MAC value for mode . The parameter is the weighing scalar for penalizing large deviation of the unknown FE model parameters from their initial values. The vector ∈ ℝ ×1 is the vector of unknown FE model parameters and is the number of unknown model parameters for modal-based model updating. The vector 0 ∈ ℝ ×1 is the initial estimates of the unknown model parameters. The ℎ entry of the vector , denoted as , , is defined as follows: In the above equation, the terms and are the modal frequencies of the ℎ identified and FE-predicted modes. The term is defined as follows: The terms ψ ID i and ψ FE i ∈ R NDOF×1 are the i th identified and FE-predicted mode shape vectors, respectively, and NDOF is the number of DOFs. The superscript T denotes matrix/vector transpose operator.
To update the FE model using modal-based model updating, the stiffness-related model parameters E (= E E = E W ), E c , K R , K V and K T are selected as unknown FE model parameters to be updated (n θM = 5). Note that parameter K L is not selected as there is no measurement in the longitudinal direction of the girders. The vector θ 0 is initiated using the initial values listed in Table 2.
The first two lateral and the first two vertical identified modes (N = 4) are used for the modal-based model updating. This is because the first two identified vertical modes are the only ones with frequencies less than 20 Hz (maximum excitation frequency). As can be seen in Figure 4, measurements are collected in 5 vertical (i.e., in direction 3) and 5 lateral (i.e., in direction 1) DOFs for each girder; therefore, NDOF = 20. The terms W f , W M and W θ are selected equal to 4, 1, and 0.01 to balance the contributions of the MAC value and frequency errors in the objective function. As a result, the difference between identified and estimated modal frequencies will be penalized more than mode shapes [61]. This is because the uncertainty in identified mode shapes is greater than frequencies. The constrained nonlinear multivariate optimization function fmincon within the Matlab Optimization Toolbox is used [62] with the interior-point algorithm [63]. The maximum iteration number is set to 30, and the process terminates when the relative difference between successive values of the objective function is lower than a given threshold of 10 −4 . Moreover, the lower and upper bounds for parameter estimation are set equal to 10 ±4 times of the initial values.
A comparison between the identified, initial, and updated modal properties is shown in Figure 9, and the updated values for the FE model parameters are shown in Table 3. It is noteworthy that as parameters E E = E W are updated, the parameters f c,E = f c,W are also updated (using E E = 2 f c,E ε c and E W = 2 f c,W ε c assuming fixed value for ε c ). These equivalent values for each girder are shown in a same column in Table 3 and are separated using '/'. A comparison between Tables 2 and 3 shows that the parameters E E and E W are estimated to a smaller value than their initial ones. This was expected as both girders are aged and have experienced severe damage. The parameter E c is estimated to a smaller value than its initial one. This is probably due to the presence of cracks in the section and could also be an indication of the fact that the connections between the girders and coupling beams are not completely rigid. The parameters K R , K V , and K T are estimated to have values greater than their initial ones, which indicates that the bearing pads are stiffer than what is considered in the initial model. The improvement in the modal frequencies is superior to the MAC values. This means that the modal-based model updating process compensates for the frequency match with MAC values. The maximum error in frequency is at the second lateral mode, which is less than 8% and is acceptable. The updated model is later used as the prior FE model for the time-domain model updating.

Second Step: Time-Domain Model Updating
To this point, the initial FE model of the studied girders is updated using modal-based model updating. It is noteworthy that concrete material has nonlinear response behavior even under small levels of excitation. In addition to that, presence of prestressing force pushes the concrete material across the section along its nonlinear response curve. Then, the applied excitation-shaker force here and traffic load in operational conditions-results in small loading/unloading of concrete material in nonlinear range of the response curve. The level of nonlinearity in the response behavior of concrete material increases as a function of deterioration and damage [64]. However, the nonlinear response behavior of the studied girders is not captured through modal-based model updating.
To account for the nonlinear response behavior of girders as well as refining the estimation of damage-related model parameters, time-domain model updating-here referred to as the second step of the model updating procedure-is carried out. The cumulative damage effects in a reinforced concrete section can be modeled by altering the stress-strain response behavior of the concrete material, e.g., reduction of the effective compressive strength [65][66][67]. Based on this, it is intended to estimate the concrete compressive strengths of girders using time-domain model updating. Moreover, the acceleration measurements contain information regarding the dynamic behavior of the testbed structure. Hence, the energy-dissipation-related model parameters can also be estimated using time-domain model updating. For this purpose, first, the most identifiable FE model parameters are selected using an information-theoretic identifiability analysis [68]. Then, these parameters are updated using the Bayesian model updating process.

Identifiability Analysis for Time-Domain Model Updating
The identifiability analysis is an approach to determine the most identifiable unknown model parameters using the sensitivity of FE-predicted responses with respect to the unknown model parameters. In this method, the relative information that each candidate model parameter gains from the responses and the mutual information gain between these parameters are calculated. The model parameters with high relative information gain and little dependency on other parameters are likely more identifiable than others and are selected to be updated in the model updating process.
The data used for identifiability analysis are noted as Data 1-2 in Table 4, which corresponds to layout 1 and the target load amplitude of 2225 N. It is noteworthy that no significant difference is expected in the identifiability analysis results using different layouts and target load amplitudes. However, the target load amplitude for identifiability analysis is selected high enough to have a moderate level of loading/unloading response behavior in the concrete material of the testbed structure. The data in Table 4 are later used for model updating. As can be seen in Table 4, the identifiability analysis and model updating process are performed using experiments with excitation frequencies of 9.73 Hz and 12.86 Hz. These excitation frequencies are selected as they are the closest ones to the identified modal frequencies (9.73 Hz and 12.79 Hz) and are expected to provide useful information on the dynamic behavior of the testbed structure. Moreover, the measured signal-to-noise ratio in the experiments with these frequencies is higher than the similar ratio in experiments with excitation frequencies far from the identified modal frequencies. A higher measured signal-to-noise ratio results in more stable parameter estimation [69].
Each data set in Table 4 augments two experiments with excitation frequencies of 9.73 Hz and 12.86 Hz (while the layout and target load amplitude are similar). This is shown in the following equations: The terms y ∈ R 20×t e and u ∈ R 2×t e denote the acceleration measurements and input excitations that are being used for identifiability analysis and time-domain model updating. The terms y 9.73 Hz ∈ R 10×t e and y 12.86 Hz ∈ R 10×t e refer to the collected acceleration measurements (at 10 measurement channels in direction 3) from experiments with input excitation frequencies of 9.73 Hz and 12.86 Hz. The terms u 9.73 Hz ∈ R 1×t e and u 12.86 Hz ∈ R 1×t e refer to the input excitations with frequencies of 9.73 Hz and 12.86 Hz. The term t e is the total length of collected data.
One of the intentions of identifiability analysis/time-domain model updating is to analyze the identifiability of/update the unknown model parameters ξ 1 and ξ 2 . For this purpose, it is required to implement measurements from both excitation frequencies equal to 9.73 Hz and 12.86 Hz. The augmenting process provides measurements that contain dynamic response behavior of the testbed from both its vertical modes.
Parameters f c,W , f c,E , E c , K R , K V , K L , C D , ξ 1 , and ξ 2 are selected as candidate unknown model parameters. It is intended to refine the damage state estimation of girders for the west and east girder separately. Therefore, the concrete compressive strengths of the west and east girders are modeled using two different parameters f c,W and f c,E , respectively. While the main intention of the time-domain model updating is to update the model parameters related to the energy dissipation and nonlinear behavior of the structure, the parameters E c , K R , K V , and K L are also considered as candidate unknown parameters. As will be seen later, this is to highlight how the introduced two-step model updating helps with unidentifiability and mutual dependencies between the model parameters. Since the input load mainly excites the first two vertical modes of the girders, the modal damping ratios of these two modes are included in the list of candidates for unknown parameters.
The value for parameters f c,W , f c,E , E c , K R , and K V are set equal to their final estimates in modal-based model updating (see Table 3). As the input load is in the vertical direction, the parameter K T is most probably not identifiable and is excluded from the list of candidate model parameters. The value of parameter K L is set equal to the updated value of K T (after modal-based model updating) as the properties of bearing pads are assumed to be the same in the transverse and longitudinal directions. The value of parameters C D , ξ 1 , and ξ 2 are set equal to their initial values (see Table 2), as these parameters were not updated using modal-based model updating.
The relative information gain of the candidate unknown model parameters and relative mutual information gain between the candidate unknown model parameter pairs are shown in Figure 10. As can be seen in this figure, although the parameters E c , K R , K V , and K L have considerable levels of relative information gain, they are highly dependent on the parameters f c,W and f c,E . This dependency is likely because all these parameters contribute to the stiffness of the testbed structure. Based on this, parameters E c , K R , K V , and K L are most probably not identifiable together with parameters f c,W and f c,E . Parameters E c , K R , K V , and K L reflects the linear-elastic response behavior of the testbed structure and are already estimated using modal signatures. Moreover, the parameters f c,W and f c,E have relatively large information gain, and their estimation is of main interest in time-domain model updating as their final estimates help to refine the damage estimation in girder level and reflect the cumulative damage status of each girder. Hence, parameters f c,W and f c,E are selected to be estimated using time-domain model updating while parameters E c , K R , K V , and K L are fixed at their corresponding values in Table 3 obtained from the modal-based model updating. This reduces the challenges of model updating due to the unidentifiability and/or mutual dependencies between model parameters. It is noteworthy that the final estimates of parameters f c,W and f c,E using time-domain model updating will inherently be dependent on the fixed values selected for parameters E c , K R , K V , and K L . Parameters ξ 1 , ξ 2 , and C D have relatively moderate information gain and have a negligible dependency on other parameters. However, parameter C D is dependent on parameters ξ 1 and ξ 2 as all of them contribute to the viscous damping energy dissipation of the structure. As the initial value for parameter C D is selected based on judgment, it is of interest to update this parameter using the model updating process. However, the final estimates of parameters ξ 1 , ξ 2 , and C D are expected to vary between different case studies and depend on each other. Therefore, the overall damping of the testbed will be calculated at the end. Based on the above discussion, parameters f c,W , f c,E , ξ 1 , ξ 2 , and C D are selected to be estimated using the time-domain model updating.

Bayesian Inference
In the Bayesian model updating, the unknown FE model parameter vector is considered as a random vector with joint probability density function (PDF) whose mean (referred to as estimate hereafter) and covariance are updated recursively through the integration of the FE-predicted and measured responses. The unknown FE model parameters are updated using measured responses in successive overlapping windows. In this approach, known as the sequential estimation window approach [70], the a priori estimates of the unknown FE model parameters at each estimation window are updated to the a posteriori estimates. A brief review of this method is provided in the following. 1 2 mation gain and have a negligible dependency on other parameters. However, parameter is dependent on parameters 1 and 2 as all of them contribute to the viscous damping energy dissipation of the structure. As the initial value for parameter is selected based on judgment, it is of interest to update this parameter using the model updating process. However, the final estimates of parameters 1 , 2 , and are expected to vary between different case studies and depend on each other. Therefore, the overall damping of the testbed will be calculated at the end. Based on the above discussion, parameters , ′ , , ′ , 1 , 2 , and are selected to be estimated using the time-domain model updating.
(a) (b) Figure 10. Identifiability analysis results: (a) relative information gain of the candidate unknown model parameters; (b) relative mutual information gain between the candidate unknown model parameter pairs. Each row is normalized to its diagonal value. Then, the diagonal values are nullified.

Bayesian Inference
In the Bayesian model updating, the unknown FE model parameter vector is considered as a random vector with joint probability density function (PDF) whose mean (referred to as estimate hereafter) and covariance are updated recursively through the integration of the FE-predicted and measured responses. The unknown FE model parameters are updated using measured responses in successive overlapping windows. In this approach, known as the sequential estimation window approach [70], the a priori estimates of the unknown FE model parameters at each estimation window are updated to the a posteriori estimates. A brief review of this method is provided in the following.
Assuming zero-mean Gaussian white noise processes for the measurement noise ( ) and the process noise ( ), the state-space model at each estimation window is set up as follows: Figure 10. Identifiability analysis results: (a) relative information gain of the candidate unknown model parameters; (b) relative mutual information gain between the candidate unknown model parameter pairs. Each row is normalized to its diagonal value. Then, the diagonal values are nullified.
Assuming zero-mean Gaussian white noise processes for the measurement noise (λ) and the process noise (q), the state-space model at each estimation window is set up as follows: As can be seen in Equation (7), the unknown FE model parameter vector for Bayesian model updating, ϑ ∈ R n ϑB FE ×1 , evolves linearly in the state equation by a random walk process. The term n ϑB FE is the number of unknown FE model parameters in the Bayesian model updating. In this study, the term n ϑB FE is equal to 5 (which counts the number of parameters f c,W , f c,E , ξ 1 , ξ 2 , and C D ) and the initial values for the unknown model parameters are equal to those described in the previous section. In Equation (8), the term y w ∈ R (n y ×t w )×1 denotes the measurement vector and h w (.) ∈ R (n y ×t w )×1 denotes the nonlinear FE-predicted response function at the w th estimation window, which spans between time steps t w,1 and t w,2 with t o time steps overlap with the previous estimation window. The parameter n y (here 20) is the number of collected signals and t w is number of time steps at the w th estimation window. The term u w ∈ R t w,2 ×1 is the deterministic input vector at the w th estimation window. In this study, estimation windows have the length of 150 time steps with 50 time steps overlap, i.e., t 1,2 = 50 and t w = 150 and t w−1,2 − t w,1 = 50 ∀ w ≥ 2.
The measurement equation (Equation (8)) is linearized using the first-order Taylor series expansion with linearization point at the a priori estimate. The following equation is obtained: in which the superscripts -/+ denote the a priori/posteriori estimates, and derivation of h w (ϑ w , u w ) with respect to ϑ w is referred to as the sensitivity matrix and is calculated using the finite difference approach. The estimation method that is used in this study is referred to as the Extended Kalman filter (EKF) for parameter-only estimation [71]. In this method, the a priori estimates of the mean vector and covariance matrix of the unknown model parameters at each estimation window are considered equal to the a posteriori estimates at the previous estimation window. This is depicted in Equations (10) and (11). Moreover, the estimation problem is solved across each estimation window iteratively, and the subscripts 0 and k in the following equations denote the iteration number.
The termP −/+ ϑϑ,w is the a priori/posteriori estimate of the covariance matrix of the unknown model parameters at the w th estimation window. For the sake of brevity, the details for the derivation of the Bayesian inference based on the EKF method is not shown here, and only the recursive equations at each estimation window are presented. After completing the iteration process at each estimation window, the estimation moves to the next window, and the process is repeated. To complete the iteration process at each estimation window, the estimates of unknown model parameters need to be converged. However, to improve the efficiency of the process, the maximum number of iterations is limited. More details are available in [71] The term Q w ∈ R n ϑB ×n ϑB is the covariance matrix for the process noise (q w ∼ N(0, Q w )). The matrix Q w is a diagonal matrix, and its j th diagonal entry is equal to q times the j th entry in vectorθ − w . The term q is set equal to 0.002. The term R w ∈ R (t w ×n y )×(t w ×n y ) is the measurement noise ( λ w ∼ N(0, R w )) and is modeled as a block diagonal matrix with the simulation error covariance matrix-including measurement noise-on the diagonal blocks. In this study, the diagonal entries of the matrix R w are set equal to (0.32%g) 2 at all measurement channels. The value 0.32%g is approximately equal to the average root-mean-square of ambient measurements. The termŷ − w,k+1 is the a priori FE-predicted response calculated atθ − w,k+1 . The matrix K w,k+1 ∈ R n ϑB ×(t w ×n y ) is the Kalman-gain matrix at the (k + 1) th iteration in the w th estimation window. The matrix C w,k+1 ∈ R (t w ×n y )×n ϑB is the FE response sensitivity matrix-with respect toθ − w,k+1 -at the (k + 1) th iteration in the w th estimation window. The term I ∈ R n ϑB ×n ϑB denotes the identity matrix. The matrixP − yy,w,k+1 ∈ R (t w ×n y )×(t w ×n y ) is a priori estimate of the covariance matrix ofŷ − w,k+1 , andP − ϑy,w,k+1 ∈ R n ϑB ×(t w ×n y ) is the a priori estimate of the cross-covariance matrix ofθ initial variance of the initial estimate of the unknown model parameters. The j th diagonal entry inP + ϑϑ,0 is equal to the square of p ϑ times the j th entry in vectorθ + 0 . In this study p ϑ is set equal to 0.1. The recursive Bayesian model updating process in each estimation window is completed after 10 iterations or meeting the following convergence criteria in the posterior estimates of unknown model parameters.

Results
This section presents the results of the second step of the introduced two-step model updating approach as well as its application for damage identification of the testbed structure. First, the Bayesian model updating is carried out, and its performance is discussed through the updating process of unknown model parameters and the fit between measurements and posterior estimates of FE-predicted responses. Subsequently, the final estimates of unknown model parameters are used to infer damage in the girders and calculate the overall damping of the testbed structure. To assess the application of Bayesian FE model updating, the data sets presented in Table 4 are used. As explained in Section 3.2.1, while the excitation frequencies close to the first two modes are lumped together, the target load amplitudes and layouts differ from one data set to another one. As mentioned before, lumping the frequencies together increases the identifiability of model parameters and the parameter estimation stability by using the dynamic response behavior of testbed structure from its two first vertical modes.
The updating process for the posterior estimates of the unknown model parameters f c,W , f c,E , ξ 1 , ξ 2 , and C D using data from Table 4 is shown in Figure 11. In this figure, the estimates of the unknown model parameters are normalized to their corresponding initial values in time-domain model updating process. As can be seen in this figure, all the unknown model parameters are updated from their initial values and smoothly converged to their final estimates.  Table 4.
Next, the prior estimates of responses-simulated responses using the prior FE model-and the posterior ones are compared with the field measurements. For the sake of brevity, only one second of data obtained using Data 1-2 are shown in Figure 12. In this figure, the two top rows correspond to channels 1-10 and an excitation frequency of 9.73 Hz, and the two bottom rows correspond to the same channels and an excitation frequency of 12.86 Hz. For the case of excitation frequency of 9.73 Hz, the prior FE responses mostly underestimate the measurements in channels located between the supports (channels 2, 3, 4, 7, 8, and 9) and slightly overestimate the responses in channels located on the  0  20  40  10  30  0  20  40  10  30  50 0  20  40  10  30  50   0  20  40  10  30  0  20  40  10  30  50  60 0  20  40  10 30 50 Figure 11. Updating process for the posterior estimates of the unknown model parameters using data in Table 4.
Next, the prior estimates of responses-simulated responses using the prior FE model-and the posterior ones are compared with the field measurements. For the sake of brevity, only one second of data obtained using Data 1-2 are shown in Figure 12. In this figure, the two top rows correspond to channels 1-10 and an excitation frequency of 9.73 Hz, and the two bottom rows correspond to the same channels and an excitation frequency of 12.86 Hz. For the case of excitation frequency of 9.73 Hz, the prior FE responses mostly underestimate the measurements in channels located between the supports (channels 2, 3, 4, 7, 8, and 9) and slightly overestimate the responses in channels located on the overhang parts (channels 1, 5, 6, and 10). However, the responses in all channels are initially overestimated for the excitation frequency of 12.86 Hz. This shows that although the prior FE model matches the main modal signatures of the real structure, it cannot correctly predict the measurements in the time domain. As can be seen in Figure 12, after the application of Bayesian model updating, the updated FE model better fits the measurement responses in the time domain. This improvement is noticeable in various channels and for both excitation frequencies of 9.73 Hz and 12.86 Hz.   To quantify the discrepancies between two signals, the relative root mean square error (RRMSE) is calculated as follows: In the above equation, s i andŝ i denotes the measured and estimated responses at the i th time step. The closer the RRMSE gets to zero, the better signals s andŝ match. The RRMSEs are calculated between the measured responses and their prior/posterior FE-predicted responses and are listed in Table 5. As can be seen in Table 5, RRMSE values are reduced from prior to posterior values for all measurement channels. This shows that the time-domain model updating reduces the discrepancies between FE-predicted and measured responses successfully. However, the highest posterior RRMSEs are due to channels 1, 5, 6, and 10. These channels are located on the overhang parts of the girders (see Figure 4), and the measurements are likely less reliable due to the low signal-to-noise ratio. It is understandable from Figure 12 and Table 5 that the FE-predicted responses at these channels are promisingly updated to match the measurements.
The final estimates of unknown model parameters using all data are shown in Figure 13. The average final estimates of parameters f c,W and f c,E are approximately 30% and 26% less than their initial values in Table 2. Previous studies [72,73] have shown more than 25% reductions in concrete compressive strength of deteriorated concrete structures-e.g., abandoned structures and structures in acidic environments. Considering that the operating environment of girders caused them to be prone to various damage mechanisms (including concrete degradation, as well as steel corrosion as discussed in Section 2.1.1), the final estimates of parameters f c,W and f c,E are reasonable. Moreover, as mentioned in Section 2.1.1, after being decommissioned, the west girder has been under more intense environmental testing conditions than the east girder [45]. This can be understood from the model updating results as the final estimates of parameter f c,W is smaller than f c,E . The only exception is the case study with Data 2-1 in which the final estimate of f c,E is smaller than f c,W . This is most probably an estimation error due to measurement noise, experimental error, etc. Moreover, based on Figure 13 and ignoring the case study with Data 2-1, the maximum difference in the estimation of parameter f c in layouts 1 and 2 is less than 7%, which shows the consistency of Bayesian model updating results regardless of the shaker location and target load amplitude.
from the model updating results as the final estimates of parameter , ′ is smaller than , ′ . The only exception is the case study with Data 2-1 in which the final estimate of , ′ is smaller than , ′ . This is most probably an estimation error due to measurement noise, experimental error, etc. Moreover, based on Figure 13 and ignoring the case study with Data 2-1, the maximum difference in the estimation of parameter ′ in layouts 1 and 2 is less than 7%, which shows the consistency of Bayesian model updating results regardless of the shaker location and target load amplitude. Figure 13. The final estimates of unknown model parameters using data sets in Table 4. The information on each set of data is shown in Table 4. The estimates for parameters , ′ , , ′ and are normalized to their corresponding initial values in time-domain model updating process.
The final estimates for modal damping parameters vary between different data sets. The identifiability analysis showed dependencies between parameters 1 , 2 , and . This dependency is most probably a reason for the discrepancy between the final estimates of these parameters using different data sets. Aside from this, it is known that structural damping roots in various factors, including opening and closing of micro cracks, friction in structural joints, level of vibration, and other sources of energy dissipation, etc. [74,75]. These conditions could be varied from one experiment to another and result in different levels of damping. To have insight into the overall damping of the testbed structure, the modal damping ratios are calculated using the state matrix of the system. For this purpose, the stiffness, mass, and damping matrices of the system are developed for the posterior models. The stiffness (K) matrix is recorded as the current global system matrix (using OpenSees printA command) while a static analysis using LoadControl integrator is carried out, and damping is removed from the model. The mass matrix (M) is calculated based on the current global system matrix while a static analysis using LoadControl integrator and a transient analysis using CentralDifference integrator is carried out, and damping is removed from the model. Using the Central Difference formulation presented in Figure 13. The final estimates of unknown model parameters using data sets in Table 4. The information on each set of data is shown in Table 4. The estimates for parameters f c,W , f c,E and C D are normalized to their corresponding initial values in time-domain model updating process.
The final estimates for modal damping parameters vary between different data sets. The identifiability analysis showed dependencies between parameters ξ 1 , ξ 2 , and C D . This dependency is most probably a reason for the discrepancy between the final estimates of these parameters using different data sets. Aside from this, it is known that structural damping roots in various factors, including opening and closing of micro cracks, friction in structural joints, level of vibration, and other sources of energy dissipation, etc. [74,75]. These conditions could be varied from one experiment to another and result in different levels of damping. To have insight into the overall damping of the testbed structure, the modal damping ratios are calculated using the state matrix of the system. For this purpose, the stiffness, mass, and damping matrices of the system are developed for the posterior models. The stiffness (K) matrix is recorded as the current global system matrix (using OpenSees printA command) while a static analysis using LoadControl integrator is carried out, and damping is removed from the model. The mass matrix (M) is calculated based on the current global system matrix while a static analysis using LoadControl integrator and a transient analysis using CentralDifference integrator is carried out, and damping is removed from the model. Using the Central Difference formulation presented in Equation (22), the mass matrix is equal to the recorded current global system matrix (K CDF ) times ∆t 2 . In Equation (22), the term ∆t is the size of the time step increment and is set to a very small value (here 10 −6 s) to record the current global system matrix with high precision. The term C is the damping matrix.K Including damping in the model and carrying out a transient analysis using Newmark integrator with γ = 0.5 and β = 0.25, the matrixK Newmark can be recorded using printA. The damping matrix is calculated as it is shown in Equation (23) using the Newmark formulation.
The matrices M, C, and K are condensed at the dynamic DOFs, and the state matrix is calculated [74]. Using the state matrix of the posterior models, the modal damping ratios are obtained for all data sets. These damping ratios are reported in Figure 14. As can be seen, the testbed structure dissipates a greater level of energy in its first mode than in the second mode, which agrees with the identified modal damping ratios (see Table 1). Moreover, the overall modal damping ratios estimated from experiments with larger target load amplitudes are smaller than those estimated from experiments with smaller target load amplitudes. Although this has been previously observed in a few system identification studies [76,77], it is in contrast with the literature in which damping increases as load amplitude increases. = ∆ (̂− − (∆ ) 2 ) (23 The matrices M, C, and K are condensed at the dynamic DOFs, and the state matri is calculated [74]. Using the state matrix of the posterior models, the modal damping ratio are obtained for all data sets. These damping ratios are reported in Figure 14. As can b seen, the testbed structure dissipates a greater level of energy in its first mode than in th second mode, which agrees with the identified modal damping ratios (see Table 1). More over, the overall modal damping ratios estimated from experiments with larger targe load amplitudes are smaller than those estimated from experiments with smaller targe load amplitudes. Although this has been previously observed in a few system identifica tion studies [76,77], it is in contrast with the literature in which damping increases as load amplitude increases. Figure 14. Overall damping ratio of the testbed structure for mode 1 and mode 2 using differen data sets.

Conclusions
This paper proposed a two-step FE model updating process for operational health monitoring and damage identification of bridge structures. In the proposed approach first, modal-based model updating was carried out to calibrate the initial FE model of th bridge. In this step, the stiffness-related model parameters-mainly related to boundary conditions-as well as the initial stiffness of concrete material were updated to fit the FE predicted and identified modal signatures of the bridge. Then, to account for the nonlinea response behavior of the bridge as well as refinement of model parameter estimation Bayesian time-domain model updating was carried out in the second step. In this step material properties that reflect the cumulative damage in the bridge, e.g., effective con crete compressive strength, as well as damping energy-dissipation-related model param eters, were estimated. The linear-elastic boundary conditions were fixed at their final es timates obtained from the modal-based model updating. To prevent convergence of th model updating algorithm to the local solution, the initial values for concrete compressiv strength were selected using the final estimates of concrete initial stiffness from modal based model updating.

Conclusions
This paper proposed a two-step FE model updating process for operational health monitoring and damage identification of bridge structures. In the proposed approach, first, modal-based model updating was carried out to calibrate the initial FE model of the bridge. In this step, the stiffness-related model parameters-mainly related to boundary conditions-as well as the initial stiffness of concrete material were updated to fit the FE-predicted and identified modal signatures of the bridge. Then, to account for the nonlinear response behavior of the bridge as well as refinement of model parameter estimation, Bayesian time-domain model updating was carried out in the second step. In this step, material properties that reflect the cumulative damage in the bridge, e.g., effective concrete compressive strength, as well as damping energy-dissipation-related model parameters, were estimated. The linear-elastic boundary conditions were fixed at their final estimates obtained from the modal-based model updating. To prevent convergence of the model updating algorithm to the local solution, the initial values for concrete compressive strength were selected using the final estimates of concrete initial stiffness from modal-based model updating.
The application of the two-step model updating approach was presented using a pair of full-scale precast prestressed deteriorated bridge I-girders as the testbed structure. For this study, a series of forced-vibration experiments were planned, and the testbed structure was subjected to sinusoidal force excitations through frequency sweeps at three different amplitudesusing a small shaker. The input excitation was measured using load cells, and the acceleration responses were collected using a wireless sensing network. The findings confirm the following conclusions: • Identifiability analysis showed significant mutual dependency between different model parameters. This mutual dependency could lead to weak identifiability of model parameters in the traditional FE model updating process. The proposed twostep model updating helped with this challenge to update the most sensitive model parameters separately using modal-based and time-domain model updating. • Sequential application of modal-based and time-domain model updating reduced the challenges due to ill-conditioning and modeling errors.

•
It was demonstrated that the updated FE model using only modal-based model updating was not capable of reflecting true response behavior of the structure in time domain.
• Concrete compressive strength was correlated with damage/deterioration in the monitored structure and could be used to assess the health condition of the structure. In this study, a 30% reduction in concrete compressive strength from its nominal value correctly showed significant deterioration in the studied girders.
It is noteworthy to mention that although the input load in this study was different from the traffic load during the operation of bridges, it provided a real-world exercise to validate the capability of the two-step model updating approach for damage identification of bridge structures.