Estimation of Soil–Structure Model Parameters for the Millikan Library Building Using a Sequential Bayesian Finite Element Model Updating Technique

We present a finite element model updating technique for soil–structure system identification of the Millikan Library building using the seismic data recorded during the 2002 Yorba Linda earthquake. A detailed finite element (FE) model of the Millikan Library building is developed in OpenSees and updated using a sequential Bayesian estimation approach for joint parameter and input identification. A two-step system identification approach is devised. First, the fixed-base structural model is updated to estimate the structural model parameters (including effective elastic modulus of structural components, distributed floor mass, and Rayleigh damping parameters) and some uncertain components of the foundation-level motion. Then, the identified structural model is used for soil–structure model updating wherein the Rayleigh damping parameters, the stiffness and viscosity of the soil subsystem (modeled using a substructure approach), and the foundation input motions (FIMs) are estimated. The identified model parameters are compared with state-of-practice recommendations. While a specific application is made for the Millikan Library, the present work offers a framework for integrating large-scale FE models with measurement data for model inversion. By utilizing this framework for different civil structures and earthquake records, key structural model parameters can be estimated from the real-world recorded data, which can subsequently be used for assessing and improving, as necessary, state-of-the-art seismic analysis and structural modeling techniques. This paper presents an effort towards using real-world measurements for large-scale FE model updating in the soil and structure, uniform soil time domain for joint parameter and input estimation, and thus paves the way for future applications in system identification, health monitoring, and diagnosis of civil structures.


Introduction
The dynamic response of a building structure to an earthquake excitation is occasionally the result of a complex interaction between the structural system and the underlying and surrounding geology. Since modeling the physics of a coupled soil-structure system in detail is an intricate undertaking, the state of practice has adopted simplified modeling and analysis procedures, such as the substructure approach (e.g., [1,2]). In this approach, the soil flexibility and energy dissipation are modeled using distributed spring and dashpot elements to which the foundation input motions (FIMs) are applied as uniform base excitations [3]. The stiffness and viscosity of these elements are derived using simplified analytical methods, which are nonetheless based on idealized and restrictive assumptionse.g., linear elastic response behavior of soil and structure, uniform soil half-space (or soil profiles with stiffness gradually varying with depth [4]), etc. These simplifying assumptions and the empirical nature of the mechanical analogs (e.g., soil springs and dashpots) could potentially lead to unquantified errors in predicting the seismic responses of real-world building structures, even though the simplified models have demonstrated acceptable accuracy in idealized case studies (e.g., [5]). Furthermore, integrating the substructure method with the Rayleigh damping model-which itself has a highly empirical nature-introduces another source of uncertainty in the seismic analysis of building structures. Clearly, there is an inconsistency between the relatively sophisticated mechanics-based techniques available for structural system modeling and the crude simplicity of the underlying assumptions that guide the damping and soil-structure interaction (SSI) models. It appears clear that the efforts on damping and SSI characterization should be guided by observations distilled from real-life data.
Soil-structure interaction (SSI) analysis has been an active research subject for more than 40 years (e.g., [6]). SSI effects can be classified into two categories: kinematic and inertial interaction effects [7]. Referred to as the FIMs, the earthquake excitations experienced by a structure-foundation system are affected by the stiffness contrast between the foundation stiffness and the surrounding geology. Therefore, FIMs are generally different from the Free-Field Motions (FFMs) that would have been recorded in the absence of the foundation system. The effects of foundation stiffness and geometry, which result in discrepancies between FIMs and FFMs, are collectively referred to as the kinematic interaction effects. Inertial interaction effects are a result of the foundation-superstructure mass, which imparts inertial forces onto the surrounding soil and causes the foundation to experience a response different from the FIMs. Due to the inertial interaction effects, the vibrating structure operates as a wave source and alters the wave field around the foundation system. In the substructure modeling approach, these two effects are treated separately (e.g., [8]). That is, FIMs are analytically or numerically calculated based on the soil and foundation properties (e.g., [9]). Then, inertial effects are represented by frequency-dependent impedance functions that are estimated through analytical, numerical, or experimental studies (e.g., [10]). Finally, the superstructure on top of the impedance functions is analyzed under the estimated FIMs.
In this study, we develop a model inversion framework that can be used to backcalculate the soil-structure model parameters from the seismic responses of real-world buildings. By repeating this effort for different buildings and earthquake records, the estimation results can be summarized and compared with the state-of-practice recommendations to improve the present seismic analysis capabilities and predictive models for building structures. We apply the model inversion framework that is devised in the present study to the data recorded from the Millikan Library building during the 2002 Yorba Linda earthquake, which was a low-amplitude seismic event. The response of the building during this small event was nearly linear elastic. As such, it provides an opportunity to compare the estimation results with linear elastic seismic analysis practices.
The utilized model inversion framework is based on a finite element (FE) model updating approach formulated in the time domain, using a sequential Bayesian estimation method [11]. The FE model updating approach has been developed and verified previously using numerically simulated data (e.g., [12,13]). This paper presents an effort towards using real-world measurements for large-scale FE model updating in the time domain, and thus paves the way for future applications. Although the presented model updating approach is used along with a linear FE model in this study, the approach can be readily extended to nonlinear FE models of civil structures.
In what follows, the Bayesian FE model updating approach and its algorithmic details are briefly discussed in Section 2. Section 3 provides the modeling details of the Millikan Library building. A two-step system identification approach is introduced in Section 4, which is applied to the Millikan Library data in Sections 5 and 6. An identifiability analysis accompanies each system identification step to select the estimation parameters. Finally, the key estimation results are summarized and compared with available literature and the state of practice in Section 7.

Data Assimilation through FE Model Updating in the Time Domain
The finite element model of a building structure depends on various "model parameters". Model parameters include, but are not limited to, inertial properties, damping parameters, parameters that characterize the stiffness properties of the structural and nonstructural components, and the time histories of the unmeasured components of base excitation. If the soil-structure interaction effects are included in the FE model, then the parameters used to model the soil subsystem and the time history of FIMs are also unmeasurable and unknown. Using the recorded acceleration response time histories of the building during an earthquake, our objective is to identify the best estimates of the unknown model parameters along with the input excitation time history, in order to minimize the discrepancy between the FE-predicted and measured structural responses.
The estimation problem is tackled by updating iteratively and sequentially the joint probability distribution function (PDF) of the unknown model parameters and the discrete values of the input motion time histories using a Bayesian inference method. The basics of the sequential FE model updating method for joint system and input identification are presented in [11]. Here, we briefly review the technical background of this method and highlight the new improvements introduced for this study.

Bayesian FE Model Updating in Time Domain for Joint System and Input Identification
The predicted response of a FE model at time step i to an earthquake excitation time history can be expressed as [14]ŷ i = h i θ, .. u g 1:i (1) whereŷ i ∈ R n y ×1 denotes the vector of structural response quantities predicted by the FE model at time step i (i = 1, 2, . . . , k), θ ∈ R n θ ×1 denotes the FE model parameter vector, .. The measured response vector of the structure, y i , is related to the FE-predicted response,ŷ i , as ( [12][13][14]) v i θ, where v i ∈ R n y ×1 is the simulation error vector, which accounts for the misfit between the measured and FE-predicted responses. This misfit stems from the output measurement noise, parameter uncertainties, and model errors. Model errors stem from simplifying idealizations adopted when devising the model as well as deviations in its geometry and configuration from the real-life structure it represents, which result in an inherent discrepancy between the FE model prediction and the measured structural responses [15]. Given an unbiased estimate of θ and .. u g 1:k , it is assumed here that the simulation error can be modeled as a stationary, zero-mean, and independent Gaussian white noise process [16]. Therefore, the likelihood function can be derived as where |R| denotes the determinant of the diagonal matrix R ∈ R n y ×n y , which is the timeinvariant covariance matrix of the simulation error vector, i.e., R = E v i v i T , ∀i. The unknown model parameter vector and time histories of the unknown input ground motion can be stacked into an extended unknown parameter vector defined as ψ = θ T , In a stochastic model inversion, the uncertainty in the unknown parameter vector is characterized by a joint probability distribution function (PDF). The objective in FE model updating is to find an estimate of the unknown parameters such that their joint posterior PDF given the measured output response of the structure is maximized. This is referred to as the maximum a posteriori (MAP) estimate. Using Bayes' rule, it can be observed that where p(y 1:k | ψ) = ∏ k i=1 p(v i ), and p(ψ) is the joint prior PDF of the extended unknown parameters.

Sequential Finite Element Model Updating Using Model Linearization
Solving the MAP estimation problem shown in Equation (4) using a batch optimization algorithm can be computationally demanding [11]. Therefore, a sequential estimation approach is introduced in [11] to improve the computational efficiency and convergence rate. In this approach, the estimation time domain is divided into successive overlapping time windows, referred to as the estimation windows. The estimation problem is solved iteratively at each estimation window to estimate the posterior mean vector and covariance matrix of the unknown parameters. At each iteration, the FE model is linearized to propagate the uncertainties in the prior PDF into the FE model response. The posterior mean vector and covariance matrix of the unknown parameter vector are then derived using Equation (4), and used as prior information for the next iteration. Once the iterations converge for an estimation window, the estimated mean vector and covariance matrix are transferred to the next estimation window and are used as prior information. In what follows, the mathematical process of the sequential Bayesian estimation approach is outlined.
Assume that the m th estimation window with length t l spans from time step t m 1 to time step t m 2 , where t l = t m 2 − t m 1 . Considering that p y t m , can be derived as [11] log p( ψ m | y t m where k 0 is a constant, andψ − m andP − ψ,m are the prior mean vector and the covariance matrix of the extended parameter vector at the m th estimation window.R ∈ R (t l ×n y )×(t l ×n y ) is a block diagonal matrix, in which the diagonals denote the simulation error covariance matrix R, and. :t m 2 ∂ψ m = 0, which results in a nonlinear algebraic equation in ψ m , as shown in [11]. This equation can be solved using an iterative first-order approximation of the FE response function atψ − m , which requires the computation of FE response sensitivities with respect to the extended unknown parameter vector, evaluated at the prior mean values. The FE response sensitivity matrix is denoted by C hereafter. As a result, the following (first-order approximate) equation for the MAP estimate of ψ m can be obtained [11] aŝ whereψ + m is the updated (or the posterior) mean estimate of ψ m . The term K =P ψy (P yy ) −1 is referred to as the Kalman gain matrix [17]. The matrixP ψy =P − ψ,m C T is the crosscovariance matrix of the unknown parameter vector and measured data, and P yy = CP − ψ,m C T +R is the covariance matrix of the measured data (the reader is referred to [11,18] for further derivation details). The subscript j is appended to the estimated mean vector and covariance matrix to denote the iteration number. The updatedψ + m,j from Equation (6) is iteratively used as the new point for the linearization of the nonlinear FE model to improve the estimation. This iterative prediction-correction procedure at each estimation window is equivalent to an iterative extended Kalman filtering (EKF) method for parameter-only estimation [17]. Following the EKF procedure, the prior covariance matrix P − ψ,m,j is updated to the posterior covariance matrixP + ψ,m,j at each prediction-correction iteration as inP+ Furthermore, to improve the convergence characteristics of the iterative predictioncorrection procedure, a disturbance matrix is added to the posterior covariance matrix at each iteration to provide the prior covariance matrix for the next iteration, i.e.,P − ψ,m,j+1 =P + ψ,m,j + Q, in which Q is a diagonal matrix with small positive diagonal entries (relative to the diagonal entries of the matrixP + ψ,m,j ). The matrix Q is referred to as the process noise covariance matrix in the Kalman filtering world [17].
Accurately computing the FE response sensitivities is an integral part of the presented sequential FE model updating. In previous studies (e.g., [11,12]), the response sensitivities have been computed using a direct differentiation method (DDM) (e.g., [19]), which requires extending the FE numerical scheme and implementing new codes in the FE software to calculate the response sensitivities. Although available in OpenSees [20], this implementation is not readily available as part of other FE simulation platforms. Therefore, to extend the applicability of the sequential FE model updating procedure to other FE simulation platforms, here we use the finite difference method (FDM) to compute the response sensitivities. The FDM method is implemented by perturbing the estimation parameters one at a time based on a central difference method with unequal spacing. We use parallel computing to enhance the computational efficiency of the FDM implementation. The Bayesian FE model updating algorithms are implemented in MATLAB [21], which calls several instances of OpenSees in parallel for FE response and response sensitivity computations at each prediction-correction iteration. Table 1 summarizes the sequential FE model updating algorithm for the joint estimation of model parameters and input motion time histories that is described above. It should be noted that the presented Bayesian FE model updating framework is applicable to both linear and nonlinear FE models. In case the FE model is linear, the model response function will still be a nonlinear function of the model parameters and input motion; however, the model updating process will become more computationally efficient.

1.
Determine the start and end points of the estimation windows-i.e., t m 1 , t m 2 , where m = 1, 2, . . . , N and N is the number of estimation windows.

2.
Set the initial mean vector and covariance matrix of the model parameter vector (θ 0 ,P θ 0 ) and input ground motion time history vector across the first estimation window (.  ). The initial mean values of the ground motion time history are usually set as zero, i.e.,.
. u g 0,t 1 1 :t 1 Use the inverse of the initial values of the model parameters as the scaling factors, i.e., a n = 1 θ 0,n . Set upψ

Obtain the FE responses usingψ
, then move to the next estimation window (m = m + 1, go to step 4); otherwise, iterate again at the current estimation window (j = j + 1, go to step 4.3).

Correction for Constraints
The estimation process described earlier is not constrained and may result in the estimation of non-physical values for the model parameters. This issue is resolved by correcting the posterior estimation results for a set of predefined lower-and upper-bound constraints. Once any of the posterior mean values exceed their designated lower or upper limit, the posterior Gaussian PDF (characterized by the updated mean vector and covariance matrix) is truncated at the constraint edges. The mean vector and the covariance matrix of the truncated PDF are then calculated. The constraint correction approach is borrowed from [17,22]. The algorithmic details are provided in [23] and are not repeated here for brevity.

Adaptive Scaling of the Unknown Model Parameters
In order to improve the performance of the estimation algorithm, the FE model parameters are scaled adaptively at each iteration. The scaling allows the FE model responses to have relatively similar sensitivities with respect to different estimation parameters at each estimation point. The scaled FE model parameter vector is defined as θ s = Aθ, in which A is an n θ × n θ diagonal scaling matrix. It should be noted that only the FE model parameters are scaled, and the vector of unknown input ground motion remains unscaled. The parameter scaling factors (i.e., the diagonal entries of matrix A denoted by a n , n = 1, 2, . . . , n θ ) are calculated to result in equal corresponding diagonal entries in the Fisher Information Matrix of the scaled parameters, which is approximated as I s = (C s θ ) T (C s θ ) [14], wherein C s θ denotes the FE response sensitivity matrix with respect to the scaled FE model parameter vector. The n th diagonal entry of I s is I s n = 1 a n 2 ∂h ∂θ n T ∂h ∂θ n . The scaling factor a n is calculated so that I s n is equal to the mean of diagonal entries of the Fisher Information Matrix corresponding to the vector of unknown input ground motion, i.e., for the m th estimation window:

The Millikan Library Building
The Millikan Library is a reinforced concrete shear wall building, with a basement and nine stories above the ground, located on the California Institute of Technology (Caltech) campus ( Figure 1). The Millikan Library has been the subject of many studies, especially in the fields of soil-structure system identification, due to its unique structural and soil properties, and its relatively old and dense instrumentation, which recorded its seismic responses during numerous earthquakes [24][25][26][27][28][29][30]. The Millikan Library building is 43.9 m tall above ground, including the roof level. It has a 4.3 m deep basement below the ground level. The basement is encased by surrounding shear walls. Precast concrete claddings are installed on the north and south faces of the building from the second floor up to the roof. The lateral force-resisting system comprises shear walls in the north-south direction on the east and west sides of the building and core walls around the elevator shaft in the north-south and east-west directions. The floor system consists of concrete slabs supported by reinforced concrete beams. Lightweight aggregate concrete is used for slabs and beams at all floors, while regular aggregate concrete is used for foundations, columns, and walls. The foundation system consists of a 1.2 m deep central pad, 6.0 m below the ground level, and two foundation strips on the north and south sides of the building 5.0 m below the ground level. Four stepped beams at the four corners of the foundation connect the central foundation pad to the north and south foundation strips.
The building is instrumented using a total of 36 uniaxial accelerometers, measuring six acceleration responses at the foundation (or basement floor) level, and three translational acceleration responses on each floor. The instrumentation details are obtained from [31] and are provided in Figure A1 and Table A1 in Appendix A. Multiple earthquake records are available for the Millikan Library; a relatively recent survey of earthquake data recorded at the building can be found in [28]. For the purpose of this study, we use the 2002 Yorba Linda earthquake record, which is a low-amplitude earthquake (PGA of~0.6% g). Two notes should be made about this earthquake record of the Millikan Library. First, the location of sensors on the foundation (or basement floor) could not be determined precisely. Second, recordings of the east sensor on the second and eighth floors (measuring NS direction) were not available during the 2002 Yorba Linda Earthquake. Therefore, only 34 measurement channels are available for this record. Since the level of considered earthquakes is low, the effects of material and geometric nonlinearities in the dynamic response of the structure are expected to be minor, and therefore, the use of a linear model is justified.
gramming Interface (API). For the purpose of this study, we prepared a linear FE model of the Millikan Library building. The model comprises linear elastic beam-column elements for the buildings' beams and columns, and quadrilateral shell (DKGQ) elements with linear elastic (elastic membrane plate) sections are used for shear walls and slabs. The kinematic interaction of precast claddings (installed on the north and south faces of the building) with the structural system is modeled using diagonal brace elements, which are included in the north and south frame bays from the second story to the roof. The brace elements are assumed to be linear elastic with square cross-sections of 0.1 m × 0.1 m. By adjusting the elastic modulus of the diagonal braces, the stiffness and force contributions of the precast panels to the lateral response of the structure can be characterized. The kinematic interaction between the structure and other nonstructural components and systems (e.g., partition walls and stairs) are not modeled explicitly; they are expected to be lumped into the stiffness properties of the structural members and precast claddings. The inertial effects of nonstructural components (i.e., stairs, precast panels, and elevators) are accounted for by modeling their corresponding mass and weight contributions. The weight of each precast panel is estimated to be about 50 kN, which conforms with [24]. A uniformly distributed load is applied on the floor slabs to account for live loads and weight of partition walls, later treated as an unknown parameter to be estimated. A uniformly distributed load of 2.5 kN/m 2 is applied on the foundation to account for the live load of library shelves in the basement. Moreover, the mass of soil filling between the foundation surface and basement floor is accounted for by applying the corresponding mass on the central foundation pad and the two foundation strips. The mass of structural components is calculated based on the nominal concrete densities specified in the structural drawings: − = 1760  Figure 2 shows the geometry of the model. Different colors in this figure present different material property groups used to define the shell (shear wall and slab) and beam-column elements. More details on material parameters

FE Model Development
Using the available structural drawings, a detailed FE model of the structural system was developed. We used the graphical user interface of SAP2000 software [32] to develop the initial geometry of the model. The SAP2000 model was then transferred to OpenSees [20] using a custom-developed MATLAB script based on the SAP2000 Application Programming Interface (API). For the purpose of this study, we prepared a linear FE model of the Millikan Library building. The model comprises linear elastic beam-column elements for the buildings' beams and columns, and quadrilateral shell (DKGQ) elements with linear elastic (elastic membrane plate) sections are used for shear walls and slabs. The kinematic interaction of precast claddings (installed on the north and south faces of the building) with the structural system is modeled using diagonal brace elements, which are included in the north and south frame bays from the second story to the roof. The brace elements are assumed to be linear elastic with square cross-sections of 0.1 m × 0.1 m. By adjusting the elastic modulus of the diagonal braces, the stiffness and force contributions of the precast panels to the lateral response of the structure can be characterized. The kinematic interaction between the structure and other nonstructural components and systems (e.g., partition walls and stairs) are not modeled explicitly; they are expected to be lumped into the stiffness properties of the structural members and precast claddings.
The inertial effects of nonstructural components (i.e., stairs, precast panels, and elevators) are accounted for by modeling their corresponding mass and weight contributions. The weight of each precast panel is estimated to be about 50 kN, which conforms with [24]. A uniformly distributed load is applied on the floor slabs to account for live loads and weight of partition walls, later treated as an unknown parameter to be estimated. A uniformly distributed load of 2.5 kN/m 2 is applied on the foundation to account for the live load of library shelves in the basement. Moreover, the mass of soil filling between the foundation surface and basement floor is accounted for by applying the corresponding mass on the central foundation pad and the two foundation strips. The mass of structural components is calculated based on the nominal concrete densities specified in the structural drawings: ρ c−LW = 1760 kg m 3 for lightweight aggregate concrete material, ρ RC−LW = 1900 kg m 3 for lightweight aggregate reinforced concrete, and ρ RC = 2500 kg m 3 for regular aggregate reinforced concrete. Figure 2 shows the geometry of the model. Different colors in this Buildings 2023, 13, 28 9 of 29 figure present different material property groups used to define the shell (shear wall and slab) and beam-column elements. More details on material parameters are provided in Section 5. A sensitivity analysis is performed to ensure convergence of analysis results with respect to the mesh size. The model in the presented configuration consists of 1885 frame elements and 4043 shell elements, corresponding to a total of 27,526 degrees of freedom.
Buildings 2022, 12, x FOR PEER REVIEW 9 of 28 are provided in Section 5. A sensitivity analysis is performed to ensure convergence of analysis results with respect to the mesh size. The model in the presented configuration consists of 1885 frame elements and 4043 shell elements, corresponding to a total of 27,526 degrees of freedom. A gravity analysis was performed before dynamic time history analyses by applying constant g acceleration in the vertical direction. The damping energy dissipation for the time history analysis is defined using mass-and stiffness-proportional Rayleigh damping. The Rayleigh damping parameters are treated as unknowns to be estimated. For the purpose of time history analyses, Newmark's constant average acceleration method is used for time-stepping with increments of = 0.03 s. The time step size is selected following a convergence study to ensure a large time step size that provides adequate accuracy. All of the measured acceleration response time histories are also resampled at = 0.03 s for the model updating phase.

A Two-Step System Identification Approach
A two-step FE model updating approach is used for model and input identification of the soil-structure system. In the first step, the measured acceleration responses at the foundation level are used to calculate six components of the foundation-level motions. The foundation-level motions are then used as uniform base input excitations, and the measured responses of the structure are used as output measurements to update the FE model of the "fixed-base" structure. The objective is to estimate the model parameters characterizing the structural model, regardless of the soil subsystem. During this process, we discovered inconsistencies in the recordings obtained from two of the foundation sensors, as is discussed further in the next section. These sensors are installed inside a utility tunnel-referred to as the steam tunnel-next to the foundation slab. Therefore, the acceleration response time histories of these two sensors are assumed as unknown input motions and are estimated jointly with the structural model parameters. This step results in FE model updating with partially unknown inputs.
The second step is an output-only FE model updating procedure, wherein the identified structural model parameters are fixed at their mean estimates obtained from the first step, and the three translational components of FIM and other parameters characterizing the soil-structure model are estimated jointly. Since the Millikan Library foundation is not deep and does not have large dimensions (comparable with the lengths of incident waves), the rotational components of FIM are not considered in this study.
To decide about and select the estimation parameters, an identifiability analysis is performed, as is outlined in the next section for each identification step. Before proceeding with the real-life data, we verified the two-step FE model updating procedure using A gravity analysis was performed before dynamic time history analyses by applying constant g acceleration in the vertical direction. The damping energy dissipation for the time history analysis is defined using mass-and stiffness-proportional Rayleigh damping. The Rayleigh damping parameters are treated as unknowns to be estimated. For the purpose of time history analyses, Newmark's constant average acceleration method is used for time-stepping with increments of ∆t = 0.03 s. The time step size is selected following a convergence study to ensure a large time step size that provides adequate accuracy. All of the measured acceleration response time histories are also resampled at ∆t = 0.03 s for the model updating phase.

A Two-Step System Identification Approach
A two-step FE model updating approach is used for model and input identification of the soil-structure system. In the first step, the measured acceleration responses at the foundation level are used to calculate six components of the foundation-level motions. The foundation-level motions are then used as uniform base input excitations, and the measured responses of the structure are used as output measurements to update the FE model of the "fixed-base" structure. The objective is to estimate the model parameters characterizing the structural model, regardless of the soil subsystem. During this process, we discovered inconsistencies in the recordings obtained from two of the foundation sensors, as is discussed further in the next section. These sensors are installed inside a utility tunnel-referred to as the steam tunnel-next to the foundation slab. Therefore, the acceleration response time histories of these two sensors are assumed as unknown input motions and are estimated jointly with the structural model parameters. This step results in FE model updating with partially unknown inputs.
The second step is an output-only FE model updating procedure, wherein the identified structural model parameters are fixed at their mean estimates obtained from the first step, and the three translational components of FIM and other parameters characterizing the soil-structure model are estimated jointly. Since the Millikan Library foundation is not deep and does not have large dimensions (comparable with the lengths of incident waves), the rotational components of FIM are not considered in this study.
To decide about and select the estimation parameters, an identifiability analysis is performed, as is outlined in the next section for each identification step. Before proceeding with the real-life data, we verified the two-step FE model updating procedure using numerically simulated data to examine the effectiveness of the proposed estimation algorithm and the accuracy of its results. The verification studies are documented in [33] and are not included here for brevity.

Model Identifiability and Parameter Selection
The identifiability of model parameters is investigated using an information-theoretic approach based on [34]. In this approach, a set of candidate parameters are selected first, and then, the entropy gain of each parameter is calculated. The entropy gain is used as a quantitative metric to measure the information each model parameter receives from the measurement data. The entropy gain is compared between the estimation parameters to assess their relative identifiability. Moreover, the mutual entropy gain (or mutual gain) between parameter pairs is used to investigate the dependence between each pair of model parameters. The entropy gain and mutual gain are then used to determine the most identifiable parameter set. This is an alternative procedure to prior methods that, for example, adaptively group updating parameters based on the sensitivities of a model updating objective function with respect to the parameters (e.g., see [35,36]).
Fifteen (15) structural model parameters characterizing the material, inertial, and damping properties were initially selected as the estimation parameters for identifiability assessment at this step. These parameters were selected initially based on our judgment on their potential contribution in the structural response. For a given input motion, the entropy gain for each parameter is a function of the parameter values, which are unknown in advance. Without the knowledge of correct parameter values, their nominal values are used for prior identifiability assessment, following the recommendation in [34]. These parameters and their nominal values are listed in Table 2. The logic based on which the nominal values are estimated is further described below.
The material parameters consist of the effective elastic modulus of the floor system (beams and slabs) and the vertical/lateral system (walls and columns). The material parameters are grouped at the story level and along the height of the building. With the exception of the first three stories, the effective elastic modulus of structural components is grouped together at the upper stories. The first three stories are considered separately because previous studies have suggested minor structural damage in the first stories of the building after the 1971 San Fernando earthquake (e.g., [26,37,38]). The effective elastic modulus for floor slabs and beams is calculated by applying a 35% reduction factor on the nominal elastic modulus of lightweight aggregate concrete, which is derived based on ACI 318-14 [39] as E c = 0.043(ρ c−LW ) 1.5 f c , where f c = 27.5 MPa is the nominal compressive strength. The effective elastic modulus for columns and walls is calculated by applying a 70% reduction factor on the nominal elastic modulus calculated as E c = 4700 f c , where f c = 27.5 MPa. For foundation, a 35% reduction factor is applied on the elastic modulus, which is calculated considering f c = 20.7 MPa, to account for flexural cracking.
To estimate the elastic modulus of diagonal brace elements representing the precast claddings, the results concluded from modal analysis performed during the construction period are used. According to [6,26], the installation of precast panels resulted in a 25% increase in the lateral stiffness of the building structure in the EW direction. The initial values for the elastic modulus of the diagonal brace elements in the model are tuned iteratively and manually to yield a similar increase in the stiffness of the first mode in the EW direction. Consequently, the estimated elastic modulus of diagonal braces is found to be 20 GPa. The initial value of the distributed floor mass (to account for live load and mass of nonstructural elements) is estimated to be about 250 kg/m 2 . Estimating the accurate weight and mass contribution of mechanical equipment on the roof was not possible; therefore, the equivalent distributed roof mass is approximately estimated as 300 kg/m 2 . Finally, the initial values of the Rayleigh damping parameters are estimated by assuming a 5% damping ratio for the first and second modes. The identifiability of model parameters is assessed using the measured foundationlevel motion, including all translational and rotational components. As mentioned earlier, the recordings of two acceleration channels on the foundation are later treated as unknowns and estimated through the model updating process. Nevertheless, these measured recordings are used for the identifiability analysis as approximate inputs. This approximation is expected to have negligible effects on the overall identifiability results. As described in [34], evaluating the entropy gain and mutual gain also requires the prior covariance matrix of the unknown parameter vector (or the prior variance of parameters). A 10% coefficient of variation is assumed for all parameters to derive the a priori covariance matrix. Figure 3 displays the relative entropy gain of the fifteen candidate estimation parameters. The presented entropy gains are relative, which means that the entropy gain values are scaled with respect to the largest value. This plot presents the relative information that the model responses (which represent the measurements) carry about each model parameter. Those parameters that receive small information from the measurement data (e.g., parameters #2 to #6) are likely to be unidentifiable. Furthermore, Figure 4 shows the relative mutual entropy gain between the parameter pairs. This figure is used to investigate the dependence between each two parameters. The diagonal components in this figure represent the entropy gain of each parameter. The dark off-diagonal colors in this figure indicate strong relative dependence between parameter pairs. For example, Figure 4 suggests that there is a mutual dependence between parameters #7 and #8, and some mutual dependence between parameters #7, #8, and #10. Moreover, there are some competing effects between these parameters and parameter #14. These observations are physically meaningful since parameters #7, #8, and #10 characterize the stiffness of the building, and parameter #14 contributes to the building mass. Figure 3 along with Figure 4 can be used to guide the selection of the estimation parameters. The parameters that gain relatively small amounts of information from the measurement data or have strong dependencies on other estimation parameters can be put aside (i.e., fixed) or merged with other parameters, if physically meaningful, to enhance the identifiability. Indeed, the process is not definitive and requires some engineering judgment in selecting the final estimation parameter set.  Table 2 for parameter IDs).  Table 2 for parameter IDs).
The strong mutual gain between parameters #7 and #8 suggests merging these two parameters into a single unknown parameter. This is justifiable as these parameters characterize the effective elastic modulus of columns and walls at the basement level and the first story, which are expected to be close. The mutual gains between other stiffness-related parameters are not significant; nevertheless, their small entropy gains in Figure 3 suggest merging them together. Furthermore, since the initial estimate of mass-related parameters (#14 and #15) are close and their individual entropy gains are small, they are also merged together to improve their identifiability. Finally, while the entropy gains of parameters #1 and #12 are small, we still consider them as unknown parameters to be estimated; however, the resulting estimates are expected to be inaccurate and include large uncertainties. As a result of the outlined identifiability assessment process, we end up with six estimation parameters, as listed in Table 3.   Table 2 for parameter IDs).
Buildings 2022, 12, x FOR PEER REVIEW 12 of 28 Figure 3. Relative entropy gain of the fifteen model parameters (see Table 2 for parameter IDs).  Table 2 for parameter IDs).
The strong mutual gain between parameters #7 and #8 suggests merging these two parameters into a single unknown parameter. This is justifiable as these parameters characterize the effective elastic modulus of columns and walls at the basement level and the first story, which are expected to be close. The mutual gains between other stiffness-related parameters are not significant; nevertheless, their small entropy gains in Figure 3 suggest merging them together. Furthermore, since the initial estimate of mass-related parameters (#14 and #15) are close and their individual entropy gains are small, they are also merged together to improve their identifiability. Finally, while the entropy gains of parameters #1 and #12 are small, we still consider them as unknown parameters to be estimated; however, the resulting estimates are expected to be inaccurate and include large uncertainties. As a result of the outlined identifiability assessment process, we end up with six estimation parameters, as listed in Table 3. Table 3. Final selected model parameters for Step 1 system identification.  Table 2 for parameter IDs).

Parameter ID Description Value
The strong mutual gain between parameters #7 and #8 suggests merging these two parameters into a single unknown parameter. This is justifiable as these parameters characterize the effective elastic modulus of columns and walls at the basement level and the first story, which are expected to be close. The mutual gains between other stiffness-related parameters are not significant; nevertheless, their small entropy gains in Figure 3 suggest merging them together. Furthermore, since the initial estimate of mass-related parameters (#14 and #15) are close and their individual entropy gains are small, they are also merged together to improve their identifiability. Finally, while the entropy gains of parameters #1 and #12 are small, we still consider them as unknown parameters to be estimated; however, the resulting estimates are expected to be inaccurate and include large uncertainties. As a result of the outlined identifiability assessment process, we end up with six estimation parameters, as listed in Table 3.

Joint System and Partial Input Identification Using Yorba Linda Earthquake Data
The seismic data recorded at the Millikan Library building during the 2002 Yorba Linda earthquake are utilized for the FE model updating. Here, we present two model updating cases. In Case 1, the recordings of the six accelerometers installed on the foundation are averaged to find the time histories of translational and rotational components of the foundation-level motion. These input excitations along with the measured structural responses obtained from all the channels are used for an "input-output" FE model updating to estimate the six model parameters listed in Table 3. In Case 2, we assume that the recordings of the two sensors located in the steam tunnel, i.e., Channels #4 and #5, are unknown. This resulted in FE model updating with partially unknown inputs to jointly estimate the six model parameters and the acceleration response time history at two sensor locations. Figure 5 presents the relative root mean square error (RRMSE) between the measured and FE predicted acceleration responses for both Case 1 and Case 2. The figure also shows the RRMSE values obtained from the initial FE model (i.e., based on the initial parameter values). As can be seen in this figure, the RRMSE is consistently reduced from the initial to the updated FE model in both cases. However, Case 2 shows more reduction in RRMSE for all measurement channels compared to Case 1 except for Channels #4 and #5, in which RRMSEs are increased. This suggests that these two sensor measurements do not comply with the other measurements on the foundation level. Hence, the updated FE model in Case 2 is most likely a better representation of the real-world building. As shown in Appendix A, these two channels correspond to the sensors installed in the steam tunnel-Channel #4 measures the acceleration response in the NS-direction and Channel #5 in the Up-Down or Z-direction. The steam tunnel is a rectangular reinforced concrete utility access tunnel built next to the foundation level and is not part of the basement structure. After the 1971 San Fernando Earthquake, a 0.5 mm crack around the steam tunnel, at the point where the tunnel connects to the basement wall, was reported [37]. Based on these observations, it is concluded that the tunnel is most likely not fully connected to the structural system, and thus, the readings of the accelerometers located in there may not correctly represent the foundation level motions. Figure 6 compares the measured and estimated acceleration response time histories at Channels #4 and #5. The differences between the measured and estimated acceleration time histories are minor and more prominent for Channel #5, which is the Z-direction measurement.
The initial and final estimates of the six model parameters, along with their final estimated coefficients of variation (COVs), are listed in Table 4. The coefficient of variation quantifies the estimation uncertainty. In an unbiased estimation, the COV should merge to zero. The non-zero values for the COV mean that the estimation might be uncertain. The larger the COV, the more uncertainty there is in the estimated model parameter values. As can be seen in Table 4, the COV of the two parameters with the highest COV in case 1, i.e., E Clad and m, is significantly reduced in Case 2. It is also seen that the estimate of these parameters is unrealistically low in Case 1, while their estimates in Case 2 seem to be more realistic. These results, along with the RRMSE results shown in Figure 5, indicate that Case 2 is likely to be more reliable than Case 1, and hence, the parameter estimates from this case will be used for the next step.
pendix A, these two channels correspond to the sensors installed in the steam tunnel-Channel #4 measures the acceleration response in the NS-direction and Channel #5 in the Up-Down or Z-direction. The steam tunnel is a rectangular reinforced concrete utility access tunnel built next to the foundation level and is not part of the basement structure. After the 1971 San Fernando Earthquake, a 0.5 mm crack around the steam tunnel, at the point where the tunnel connects to the basement wall, was reported [37]. Based on these observations, it is concluded that the tunnel is most likely not fully connected to the structural system, and thus, the readings of the accelerometers located in there may not correctly represent the foundation level motions. Figure 6 compares the measured and estimated acceleration response time histories at Channels #4 and #5. The differences between the measured and estimated acceleration time histories are minor and more prominent for Channel #5, which is the Z-direction measurement.  The initial and final estimates of the six model parameters, along with their final estimated coefficients of variation (COVs), are listed in Table 4. The coefficient of variation quantifies the estimation uncertainty. In an unbiased estimation, the COV should merge to zero. The non-zero values for the COV mean that the estimation might be uncertain. The larger the COV, the more uncertainty there is in the estimated model parameter values. As can be seen in Table 4, the COV of the two parameters with the highest COV in Figure 6. Comparison of the measured and estimated acceleration response time histories at Channels #4 and #5. The right-hand-side plots magnify the response time history between 1 and 7 s. Each figure's title shows the sensor location and recording direction, e.g., "Basement, North West Corner-Z" means acceleration response time history in the Z-direction, recorded by the sensor located on the north-west corner of the basement.  Figure 7 shows how well the updated model predictions match the measurement records for both model updating cases, wherein selected measured acceleration response time histories are shown together with those estimated using the initial and final updated FE models. The figure shows that while the initial FE model response predictions have remarkable discrepancies with the measured responses, the updated model predictions have good agreement with the measurements. It also shows that the updated model predictions in Case 2 have better agreement with the measurements than Case 1, which further proves the reliability of Case 2 over Case 1. In the next paragraph, we describe the details of the estimation algorithm setup. The initial values of model parameters (̂0) are listed in Table 3, and their initial coefficient of variation is selected as 10%. The estimation constraints for model parameters are selected as 0.1̂0 ≤ ≤ 5̂0. The initial estimate of the discrete values of the unknown ground motion time history is selected to be zero, and their initial standard deviation is selected as 0.05 Rad s 2 . A variable estimation window configuration (length and overlap) is considered. During the first 6 s of the earthquake (which corresponds to the strong motion part), a small window length (60 time steps = 1.8 sec) with 50% overlap is used. As we The initial values of model parameters (θ 0 ) are listed in Table 3, and their initial coefficient of variation is selected as 10%. The estimation constraints for model parameters are selected as 0.1θ 0 ≤ θ ≤ 5θ 0 . The initial estimate of the discrete values of the unknown ground motion time history is selected to be zero, and their initial standard deviation is selected as 0.05 Rad s 2 . A variable estimation window configuration (length and overlap) is considered. During the first 6 s of the earthquake (which corresponds to the strong motion part), a small window length (60 time steps = 1.8 s) with 50% overlap is used. As we march through the time history, the estimation window length is gradually increased to 100 time steps (3 s) with 20% overlap. Selecting smaller lengths and larger overlaps for the estimation windows at the beginning of the earthquake ensures the incremental absorption of the information when the model has large uncertainties and the measurements have large information content (due to strong motion). The window length is subsequently increased and the overlap is reduced to improve the computational efficiency. The estimation window configuration can affect the computation demand and the estimation accuracy and should be selected based on experience and adjusted by observing the performance of the estimation algorithm.
The process noise covariance matrix Q (see Table 1) is selected as a diagonal matrix.
The diagonal entry corresponding to the i th model parameter is selected as qθ i 2 , where q = 0.001 andθ i is the mean estimate of the i th model parameter. This means that the RMS of the process noise corresponding to the model parameters is taken as 0.1% of the mean estimate of the associated model parameter. Thus, this part of the matrix Q would be time-varying. The diagonal entries corresponding to the unknown ground motion time histories are constant, and are taken as (0.001g) 2 . This means that the root mean square (RMS) of the process noise corresponding to the input excitation is time-invariant and is equal to 0.1% g. This selection is performed based on our previous experience with the estimation algorithm to ensure stability and proper convergence rate (more discussions are provided in [11]). Finally, the simulation error covariance matrix R is also selected as a time-invariant diagonal matrix, whose diagonal entries are selected as (0.001g) 2 . In other words, the simulation error is modeled as a 0.1% g RMS Gaussian white noise process.

Step 2 System Identification: FE Model Updating of the Soil-Structure System with Unknown Inputs
At this stage, we use the updated model of the superstructure (Case 2) for soilstructure system identification. The objective is to estimate the FIMs and the parameter characterizing the soil-structure system, including the stiffness and viscosity of the soil springs and dashpots used for modeling the inertial soil-structure interaction effects. For this purpose, and similar to the previous step, we start with an identifiability analysis by evaluating the information contained in the model responses. Then, we use the real measurement data to estimate the soil-structure model parameters.

Model Identifiability and Parameter Selection
Distributed linear soil springs and dashpots are included underneath the foundation slab of the updated FE model, obtained from the previous step. Three linear springs and three linear viscous dashpots are modeled independently in the x-, y-, and z-direction (corresponding to EW, NS, and Up-Down directions, respectively) at each nodal point of the foundation slab. The stiffness of soil springs and viscosity of dashpots are computed using the subgrade modulus (i.e., soil stiffness per unit area) and viscosity modulus (i.e., soil viscosity per unit area), respectively. The spring stiffness and dashpot viscosity at each nodal point are calculated by multiplying the tributary area of the nodal point by the corresponding subgrade modulus and viscosity modulus, respectively.
The Millikan Library has a two-level foundation system consisting of a central pad and two north and south foundation strips, as shown in Figure 8. Six unknown subgrade moduli-namely k x , k y 1, k y 2, k z 1, k z 2, and k z 3-are defined for different foundation regions, where, for example, the parameter k z 2 characterizes the vertical subgrade modulus for the interior area of the central pad. Two strips along the east and west edges of the central pad, which are underneath the two box shear walls at the east and west sides of the building, are assigned a different vertical subgrade modulus (k z 3) based on the recommendation provided in [3] to account for the rotational stiffness of the soil subsystem. Similarly, a different vertical subgrade modulus (k z 1) is assigned to the two north and south foundation strips. Likewise, six (unknown) parameters-namely c x , c y 1, c y 2, c z 1, c z 2, and c z 3-are used to define the viscosity modulus. Different foundation regions and their corresponding subgrade and viscosity modulus are shown in Figure 9. The twelve unknown subgrade and viscosity moduli along with the Rayleigh damping parameters (α and β) and the effective elastic modulus of foundation slabs (E Found ) are selected as parameters for the identifiability analysis. Although estimated in the previous step, the Rayleigh damping parameters are included for the soil-structure system identification since the damping model behavior of the fixed-based structure is expected to be different from the flexible-base structure. To assess the identifiability of these fifteen (15) unknown parameters, we pursue the same identifiability analysis approach presented before. Initial estimates of soil model parameters are assigned using the NIST standard guidelines [3]. The FIMs are unknown at this stage and will be estimated through the model updating process. In the absence of the correct FIMs, the three translational components of the foundation-level motion, estimated from the previous step, are used as the FIMs for identifiability assessment. The parameters used for identifiability analysis and their corresponding initial values are listed in Table 5. The initial values of the Rayleigh damping parameters are selected based on the estimation results obtained in the previous step. Moreover, a 20% coefficient of variation is assumed for all parameters to derive the a priori covariance matrix. Figure 10 displays the relative entropy gain of the fifteen estimation parameter candidates. Figure 11 displays the mutual entropy gain between parameter pairs. In Figure  11, the diagonal terms (i.e., entropy gain of estimation parameters) are found to be much larger than the off-diagonal terms (i.e., mutual entropy gain between parameter pairs); hence, including both of them in the figure would underrepresent the off-diagonal terms since they would appear with a hard-to-see, light color. Therefore, the diagonal terms are excluded from this figure. These two figures can be used together for selecting the estimation parameters. Figure 10 shows that some of the parameters (e.g., parameters #2, #8, #11, and #12) receive little information from the measurements. These parameters are, therefore, unlikely to be identifiable. Based on this figure, parameters #2 and #3, which repre- To assess the identifiability of these fifteen (15) unknown parameters, we pursue the same identifiability analysis approach presented before. Initial estimates of soil model parameters are assigned using the NIST standard guidelines [3]. The FIMs are unknown at this stage and will be estimated through the model updating process. In the absence of the correct FIMs, the three translational components of the foundation-level motion, estimated from the previous step, are used as the FIMs for identifiability assessment. The parameters used for identifiability analysis and their corresponding initial values are listed in Table 5. The initial values of the Rayleigh damping parameters are selected based on the estimation results obtained in the previous step. Moreover, a 20% coefficient of variation is assumed for all parameters to derive the a priori covariance matrix. Figure 10 displays the relative entropy gain of the fifteen estimation parameter candidates. Figure 11 displays the mutual entropy gain between parameter pairs. In Figure  11, the diagonal terms (i.e., entropy gain of estimation parameters) are found to be much larger than the off-diagonal terms (i.e., mutual entropy gain between parameter pairs); hence, including both of them in the figure would underrepresent the off-diagonal terms since they would appear with a hard-to-see, light color. Therefore, the diagonal terms are excluded from this figure. These two figures can be used together for selecting the estimation parameters. Figure 10 shows that some of the parameters (e.g., parameters #2, #8, #11, and #12) receive little information from the measurements. These parameters are, there- To assess the identifiability of these fifteen (15) unknown parameters, we pursue the same identifiability analysis approach presented before. Initial estimates of soil model parameters are assigned using the NIST standard guidelines [3]. The FIMs are unknown at this stage and will be estimated through the model updating process. In the absence of the correct FIMs, the three translational components of the foundation-level motion, estimated from the previous step, are used as the FIMs for identifiability assessment. The parameters used for identifiability analysis and their corresponding initial values are listed in Table 5.
The initial values of the Rayleigh damping parameters are selected based on the estimation results obtained in the previous step. Moreover, a 20% coefficient of variation is assumed for all parameters to derive the a priori covariance matrix. Figure 10 displays the relative entropy gain of the fifteen estimation parameter candidates. Figure 11 displays the mutual entropy gain between parameter pairs. In Figure 11, the diagonal terms (i.e., entropy gain of estimation parameters) are found to be much larger than the off-diagonal terms (i.e., mutual entropy gain between parameter pairs); hence, including both of them in the figure would underrepresent the off-diagonal terms since they would appear with a hard-to-see, light color. Therefore, the diagonal terms are excluded from this figure. These two figures can be used together for selecting the estimation parameters. Figure 10 shows that some of the parameters (e.g., parameters #2, #8, #11, and #12) receive little information from the measurements. These parameters are, therefore, unlikely to be identifiable. Based on this figure, parameters #2 and #3, which represent the subgrade modulus in the north direction for the foundation beams and central pad, respectively, are merged together. This is due to the small relative entropy gain of parameter #2. Similarly, parameters #8 and #9, and #10 to #12 are merged together to improve their identifiability. Moreover, Figure 11 shows mutual dependence between parameter #15 and parameters #4 to #6. This implies that the estimation of the effective elastic modulus of foundation slabs will likely affect the estimation accuracy of the vertical and rocking soil stiffnesses. Hence, parameter #15 is fixed at its nominal value given the relatively large thickness of foundation slab, which is expected to result in linear elastic behavior. Based on this identifiability study, the estimation parameters characterizing the soil-structure model are reduced to ten, which are k x , k y , k z 1, k z 2, k z 3, c x , c y , c z , α, and β. relatively large thickness of foundation slab, which is expected to result in linear elastic behavior. Based on this identifiability study, the estimation parameters characterizing the soil-structure model are reduced to ten, which are kx, ky, kz1, kz2, kz3, cx, cy, cz, α, and β.  Figure 10. Relative entropy gain of the fifteen model parameters (see Table 5 for parameter IDs). Figure 11. Relative mutual entropy gains between the parameter pairs (see Table 5 for parameter IDs).  Table 5 for parameter IDs). Figure 10. Relative entropy gain of the fifteen model parameters (see Table 5 for parameter IDs). Figure 11. Relative mutual entropy gains between the parameter pairs (see Table 5 for parameter IDs).

Joint System and Input Identification Using the Recorded Yorba Linda Earthquake Data
The seismic data recorded at the Millikan Library building during the 2002 Yorba Linda earthquake are now utilized in an output-only FE model updating procedure to estimate the eight soil parameters, two Rayleigh damping parameters, and the three components of the FIM in EW, NS, and Up-Down directions (no rotational component is considered for the FIM). The model parameters characterizing the superstructure (except the Rayleigh damping parameters) are kept constant based on the estimated results obtained from the previous system identification step. Except the initial coefficient of variation of Figure 11. Relative mutual entropy gains between the parameter pairs (see Table 5 for parameter IDs).

Joint System and Input Identification Using the Recorded Yorba Linda Earthquake Data
The seismic data recorded at the Millikan Library building during the 2002 Yorba Linda earthquake are now utilized in an output-only FE model updating procedure to estimate the eight soil parameters, two Rayleigh damping parameters, and the three components of the FIM in EW, NS, and Up-Down directions (no rotational component is considered for the FIM). The model parameters characterizing the superstructure (except the Rayleigh damping parameters) are kept constant based on the estimated results obtained from the previous system identification step. Except the initial coefficient of variation of estimation parameters, which is selected as 20%, and the estimation constraints for model parameters, which are selected as 0.1θ 0 ≤ θ ≤ 50θ 0 , the setup of the estimation algorithm and its parameters are similar to the previous step (see Section 5.2). Figure 12 shows the time history of the posterior mean and standard deviation (SD) of the three components of the FIM. The initial and final estimates of the ten soil-structure model parameters along with the estimated coefficient of variation (COV) are listed in Table 6. The estimated parameter values with larger COV include larger estimation uncertainties. Figure 13 compares the measured acceleration response time histories at the selected measurement channels with those estimated using the final estimates of the model parameters and FIMs. This figure indicates that there is a remarkable agreement between the estimated and measured acceleration responses. Furthermore, Figure 14 presents the RRMSE of the updated FE model responses at different measurement channels. Comparing Figure 14 with Figure 5, the discrepancies between the FE predictions and measurements in the NS direction are generally less than those in the EW direction. This is predictable since the soil-structure interaction effects are more dominant in the NS direction of the Millikan Library [28], and hence, the flexible-base model uncertainty is expected to be less in the NS direction than in the EW direction.

Effective Structural Stiffness Parameters
The estimated effective elastic moduli of columns and walls (EW&C1 = 19.6 GPa, EW&C2 = 26.1 GPa) are larger than the nominal value that is based on code recommendations, which is 17.3 GPa. Several factors can collectively account for this difference. (i) The effects of lateral stiffness of the backfill soil at the basement level are not modeled in the FE model. This is likely the reason behind the larger estimated value for EW&C1. It can also be a source of modeling error because the columns and walls at the first-story and basement level are parameterized with the same effective elastic modulus. (ii) Concrete aging, which will increase the stiffness of concrete material, is not accounted for in the nominal estimations. (iii) The lateral stiffness of nonstructural components and systems (i.e., partition walls and stairs) are not accounted for in the FE model. This can result in estimating larger values for the effective elastic moduli of the columns and walls, which characterize the lateral stiffness of the building. (iv) Finally, the nominal elastic modulus is calculated for the concrete material without accounting for the effects of steel reinforcement, which can partially explain the larger estimated effective elastic moduli.

Effective Structural Stiffness Parameters
The estimated effective elastic moduli of columns and walls (E W&C 1 = 19.6 GPa, E W&C 2 = 26.1 GPa) are larger than the nominal value that is based on code recommendations, which is 17.3 GPa. Several factors can collectively account for this difference. (i) The effects of lateral stiffness of the backfill soil at the basement level are not modeled in the FE model. This is likely the reason behind the larger estimated value for E W&C 1. It can also be a source of modeling error because the columns and walls at the first-story and basement level are parameterized with the same effective elastic modulus. (ii) Concrete aging, which will increase the stiffness of concrete material, is not accounted for in the nominal estimations. (iii) The lateral stiffness of nonstructural components and systems (i.e., partition walls and stairs) are not accounted for in the FE model. This can result in estimating larger values for the effective elastic moduli of the columns and walls, which characterize the lateral stiffness of the building. (iv) Finally, the nominal elastic modulus is calculated for the concrete material without accounting for the effects of steel reinforcement, which can partially explain the larger estimated effective elastic moduli. Table 7 compares the estimated lumped stiffness and viscosity of the soil subsystem with the NIST recommendations [3] based on [40]. The lumped soil stiffness and viscosity should be evaluated at the fundamental frequency of soil-structure system using the flexible-base dimensionless frequency a 0 = ωb /V s , where ω is the flexible-base fundamental circular frequency, b = 11.5 m is the foundation half-width, and V s is the shear wave velocity, which is averaged over an effective soil profile depth based on the recommendations provided in [3]. Using the shear wave velocity profile of the Millikan Library site presented in [41], the average shear wave velocity is calculated as 398.1 m/s for translation, 392.1 m/s for rocking in the NS direction, 407.7 m/s for rocking in the EW direction, and 422.3 m/s for torsion. The flexible-base frequency is a function of the soil stiffness coefficients as follows [3].

Soil Subsystem Stiffness and Viscosity
Equation (9) presents the flexible-base to fixed-base fundamental period ratio (i.e., period elongation) in the x-direction, where T x is the flexible-base fundamental period, T x is the fixed-base fundamental period, k x is the fixed-base structural stiffness, h is the effective fundamental modal height, K x is the soil horizontal stiffness, and K yy is the soil rocking stiffness. k x is calculated as 411.4 MN/m given the effective fundamental modal mass (m x = 5252.9 ton) and T x = 0.71 s. Likewise, k y is calculated as 1129.9 MN/m given m y = 5541.2 ton and T y = 0.44 s. The effective fundamental modal height is taken as 0.7 times the building height [1] (h = 0.7 × 43.9 = 30.73 m). The process to calculate T is iterative [42]. T is intially utilized to evaluate the soil stiffnesses. Then, T is calculated using Equation (9), and the soil stiffnesses are reevaluated using T. Then, T and the soil stiffnesses are iteratively updated until no considerable change in T is observed. This procedure is used for calculating the soil parameters corresponding to the two translational vibration modes in the NS and EW directions. For the torsional mode, the fixed-base torsional period (T zz = 0.33 s) is used for evaluating the torsional soil parameters. The estimated values in Table 7 are the integrated stiffness and viscosity of the soil subsystem with respect to the center of the foundation pad using the mean estimates of subgrade and viscosity moduli obtained in Section 6.2. To maintain compatibility with previous studies [28], we normalize the sway stiffness and viscosity by GL, where G = 268 MPa and 2L = 27.4 m are the soil shear modulus and a reference foundation length, respectively. For rocking and torsional components, the normalizing factor is GL 3 . Table 7 indicates that the estimated stiffness and viscosity coefficients have a nonnegligible difference with the NIST-based recommended values. One potential reason for the observed difference is the assumption of rigid foundation used in deriving the approximate formulas of the impedance function reported in [3], which was shown to be violated in the Millikan Library, especially in the EW direction [26]. In contrast, the foundation flexibility is explicitly modeled in this study. Furthermore, the estimated impedance function in this study is obtained using the building's seismic response that contains effective high-frequency components, as shown in the recording of Channel #16 presented in Figure 13, which is typical for the response of mid-level and lower-level floors. Since the impedance function is frequency-dependent, the presence of multiple frequencies in the response will likely make it difficult to exclusively tune the estimated impedance function to the flexible-base fundamental frequency. Therefore, some deviation of the estimated impedance function using multiple-frequency responses from the NIST recommendations, in which the impedance function is evaluated at a single frequency, is expected.
Comparing the stiffness coefficients, an especially significant discrepancy can be observed between the estimated and NIST-based recommended values of the translational stiffness in the NS direction and Up-Down direction. A similarly significant discrepancy between the theoretical and experimental estimates of soil translational stiffness in the NS direction was reported in [41], and it was attributed to either an error in the measured foundation translational motion or the crude simplification of the foundation model. To our knowledge, no studies have reported the soil stiffness in the Up-Down direction; therefore, we could not compare our estimated values to other reference values.
Comparing the viscosity coefficients, almost all estimated viscosity coefficients deviate considerably from the corresponding NIST-based values, which is possibly a manifestation of the high uncertainty with which these parameters are estimated, as shown in Table 6. Estimating damping experimentally has always been a key problem in system identification of building structures [43]. This is likely the case with the Millikan Library since there were challenges with accurately determining the phase of the response (with respect to the applied force) in forced vibration tests [41]. Moreover, the estimated viscosities are almost always less than their NIST-based counterparts. This can be explained by examining the soil shear wave velocity profile of the Millikan Library site reported in [41]. Looking at this profile, the shear wave velocity increases almost linearly with depth until it reaches 944.8 m/s at 118.57 m depth. This rapid increase in the shear wave velocity with depth most likely causes multiple wave reflections at the different soil layer interfaces, hence the reduction in soil radiation damping. We can also observe that the reduction in radiation damping due to the nonuniform soil profile is more significant in rocking than in horizontal translation, which is reported in [3]. Thus, we can state that although using an average shear wave velocity can fairly simulate the nonuniform soil stiffness, it cannot simulate the nonuniform soil radiation damping at the same level of accuracy.
The lumped stiffnesses of the soil subsystem for the Millikan Library have also been estimated using other system identification methods in the literature. Table 8 presents the estimated stiffness coefficients in six different previous studies compared to those identified in this study. As can be seen, the translational stiffness in the EW direction is close to the value identified by Luco et al. [26]; the NS translational stiffness is not comparable to any other study; the EW and NS rocking stiffnesses are close to the static values used by Chen et al. [45] and adopted from Balendra et al. [46]; and the torsional stiffness is close to the value estimated by Ghahari et al. [28]. The observed discrepancy between the identified model parameters and the corresponding parameters identified in previous studies highlights the challenging nature of the problem at hand. It also draws attention to the importance of widening the application of structural sensing to real-life structures to develop accurate numerical models for structural health monitoring, condition assessment, and damage prognosis, especially in high-seismicity areas [47]. Table 9 compares the period of the first five structural modes between the updated flexible-base and fixed-base models. The table, moreover, reports the period elongations. The estimated modal periods of the flexible-base structure have good agreement with the results reported in [28].  Figure 15 compares the Rayleigh damping ratio as a function of frequency between the initial, updated fixed-base, and updated flexible-base FE models. Table 6 shows that the estimated mass-proportional Rayleigh damping coefficient for the flexible-base structure is relatively small, which results in a nearly stiffness-proportional linear Rayleigh damping curve for the updated flexible-base model. This is a well-known effect of SSI on effective modal damping [48].

Period Elongation
prognosis, especially in high-seismicity areas [47]. Table 9 compares the period of the first five structural modes between the updated flexible-base and fixed-base models. The table, moreover, reports the period elongations. The estimated modal periods of the flexible-base structure have good agreement with the results reported in [28].  Figure 15 compares the Rayleigh damping ratio as a function of frequency between the initial, updated fixed-base, and updated flexible-base FE models. Table 6 shows that the estimated mass-proportional Rayleigh damping coefficient for the flexible-base structure is relatively small, which results in a nearly stiffness-proportional linear Rayleigh damping curve for the updated flexible-base model. This is a well-known effect of SSI on effective modal damping [48]. Figure 15. The Rayleigh damping ratio as a function of frequency in the initial, updated fixed-base, and updated flexible-base FE models.

Possible Sources of Estimation Uncertainty
The estimated model parameters include uncertainties, which are quantified relatively using the estimated coefficient of variations (COVs). Aside from parameter identifiability issues, which are discussed earlier in the paper, modeling errors are most likely the most important source of estimation uncertainties. Modeling errors result from the Figure 15. The Rayleigh damping ratio as a function of frequency in the initial, updated fixed-base, and updated flexible-base FE models.

Possible Sources of Estimation Uncertainty
The estimated model parameters include uncertainties, which are quantified relatively using the estimated coefficient of variations (COVs). Aside from parameter identifiability issues, which are discussed earlier in the paper, modeling errors are most likely the most important source of estimation uncertainties. Modeling errors result from the inherent imperfections, approximations, and idealizations in the mathematical FE model. As a result, the selected class of FE models does not contain the real-world structure, and the estimated parameters characterize the closest possible model, in the model class, to the real structure [49]. Geometric approximations, material nonlinearities (even during low-amplitude earthquake excitations), effects of nonstructural components and systems, and modeling of damping energy dissipation mechanisms are some examples of imperfections in the numerical FE model. Model parameterization (i.e., selecting the estimation parameters in the model) is another source of modeling error. Modeling errors lead to biased estimation results. This means that by changing the earthquake event, the estimation result may vary, preferably slightly. Addressing the effects of modeling error requires repeating this study using more refined FE models and/or different earthquake data, which is beyond the scope of the present paper and can be the subject of future studies.

Conclusions
In this study, we developed a finite element (FE) model updating framework using a sequential Bayesian estimation approach to identify the soil and structural model parameters as well as the input motions from recorded seismic response of the Millikan Library building during the 2002 Yorba Linda earthquake. The newly developed Bayesian FE model updating procedure was used in a two-step system identification approach. In the first step, a fixed-base structural model was identified using partially unknown inputs. The structural model parameters included the effective stiffness of columns and shear walls (grouped along the height of the structure), the Rayleigh damping parameters, and the distributed floor mass. In the second step, the parameters characterizing the soil-structure model and the time histories of the Foundation Input Motions (FIMs) were estimated. The soil-structure model parameters included the soil subgrade and viscosity moduli at different foundation regions and the Rayleigh damping parameters. For each step, an identifiability analysis based on an information theoretic approach was performed to study the dependencies between model parameters and select the most identifiable parameter sets to be estimated. In Section 7 of the paper, the estimated model parameters were compared with the state-of-practice recommendations and other previous studies on the Millikan Library, and the differences were highlighted.
Minor differences were observed between the measured and FE-predicted response time histories obtained from the updated FE models. Considering the various sources of modeling error in this identification problem, the agreement between the updated FE model responses and the measurements was remarkable. Future studies are suggested to extend this process with more refined FE models (e.g., including material nonlinearities) and different earthquake data to further quantify the effects of modeling error. Applying this process to nonlinear SSI models using the recorded structural responses to strong ground motions may introduce additional complications due to the damage in nonstructural components, the nonlinear kinematic interaction between the structural and nonstructural systems, and the nonlinearities in the soil response. As a result, the identification of the nonlinear structural model parameters and the Raleigh damping model parameters is expected to include uncertainties and estimation biases. Furthermore, the identified soil-foundation subsystem model parameters, i.e., the impedance function, are expected to deviate more considerably from the state-of-practice recommendations. Despite these difficulties, undertaking such a research effort is crucial for developing a better understanding of the nonlinear interaction between the non-structure, structure, and soil subsystems. By utilizing the presented model inversion framework for different civil structures and earthquake records, key structural model parameters can be estimated from real-world seismic data. These are valuable information that can guide-along with theoretical results-the selection of engineering analysis and design parameters to improve state-of-the-art seismic analysis and design procedures.  Data Availability Statement: All data, models, or codes that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.