Optimal implementation of frequency domain impedances in time domain simulations of building structures on embedded foundations

Soil-Structure Interaction (SSI) have been studied the last decades, and proper analysis for the linear elastic case in frequency domain has been established success-fully. However, SSI is rarely considered in the seismic design of building structures. Regardless of its importance as a signiﬁcant source of ﬂexibility and energy dis-sipation, buildings are analyzed using a rigid base assumption, and the design is based on a response spectrum analysis, for which not only the soil, but also time are totally ignored. In a ﬁrst attempt to improve and to incentivize time domain analyzes compatible with standard ﬁnite element packages for the engineering com-munity, the state-of-practice introduces two major simpliﬁcations to transform the frequency domain analysis into a time domain analysis: (a) it assumes the frequency at which the impedance value should be read is the ﬂexible-base frequency, and (b) it also assumes that the foundation input motion preserves the phase of the free ﬁeld motion. Upon these simpliﬁcations, the following questions may arise: How does NIST recommendations perform in overall against a full ﬁnite element model? Are the embedment eﬀects for shallow foundation not important so that the phase angle can be neglected? What is the best dimensionless frequency to estimate the soil impedance? Is it possible to make a better estimation of the dimensionless frequency to increase the NIST accuracy? In this study, we attempt to address these questions by using an inverse problem formulation. sion, Period Elongation, Radiation Damping, Inverse problem, NIST, Reduced-Order Model, Mathematical modeling.


Introduction
Accurate quantification of dynamic soil-structure interaction (SSI) effects is critical in design of earthquake resistant structures. If SSI effects are poorly estimated or even neglected, then an unsafe or overly conservative designs is produced because of an overor under-estimation of the critical response of the structure [1,2]. In general, SSI analyses are carried out by means of either the direct or the substructure methods [3][4][5][6][7]. In the direct method, the near-field soil and the structure are modeled explicitly and artificial boundary conditions, which ideally introduce transparency to both incoming and outgoing waves, are employed to truncate the semi-infinite extent of the far-field soil [8][9][10].
The substructure method [4,11,12] , on the other hand, divides the SSI problem into two sub-systems: the structure subsystem and the soil subsystem represented by impedance functions [13][14][15]. These functions are complex-valued, the real part is the representation of the soil stiffness and inertia while the imaginary part is the representation of geometric damping within the soil. To excite the model properly, one needs to modify the free field motion (FFM) to map the ground motion in-coherency effects along the interface of the structure foundation and the soil. Despite the computational efficiency of the substructure method, one can only take full advantage of the method accuracy in frequency domain analysis which is only suitable for linear elastic problems. This is mainly because the soil impedance functions are nonlinear function of frequency and, in time domain analysis, they need to be convoluted with the so-called foundation input motion (FIM) to map the correct tractions along the structure foundation and the soil.
To make the use of the substructure method practical in time domain analysis, researchers and practitioners have proposed different models to approximate the soil impedance functions [16][17][18]. Among them, and for building structures with shallow foundation, the National Institute of Standards and Technology (NIST) under the project entitled Improved Procedures for Characterizing and Modeling Soil-Structure Interaction for Performance-Based Seismic Engineering recommends using single-valued functions at a representative frequency, which can be modeled as constant-valued springs and dashpots along the soilstructure interface (see Figure 1d). The frequency at which the soil impedance is read is computed using an iterative method proposed by Bielak and Veletsos [1,19]. In this method, one starts assuming that the period of the interaction is the fixed-base building period (i.e., T ). This value and geometry of the foundation allow to compute the dimensionless frequency a 0 = ω B Vs as shown in Figure 1b. Such frequency enables to compute the values of the lumped translational and rotational spring and dashpot coefficients as represented in Figure 1c. Provided with the lumped soil impedance coefficients, one can compute the period elongation of the system using where the variableT is the flexible-base period, T is the fixed-base period, h is the firstmodal height, k is the fixed building stiffness, k xx and k θθ are the lumped translational and rotational springs coefficients that account for the flexibility of the surrounding soil. The process is repeated until there is no variation of the period elongation. Finally, provided with the flexible-base period, one can read from the impedance function the lumped soil spring k xx , k θθ as well as the lumped soil dashpots c xx , c θθ and -if necessary -distribute them along the foundation perimeter as shown schematically in Figure 1d. Furthermore, for modifying the FFM due to the foundation embedment, the NIST recommends using zero-phase transfer functions H u and H θ , which reduce the translational and introduce the rotatational motion, respectively.
where the variable D represents the foundation embedment, V s the shear wave velocity, and ω the frequency. This mapping to determine the FIMs are shown schematically in Figure 2 where the FFM signal in Figure 2a is first transformed in Fourier space (dashed blue line), then the Fourier coefficients are scaled by the translational and rotational transfer functions as shown in Figure 2b (in solid black line), and the results are transformed back to the real space (using the inverse Fourier transform) to obtain the FIM as shown in Figure 2c in red solid lines. It is evident that the NIST recommendations make two important simplifications for using the substructure method in time domain analysis: (a) the frequency at which the impedance value should be read is the flexible-base frequency, and (b) it assumes that the FIM preserves the phase of the FFM and it does not change as a function of foundation embedment ratio D/B. Upon these simplifications, the following questions may arise: 1. How does the NIST recommendations and the proposed model perform in overall?
2. Are the embedment effects for shallow foundation not important so that the phase angle can be neglected?
3. What is the best dimensionless frequency a 0 to read the soil impedance?
4. Is it possible to make a better estimation of the dimensionless frequency a 0 to increase the NIST accuracy?
The main objective of this study is to address the above listed questions, and to the best of our knowledge, this is the first comprehensive study in this regard. In the rest of this manuscript, we provide details of the approach we used to assess the performance and optimality of the NIST recipe for modeling SSI effects for building structures. This calls for (1) computing the ground truth (reference) solution, and searching the parameter space of the reduced order model to find the optimal impedance function values or the frequency at which one needs to read them. For the former, we model the problem using the direct method and for the latter we use a Bayesian approach based on the ensemble Kalman filtering; details of each are provided in §2- §4. Then, in §5, we use the developed framework to comprehensively assess the predictive capability of the NIST recommendations for a series of building structures embedded in elastic half-space. We summarize the findings of this study in §6.
2 Direct modeling of the problem For direct modeling of the building structure response on elastic half-space -see Figure 1a we use the procedure shown in Figure 3. We truncate the semi-infinite half-space along the red line shown in Figure 3a, and introduce Lysmer dashpots accompanied with prescribed nodal forces to artificially approximate an interface transparent to both incoming and outgoing waves [8,20]. The Lysmer dashpot coefficient shown in Figure 3b for absorbing outgoing waves are computed as: Similarly, the introduced nodal forces, shown in Figure 3c and d, for translation of the vertically propagating incoming shear waves within the domain of interest are In Equation (4) and (5) the variable t h = 1.0 [m] is the out-of-plane thickness of the truncated soil domain, ∆l i is either ∆x i or ∆z i and represent the tributary i-th side length of the element which are different for corner and inner boundary nodes, V s is the shear wave velocity of the soil, and V p is the compressive wave velocity in the soil. It worth to point out that the boundary conditions enforced in this manner can almost perfectly absorbs body waves with angles of incidence greater than 30 • degrees with respect to the vertical axis, and will lose their effectiveness for lower angles of incidence or for surface waves. Moreover, in Equation (6)-(8) the signalsu g (t) andu s i (t) represent the incident and the total free field motion for which a Ricker-wavelet is prescribed at the bottom of the model and given in Equation (9) as: where β = (πf 0 ) 2 , f 0 is the characteristic frequency, t 0 is the time position where the velocity will become maximum. In all simulations, we consider a characteristic frequency f 0 = 2.0 [Hz], and a peak time velocity t 0 = 1.0 [s]. The node velocityu s i (t) andu r i (t) needed to compute the horizontal and vertical forces on the left and right boundary areu and,u where z i is the vertical coordinate of the i-th node measured with respect to the bottom of the soil domain, H is the total vertical height of the truncated soil domain, and ± sign in Equation (8) Figure 4 for both components of the velocity field. This figure makes evident that the implemented absorbing boundary interface works perfectly for this case. Here, one can observe that there is no reflection coming in from the boundaries, and as expected no velocities in the vertical direction are developed throughout the simulation.

Reduced order modeling of the problem
For modeling the problem using the substructure method, we consider the reduced order model (ROM) shown in Figure 5c. The ROM is constructed in such a way that it preserves the modal information of the fixed base building. To this end, we employ frame elements with three degree of freedom per node to represent the structure's geometry as shown in Figure 5a. We assume that each floor acts as a rigid diaphragm, so that the buildings mass can be lumped at the floor levels as shown in Figure 5b. We then use static condensation [21,22] on the fixed-base building to compute the M s , C s , K s ∈ R n×n matrices where n is the number of floors; therefore, the vector x ∈ R n represents the horizontal degree of freedom of the building, since there is only one translational degree of freedom on each floor. In addition, we assume a rigid rectangular foundation of half-width B and depth D sitting on an elastic half-space, for which two additional degree of freedom u, θ ∈ R appear because of the flexibility of the soil. The foundation has a total mass m f , and a moment of inertia I 0 . Then, depending on the motion that one can use to excite the model, we consider two variants for forward modeling: (1) excited by FIM and (2) excited by FFM. The semi-discrete equations of motion for both cases are shown below.
where the external force vector for the free field motion is and for the foundation input motion is andẍ,ẋ, x ∈ R n are the acceleration, velocity and displacement vector for the condensed horizontal degrees of freedom of the building, alsov,v, v ∈ R andθ,θ, θ ∈ R are the foundation horizontal and rotational acceleration, velocity and displacement. Note that we define 1 ∈ R n the vector of ones, this is 1 = (1, 1, . . . , 1), and h ∈ R n the vector of height, this is h = (h 1 + D, . . . , h n + D). We finally define O n×m ∈ R n×m the matrix of zeros. In addition, it is important to highlight that in the ROM the horizontal and vertical springs and dashpots elements are distributed uniformly over the foundation perimeter. However, the contribution can be lumped as follows for both stiffness and damping components: Distributing the spring and dashpot elements as represented in Figure 5 allows us to take the coupling effects in the stiffness and damping matrices into account approximately. It worth mentioning that the NIST recommendations assume the coupling terms are zero.

Ensemble Kalman inversion for parameter estimation
In order to find the optimal spring and dashpot coefficients of the ROM, we use the Bayesian approach based on the ensemble Kalman inversion (EnKI) [23][24][25], which can be also considered as a particle-based derivative-free sequential optimization method. In the inversion setting, we consider the problem of finding u ∈ R n from y ∈ R m where The variable u ∈ R n consists of all the unknown parameters that we want to estimate, the variable y ∈ R m consists of the ground truth quantities of interest, and η is a zeromean Gaussian noise with covariance Γ. The nonlinear function (a.k.a. forward model) G : R n → R m maps the parameter space to the data space. In this paper, we work with one type of data-sets; the displacement time series recorded at the floors and foundation levels computed using the model described in §2. Given N particles u (n) j , n = 1, . . . , N within the ensemble and at each iteration j, we use the predictions G(u (n) j ) by the forward model and the observation data y to update the particles for iteration j + 1 (see Figure 6 for schematic illustration). That is, j )) for n = 1, . . . , N .
In Equation (18), y (n) j+1 can be either identical to y (the observation data) or randomly perturbed using zero-mean Gaussian noise η; C uw j+1 and C ww j+1 are empirical covariance matrices that can be computed at each iteration based on predictions and the ensemble meanū j+1 using the following equations: It worth mentioning here the stopping criterion is defined as either reaching a relative change of 0.001 in all the parameter in two consecutive EnKI iterations or a maximum number of 500 iterations. Furthermore, the initial ensemble mean is defined to be equal to the lumped stiffness and dashpot coefficients obtained using the NIST [18] procedure. In addition, since the prior distribution of the soil parameter is unknown, we use a uniform distribution centered at the computed NIST value with lower and upper limit ten times smaller and larger to represent this uncertainty. Finally, the positiveness of the stiffness and dashpot coefficients are enforced through a change of variables u ∈ R, exp(u) : u → [0, ∞] on which the EnKI is applied.

Model specifications
Three different topologies of buildings are considered, which are illustrated from the shortest in Figure 7a to the tallest in Figure 7c.

Estimated soil impedances
To estimate the constant valued soil impedances, the ROM is subjected to the FIM and the EnKI is employed to estimate the soil spring and dashpot coefficients. Then, they are compared against the true impedance functions derived numerically using the method proposed by [15]. These results, named hereafter as EnKI-FIM, are displayed in Figure 8 and 9 in which the vertical axis has been normalized such thatk xx = π G s (k xx + ω c xx i), k xθ = π G s B(k xθ + ω c xθ i), andk θθ = π G s B 2 (k θθ + ω c θθ i), where i = √ −1, B is the half-width foundation, G s = ρ s V 2 s is the soil shear modulus, ρ s the soil's density, and V s the soil shear wave velocity.
We also employ the EnKI to re-estimate the soil springs and dashpots when the system is subject to FFM; the results named EnKI-FFM hereafter. Consdiering the EnKI-FIM results as the optimal solution for the considered ROM, we assessed the accuracy of the EnKI-FFM results as well as those recommended by NIST. In Figure 10, the horizontal axis is the EnKI-FIM results while the vertical axes are the EnKI-FFM and NIST results. This comparison study delineated that 1. The embedment effect is almost negligible as the EnKI-FIM and EnKI-FFM results are highly correlated with each other, which in turn confirms that the embedment  ratios considered in this study are shallow.
2. The NIST results capture very well the total lumped horizontal k xx and rotational k θθ spring coefficients. However, the total lumped horizontal c xx and rotational c θθ dashpot coefficients are poorly estimated. In particular, the lumped rotational dashpot coefficients for the case of building with T = 1.5 [s] are the worst. This problem is due to the fact that the imaginary component of the lumped rotational impedance c θθ becomes small as the fixed building period becomes larger, and then when such dimensionless frequency is used to compute the lumped rotational impedances generates a deviation about one order of magnitude smaller when it is compared against EnKI-FIM results, see Figure 10 for more details. This suggests that the frequency at which flexible buildings irradiates energy due to rocking mechanism is different from the sway (flexible base) mechanism proposed by NIST. Figure 10: Identified impedance component comparison between NIST and EnKI using FIM and FFM. The horizontal axis correspond to the EnKI once the system is subject to FIM. The vertical axis represents the impedance components using NIST (in blue solid dots) and EnKI subjected to FFM (in red solid dots). The black solid/dashed lines represent a linear regression for either case.
The previous analysis made clear that the flexible-base period may not be the best frequency to use to estimate the imaginary component of the soil impedances. In particular, if the flexible-base period is used to compute the lumped rotational dashpot coefficient for flexible buildings, the radiation due to damping is underestimated. Therefore, we would like to investigate what practical dimensionless frequency (based on the flexible-base pe-riodT , fixed-base period T , or the input signal dominant frequency f 0 ) one should use to read the impedance functions. In this section, as before we perform similar analyses, but now we consider only the NIST recommendation using these three frequencies. As shown in Figure 11, we compare on the vertical axis the soil spring and dashpot coefficients obtained with NIST for these three cases, against EnKI-FIM results on the horizontal axis. Once again, we notice that the horizontal and rotational spring coefficients are estimated properly, no matter which frequency is used. However, a smaller variability is obtained if the fixed-base period is employed. On the other hand, a similar behavior, with a large variability, for the lumped horizontal and vertical dashpot coefficients is obtained. In general, the input signal dominant frequency f 0 as well as the fixed-base period T performs better in estimating the lumped vertical dashpot coefficient. The main differences in time history responses comes from the fact (once again) that the lumped rotational dashpot Figure 11: Impedance comparison between NIST using the flexible-base periodT , fixedbase period T , or the input signal dominant frequency f 0 and EnKI using foundation input motion.
coefficient is poorly estimated. In order to show this point, we quantify the total error made in a dynamic analysis between the response of the full FEM (as described in §2) and the ROM with spring and dashpot coefficients obtained using NIST for the three different representative frequencies and using EnKI-FIM. We measure the 2 -error for each case as: where y (k) j represents the k-th observation (from the FEM) at j-th time step,ŷ (k) j represents the k-th response (from the ROM) at j-th time step, N t corresponds to the number of time steps, and N m represents the number of observations considered. We consider a time series of 15 [s] with time step ∆t = 0.01 [s] for the translation and rotation of the foundation, and the horizontal translation of each floor. Figure 12 represents the total 2 -error in the vertical axis, and the model number in the horizontal axis. This figure shows that in general the NIST model using the fixed-base frequency and the signal predominant frequency have smaller errors when they are compared with the flexible-base case. The latter is smaller only for the building with period T = 1.0 [s] and embedment D = 5.0 [m]. However, none of the configurations reached the error obtained using EnKI-FIM coefficients. The discrepancies are attributed mostly to the poor estimation of the rotational dashpot coefficients, which generates more oscillations in the transient part, and they take longer to be damped out in the free oscillation part. Since the estimation of this coefficient is slightly better for the signal dominant frequency and the fixed-base frequency when is compared to the EnKI, we obtain errors that are smaller. Figure 12: Total time history error measured in 2 -norm for the three practical dimensionless frequencies mentioned in NIST. In blue dots the flexible-base period, in red-dots the fixed-base period, and in yellow-dots the signal-predominant frequency.

Dynamic Global Properties
Provided with the soil impedances using NIST and EnKI, we are interested in computing the period elongation ratioT /T and foundation-damping β f that the ROM experiences once is supported on lumped spring and dashpot elements. Since the proposed ROM in §3 considers the coupling between the translation and rotational degree of freedom approximately, we employ Equations (28) and (39) derived in Appendix A to estimate the period lengthening and foundation damping. We then compare these results against the iterative method summarized in §1. (14) show the period elongation and foundation damping computed using the NIST and EnKI-FIM results. As shown, a good-agreement between the proposed expressions in Appendix A and NIST expressions using the iterative method is achieved. From these figures we note that Equation (28) provides with more flexible system, while Equation (39) provides with less-dissipative system when they are compared with the NIST recommendations. A very close fit is however obtained for the three-buildings when the foundation aspect ratio is small (i.e., D B = 0.10). The discrepancies must be attributed to the following facts:  1. The estimated frequency of the interaction using the EnKI is different from the one obtained using the NIST recommendations.

Figures (13) and
2. The iterative method uses the first modal information to compute the lumped spring and dashpot coefficients while the EnKI-FIM estimates implicitly compensates for higher mode contribution since the soil parameters are estimated from the observations of the full FEM.
3. As discussed in §3, we are explicitly considering the spring and dashpot coupling terms in the stiffness and damping matrices. It can be seen in both Equations (28) and (39) that mentioned coupling in the stiffness and damping matrices generate a slight increase in the period elongation ratio as well as a slight decrease in the radiation damping.
Regardless of the discrepancies presented, the results of the global behavior are very similar as shown in Figures 13 and 14. It is not a surprise that the period elongation ratios are well-captured since NIST estimates well the lumped translational and rotational springs.
On the other hand, it is interesting to see that the foundation damping for an equivalent one degree-of-freedom system is also well captured, this is due to the fact that the contribution of the lumped rotational dashpot is not important to the global foundation damping. These differences become more pronounced when the stiffer buildings rest on a deeper foundation, in which case the coupling term may become more important for the foundation damping.

Can NIST do better?
To find the optimal frequency to read the impedance functions as in NIST procedure, we use the EnKI algorithm to estimate the optimal dimensionless frequency a 0 , so that the misfit between the FEM model and the NIST ROM is minimized. In this analysis we consider two configurations: (a) only one frequency a * 0 is considered for the lumped soil impedances, and (b) two frequencies a x 0 (horizontal) and a θ 0 (rotational) for the lumped soil impedances are estimated. The results of this analysis are shown in Figure 15. Figure 15 (a) (b) Figure 15: Optimal dimensionless frequency a * 0 obtained using EnKI for (a) One dimensionless frequency, and (b) Two dimensionless frequency.
reveals two interesting results. The first confirms the hypothesis that the dimensionless frequency that controls the response of the SSI is related to the rotational soil impedance. This can be seen since the dimensionless frequency a * 0 in Figure 15a and a θ 0 in Figure 15b are very similar. The second result suggests that the optimal a * 0 is not only a function of the building itself (T ), but also a function of the embedment ratio and the signal frequency content. In this regard, the dimensionless frequency a 0 ∼ B V sT may not be good enough to estimate the frequency independent impedance, a form like a 0 ∼ h V s T , h B , D B , B f 0 V s , ν s may be more suitable to carry out this task; however, there is no clear pattern on how this dependency could be taken into account. Therefore, the iterative method to obtain the flexible-base period may generate large errors specially when the soil is stiffer, which is consistent with results presented in Figure 12. This does not mean it is a bad approximation for the engineering practice, it suggests that the frequency at which each lumped spring and dashpot may respond are completely different from those proposed by NIST.

Conclusions
The state-of-practice introduces two major simplifications to transform the frequency domain analysis into a time domain analysis. It assumes that (1) the frequency at which the impedance value should be read is the flexible-base frequency and (2) the foundation input motion is preserves the phase of the free field motion. Both simplifications are however adopted in order to perform time history analysis using standard finite element packages. Unfortunately, the recommendations proposed in NIST only work in linear elastic soils, and they are an over simplification of the substructure method (described in [3][4][5]), since do not consider the frequency dependancy of the excavated soil to perform a rigorous time domain analysis.
In this study, we found that these recommendations could capture the global dynamic behavior of building structures relatively well. However, their performance deteriorates in capturing the response of the structure such as the floor displacement time-series, and this is mainly due to the poor approximation of the lumped rotational dashpot coefficients. Consequently, the foundation of buildings with large period tends to be less dissipative since the radiation due to rocking is not properly estimated. Our results in Figure 12 show that the fixed-base period generates better results, however more test should be conducted to confirm the claim. Furthermore, the reader should be aware that the conditions in which NIST recommendations were tested are the best case scenario, in which not only the impedance functions are correctly estimated, but also the foundation input motion was prescribed so that it has the right phase; thus, the discrepancies are only attributed to the impedance frequency estimation.
On the other hand, we showed that the foundation input motion effects for negligible for cases considered here with embedment ratios smaller than 0.5. Therefore, the free field motion can be directly applied without generating a substantial error. Finally, we notice that the selection of the dimensionless frequency a 0 to read the impedance values depends on other dimensionless parameters such as the frequency of the input signal and the embedment ratio. Future work should aim to evaluate the inclusion of these parameters in determination of the system representative frequency.

A Period Elongation and Radiation Damping Derivations.
Let us first, consider two different systems. The first, depicted in Figure 16a of one floor with heighth = h + D, stiffness k ∈ R + , structural damping β i ∈ R + , supported by a distributed horizontal spring k x ∈ R + and distributed vertical springs k z ∈ R + . The second systems, depicted in Figure 16b is a fixed-base single-degree-of-freedom system with mass m ∈ R + , modifiedstiffnessk ∈ R + , modified-fundamental periodT ∈ R + , and modified damping β 0 ∈ R + as it is presented in [16,31].
In the case of the equivalent fixed-base system, the total static displacement∆ generated by an applied external load F is computed as: Figure 16: Model reduction from (a) the simplified flexible-base system into (b) the equivalent fixed-base system.
In this regard, the fixed-base period of the system becomes In addition, the total displacement of the reduced-order-model, when an external load F is applied, becomes Note that, in order to equate the stiffnessk of the fixed-base system to that of the flexible-base system, the displacements must be equal, i.e.,∆ = ∆. By comparison of equations (23) and (25), we conclude that the following relation between the stiffness of the two systems must hold Therefore, the corresponding fixed-base period will be given as In turn, the period elongation of the reduced-order-model T relative to the period of the fixedbased system,T , is defined as follows: