Optimization of frequency domain impedances for time-domain response analyses of building structures with rigid shallow foundations

The effects of dynamic soil–structure interaction (SSI) have been extensively studied in the last few decades, and proper analysis for the linear elastic case in frequency domain has been established successfully. However, SSI is rarely considered in the design of building structures, and instead, buildings are frequently analyzed using a rigid base assumption and quasi-static loading conditions that ignore SSI and its dynamic nature. Acknowledging these shortcomings, the National Institute of Standards and Technology (NIST) published in 2012 a set of recommendations on time-domain analyses of SSI for building structures compatible with standard finite element packages for consideration in engineering design. The so-called NIST GCR 12-917-21 report introduced a major simplification to enable frequency domain tools to be implemented in time domain analyses. That is, replacing the frequency-dependent soil impedance functions by a single-valued functions read at the flexible-base structure frequency; This work seeks to quantify the accuracy of this simplification considering fully coupled two-dimensional (2D) finite element models (FEM) as the reference. Using a Bayesian approach based on ensemble Kalman inversion (EnKI) and a range of numerical simulations of soil–foundation–building interaction, we estimate the optimal frequency that can be used to estimate soil impedance for time domain analyses; and we evaluate the improvement that the corresponding impedance offers relative to the full FEM results when compared to time domain analyses performed in accordance to the NIST recommendations outlined above.

The substructure method (Bielak, 1974;Parmelee, 1967;Wolf and Obernhuber, 1985), on the other hand, divides the SSI problem into two sub-systems: the structure and the soil subsystem, which is represented by impedance functions (Gazetas, 1991;Pais and Kausel, 1988;Seylabi et al., 2016;Veletsos and Meek, 1974). 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 more realistically, one needs to modify the free field motion (FFM) to map the coherency effects introduced by the differential compliance of the foundation and the surrounding 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 functions of frequency and, in time domain analyses, they need to be convolved with the so-called foundation input motion (FIM) to map the correct traction along the interface between structure foundation and soil.
To take advantage of the substructure method in time domain analyses, researchers and practitioners have proposed different models to approximate the soil impedance functions (Stewart et al., 1999a(Stewart et al., , 1999b(Stewart et al., , 2012. For building structures on shallow foundations as shown in Figure 1a, 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 recommended using single-valued functions at a representative frequency, which can be modeled as constant-valued springs and dashpots distributed along the soil-structure interface (see Figure 1b). This representative frequency is computed using an iterative method proposed by Bielak (1974); Veletsos and Nair (1975): first, as shown in Figure 1c, the fixed-base building period (i.e. T) and building structure geometry are used to compute the dimensionless frequency a 0 = vB=V s , where v is the angular frequency of the interaction, B is the half-width, and V s is the shearwave velocity. The lumped translational and rotational spring and dashpot coefficients are next evaluated at this frequency, as shown in Figure 1d. Next, the lumped soil impedance coefficients are used to compute the period elongation of the system as follows: where the variableT is the flexible-base period, T is the fixed-base period, h is the firstmodal height, k is the fixed building first modal stiffness, and k x and k yy are the lumped translational and rotational spring coefficients that account for the flexibility of the surrounding soil. The process is repeated until the change in period elongation estimated in two consecutive iterations is below a certain threshold. Finally, the foundation impedance function-the lumped soil springs k x , k yy and dashpots c x , c yy -are evaluated at the converged flexible-base period and if necessary, are distributed along the foundation perimeter as shown schematically in Figure 1b. The period elongation equation adopted in NIST has been introduced by Veletsos and Meek (1974), Bielak (1974), and Wolf (1985) and unfortunately does not include the coupled rotational-translational foundation stiffness k xy . This term becomes important for deep embedded foundations as presented by Gazetas (1996). Therefore, other expressions (e.g. Maravas et al., 2014;Zania, 2014) can be used in a similar manner to improve NIST accuracy.
To modify the FFM due to vertically propagating shear waves that take the foundation embedment effects into account, NIST recommends using the zero-phase transfer functions H u and H yy originally proposed by Kausel et al. (1978) and Day (1978), which reduce the translational and introduce rotational motion to the FIM, respectively.
In Equations 2 and 3, the variable D represents the foundation embedment, V s the shear wave velocity, and v the frequency. Figure 2 shows schematically the FFM to FIM mapping, where the FFM signal in Figure 2a is first transformed in the Fourier space (dashed gray 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 then transformed back to the real space (using the inverse Fourier transform) to obtain the FIM as shown in Figure 2c. Finally, the resulting FIM is applied to a rigid diaphragm (a.k.a bathtub) as shown in Figure 1b.
It is evident that the NIST recommendations make two important simplifications to enable the substructure method to be used in time domain analyses: (a) the impedance value is evaluated at the flexible-base frequency and (b) the FIM preserves the phase of the FFM, independent of the foundation embedment ratio D=B. Upon these simplifications, the following questions may arise: 1. How do time-domain analyses using the NIST procedure compare to the direct SSI method? 2. What is the optimal impedance and corresponding dimensionless frequency a 0 for implementation in time-domain analyses? The main objective of this study is to address the above listed questions. In the rest of this manuscript, we provide details of the approach we used to assess the performance of NIST recommendations for time domain modeling of SSI effects for two-dimensional (2D) building structures embedded in elastic half-space. We explain in detail the ground true (reference) solution; the derivation of a reduced order model (ROM) for inverse analyses; and the search of the parameter space of the ROM to the optimal impedance function values and the corresponding frequency. For the former, we model the problem using the direct method and for the latter we use a Bayesian approach based on the EnKI. Then, we use the developed framework to comprehensively assess the predictive capability of the NIST recommendations for a series of building structures. Finally, we summarize the findings of this study.

Direct modeling of the problem
For direct modeling of the 2D building structure embedded in elastic half-space, we use the procedure described in Figure 3. In the following sections, we first describe the building and soil finite element model (FEM), and then the far-field boundary conditions.

FEM
Three different topologies of buildings are considered, which are illustrated from the shortest in Figure 4a to the tallest in Figure 4c. The first reinforced concrete building has columns of rectangular cross section 0:9030:60 m, each of 3:5 m height. Beams have rectangular cross sections 0:8030:60 m and 6:0 m length. The reinforced concrete density is 2500 kg=m 3 , the reinforced concrete elasticity modulus is taken as 25 GPa, and beams are subjected to a dead load of 3600 kg=m. This configuration gives a total height of 21 m, a total mass of 6:45 Á 10 5 kg, a fixed-base fundamental period of 0:507 s, a firstmodal height of 14:91 m, and a first modal mass of 0:545 Á 10 6 kg. The solid core foundation has a half-width of 10 m, an equivalent out-of-plane thickness of 1:0 m, and three different embedment depths of 1:0, 2:5, and 5:0 m are considered, respectively. The reinforced concrete material properties for the foundation are taken such that the density provides a foundation mass 50:0 Á 10 3 kg, elasticity modulus 30 GPa, and Poisson's ratio of 0:25. In this regard, a thick rigid mat foundation is considered and the concrete density is modified so that it emulates the basement levels. The foundation rigid diaphragm behavior is enforced by adding an auxiliary node with three degrees of freedom (i.e. diaphragm node) and then impose kinematic constraints to all nodes at the foundation and columns as suggested by Bathe (2006) and Cook et al. (2007). The second and third reinforced concrete buildings only differ in terms of the column cross-sections, which for the taller buildings are 1:0030:80 m. Beams and material properties of building and foundation are same as above. The foundation is assumed to be a rigid diaphragm and to have perfect bonding with the surrounding soil. Table 1 summarizes some of the buildings' characteristics.  The soil is represented as an elastic, homogeneous, and semi-infinite medium in plane strain setting with density r s = 2000 kg=m 3 , Poisson's ratio n s = 0:25, and shear wave velocity V s = 80, 100,125,150,175,200,225,250,300,400, 500 m=s. The time step that satisfies stability constraints is Dt = 0:001 s, and the simulation time is set to be t sim = 30 s. The spatial discretization is then taken as Dx = Dz = 0:5 m, and the soil domain is truncated so that L x = 350 m and L z = 200 m. The resulting finite element mesh has approximately 280, 000 quadrilateral elements; we employed parallel computing for the simulations using an inhouse finite element code Seismo-VLab (Kusanovic et al., 2020); we divided the domain into seven sub-domains using the computer software OpenMPI 1 (Broquedis et al., 2010;Gabriel et al., 2004;Hursey et al., 2007). The average time per simulation was 65 min approximately.

Far-field boundary conditions
We truncate the semi-infinite half-space as 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 (Assimaki, 2004;Lysmer et al., 1969). The Lysmer dashpot coefficient shown in Figure 3b for absorbing outgoing waves is 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 Equations 4 and 5, the variable t h = 1:0 m is the out-of-plane thickness of the truncated soil domain, Dl i is either Dx i or Dz i and represent the tributary ith 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 is worth pointing out that the boundary conditions enforced in this manner can almost perfectly absorb body waves with angles of incidence greater than 308 with respect to the vertical axis, but lose their effectiveness for lower angles of incidence, including surface waves. Moreover, in Equations 6 to 8, the signals _ u g (t) and _ u 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 b = pf 0 ð Þ 2 , f 0 is the characteristic frequency, and 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 characteristic frequency f 0 is selected to generate an FFM with frequency content that ranges between 0 À 7:5 Hz and emulates the majority of real earthquake signals. Base slab averaging effects are not included in the simulations since the FFM consists of vertically propagating waves. To compute the horizontal and vertical forces on the left and right boundaries, the nodal velocities _ u s i (t) and _ u r i (t) are calculated as follows: and, where z i is the vertical coordinate of the ith node measured with respect to the bottom of the soil domain, H is the total vertical height of the truncated soil domain, and 6 sign in Equation 8 represents the left ( + ) or right (À) vertical boundary forces, respectively.
We verify our implementation by reproducing the free-field motion in a homogeneous half-space with length L x = 350 m and height L z = 200 m. We consider linear elastic soil elements with r s = 2000 kg=m 3 , n s = 0:25, t h = 1:0 m, for which plane-strain conditions are enforced. The shear wave velocity is set to V s = 225 m=s, leading to a wave-length of l n = V s =f 0 = 113 m corresponding to the dominant frequency of the wavelet. The results of the normalized velocity fieldv(x,ŷ) = v(x,ŷ)= max jvj as a function of the normalized coordinatesx = x=l n ,ŷ = y=l n are displayed in Figure 5 for both components of the velocity field. This figure makes evident that the implemented boundary interface works perfectly for this case. Here, one can observe that there are no reflections coming from the boundaries, and as expected, velocities in the vertical direction do not manifest.

Reduced order modeling of the problem
For the substructure method, we consider the 2D ROM shown in Figure 6c.
The ROM is constructed to preserve the modal characteristics of the fixed base building (Asimaki et al., 2019;Kusanovic et al., 2018). To this end, we employ frame elements with three degrees of freedom per node to represent the structure's geometry as shown in Figure  6a. 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 6b. We then use static condensation (Guyan, 1965;Irons, 1965) on the fixed-base building to compute the M s , C s , K s 2 R n3n matrices where n is the number of floors; therefore, the vector x 2 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 embedment depth D on an elastic half-space, for which two additional degrees of freedom u, u 2 R represent the translation and rotation of the foundation resulting from the soil compliance. The foundation has a total mass m f and a moment of inertia I 0 . Then, depending on the excitation, 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: k xy k yy 2 6 4 3 7 5 x v u where the external force vector for the free field motion is and for the foundation input motion is and € x, _ x, x 2 R n are the acceleration, velocity, and displacement vector for the condensed horizontal degrees of freedom of the building, also € v, _ v, v 2 R and € u, _ u, u 2 R are the foundation horizontal and rotational acceleration, velocity, and displacement. Note that we define 1 2 R n the vector of ones, this is 1 = (1, 1, . . . , 1), and h 2 R n the vector of height, this is h = (h 1 + D, . . . , h n + D). We finally define O n3m 2 R n3m the matrix of zeros.
We should highlight here that in the ROM, the horizontal k i x and vertical k i z 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 6 allows us to take the coupling effects in the stiffness and damping matrices into account. We next describe how the stiffness and damping coefficients are computed using the NIST procedure and the methodology proposed in this article.

NIST methodology
The NIST methodology presented in Stewart et al. (2012) only provides expressions for impedance functions and FIM transfer functions in the three-dimensional (3D) setting. However, the same procedure can equally be applied to 2D plane-strain cases using the impedance functions and FIM transfer functions for the 2D case. To make a fair comparison, we first derived the ''true'' 2D impedance and FIM transfer functions in a planestrain setting using the methodology proposed by Seylabi (2016) and Seylabi et al. (2016). The numerically computed impedance functions and FIM transfer functions are presented in Figures 7 to 9, respectively. Please note that we only use k x and k yy components within the NIST framework and read them at the representative frequency considered for analysis. However, it should be noted that we use both the real and imaginary components (or equivalently the amplitude and phase angle) of H u and H yy to compute the FIM to be used in Equation 14. In Figure 7, we only represent the amplitude of the TFs, and we verify them against results presented in Conti et al. (2018). With these considerations, we here assess the NIST procedure's errors associated with reading impedance functions at a nonrepresentative frequency as well as those associated with ignoring the coupling terms.

EnKI-based methodology
For Equation 12 to provide accurate results, that is, closest responses to the direct modeling approach, the spring coefficients k i x , k i z and dashpots coefficients c i x , c i z need to be computed within an optimization framework. 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) (Evensen, 2006;Iglesias et al., 2013;Law et al., 2016), often referred to as particle-based derivative-free sequential optimization method. In the inversion setting, we consider the problem of finding u 2 R n from y 2 R m where The variable u 2 R n consists of all the unknown parameters that we want to estimate, the variable y 2 R m consists of the ground truth quantities of interest, and h is a zero-mean  Gaussian noise with covariance G. The nonlinear function (a.k.a. forward model) G : R n ! R m maps the parameter space to the data space. In this article, we work with one type of observation: the displacement time series recorded at the floors and foundation levels computed using the direct modeling method. 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, that is, the displacement time series computed by the ROM, and the observation data y to update the particles for iteration j + 1 (see in Appendix 1Figure 18 for schematic illustration). That is, In Equation 18, y (n) j + 1 can be either identical to y (the observation data) or randomly perturbed using zero-mean Gaussian noise h; 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 u j + 1 using the following equations: and We should note here that the convergence criterion was defined as either reaching a relative change of 0:1% in all the parameters in two consecutive EnKI iterations, or a maximum number of 500 iterations. Note that the initial ensemble mean is set equal to the lumped stiffness and dashpot coefficients obtained using the NIST (Stewart et al., 2012) procedure. In addition, since the prior distribution of the soil parameters is unknown, we use a uniform distribution centered at the computed NIST value with lower and upper limits 10 times smaller and larger to represent this uncertainty. Finally, the positiveness of the stiffness and dashpot coefficients is enforced through a change of variables u 2 R, exp (u) : R ! ½0, ' on which the EnKI is applied. The ROM constructed using EnKI and subjected to FFM or FIM are named hereafter as EnKI-FFM or EnKI-FIM, respectively.

Soil impedances of ROMs
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 2D impedance functions derived numerically using the method proposed by Seylabi et al. (2016). These EnKI-FIM results are displayed in Figures 8 and  9 in which the vertical axis has been normalized such thatk x = p G s (k x + v c x i), k xy = p G s B(k xy + vc xy i), andk yy = p G s B 2 (k yy + vc yy i), where i 2 = À 1, B is the half-width foundation, G s = r s V 2 s is the soil shear modulus, r 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. Considering the EnKI-FIM results as the optimal solution for the considered ROM, we next assess the accuracy of the EnKI-FFM results as well as those obtained using NIST procedure. 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 x and rotational k yy spring coefficients. However, the total lumped horizontal c x and rotational c yy dashpot coefficients are not equally well estimated. The least well estimated are the lumped rotational dashpot coefficients for tallest building, with T = 1:5 s. This problem arises because the imaginary component of the lumped rotational impedance c yy is inversely proportional to the fixed building period; consequently, for very low dimensionless frequencies, the imaginary component of the lumped rotational impedances is approximately one order of magnitude smaller compared to the EnKI-FIM results, as shown in Figure 10. This suggests that the frequency at which flexible buildings radiate energy due to rocking is different from the swaying mechanism proposed by NIST, which is associated with the flexible base response.
The previous analysis made clear that one can improve upon the frequency recommendation used to estimate the imaginary component of the soil impedances. We specifically showed that if the flexible-base period is used to compute the lumped rotational dashpot coefficient for flexible buildings, the radiation damping is underestimated. In the remainder of this section, we seek the optimal dimensionless frequency (based on the flexible-base periodT , fixed-base period T, or the input signal dominant frequency f 0 ) that can be used to evaluate the foundation impedance functions. As before, we perform numerical analyses to answer this question, but now we consider only the NIST recommendation using these three frequencies.
In Figure 11, we compare the NIST-recommended soil spring and dashpot coefficients (vertical axis) against the 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. However, the lumped horizontal and vertical dashpot coefficients show high sensitivity to the building and foundation geometry and stiffness. In general, the input signal dominant frequency f 0 as well as the fixed-base period T perform better in estimating the lumped vertical dashpot coefficient. Here as well, the main difference in the time history response is associated with the lumped rotational dashpot coefficient that is poorly estimated. In order to show this point, we quantify the total deviation of the dynamic response of the full FEM model (as described in direct model section); and the two ROMs, one with the NIST-recommended spring and dashpot coefficients, and one with the EnKI-FIM impedances, for the three structures of interest. We measure the ' 2 -error for each case as: where y (k) j represents the kth observation (from the FEM) at jth time step,ŷ (k) j represents the k-th response (from the ROM) at jth 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 a sampling time step Dt = 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 shear wave velocity in the horizontal axis.
2 Note that in overall the ' 2 -norm error in the vertical axis is small because the displacements are of the order of centimeters; however, a difference of almost two orders of magnitude is obtained for some cases between the NIST and EnKI. This figure also shows that, in general, the NIST model using the fixed-base frequency and the signal predominant frequency yields smaller errors relative to 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. The discrepancies are attributed mostly to the poor estimation of the lumped rotational dashpot coefficients, which generate more oscillations in the transient part, and take longer to be damped out in the free oscillation part. These discrepancies are evident in Figure 13, where the NIST building responses are compared against the direct modeling approach. Note how the estimated NIST frequency (a) underestimates the inter-story drift ratio and the shear force at each floor, and (b) fails to capture the global period elongation and radiation damping for the displacement response. Figure 11. Impedance comparison between NIST using the flexible base periodT (in white circles), fixed base period T (in black solid dots), or the input signal dominant frequency f 0 (in asterisk) and EnKI using foundation input motion.

Dynamic global properties
Using the NIST and results of soil impedance obtained using EnKI-FIM, we computed the period elongation ratioT =T and foundation-damping b f that the ROM experiences when supported on lumped spring and dashpot elements (compared to fixed base conditions). Since the proposed ROM considers approximately the coupling between the translational and rotational degrees of freedoms, we employ Equations 31 and 37 derived in Appendix 1 to estimate the period lengthening and foundation damping. We then compare these results against the iterative method summarized in the ''Introduction'' section. Figures 14 and 15 show the period elongation and foundation damping computed using the two methods. Results show a good-agreement between the proposed expressions in the Appendix 1 and the NIST-proposed expressions, derived using the iterative method. From these figures, we note that Equation 31 results in a more flexible system compared to Equation 37, which yields a 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. 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 the reduced order modeling section, we are explicitly considering the spring and dashpot coupling terms in the stiffness and damping matrices. It can be seen in both Equations 31 and 37 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 global behavior of the two systems is very similar as shown in Figures 14 and 15. It is not a surprise that the period elongation ratios are well-captured given that the NIST recommendations compare well with the EnKI-FIM estimated lumped translational and rotational springs. Furthermore, the foundation damping for an equivalent one degree-of-freedom system is also well captured, which is expected  80,100,125,150,175,200,225,250,300,400, 500 m=s.  80, 100, 125, 150, 175, 200, 225, 250, 300, 400, 500 m=s. since the contribution of the lumped rotational dashpot to the global foundation damping is minimal. These differences become more pronounced for the stiffer buildings, whose response has a strong rocking component 3 that decreases the contribution of the foundation damping due to the coupling impedance term.

Frequency estimation of optimal EnKI-FIM impedance
Motivated by the inconsistencies that we demonstrated in the NIST definition of the most representative frequency to evaluate single valued stiffness and damping coefficients, we next used the EnKI algorithm to estimate the frequency at which one should read the numerically computed ''true'' impedance functions to determine the spring and dashpot values in the EnKI-FIM ROM. Again, we note that both the FEM and ROM analyses performed in this study are 2D plane-strain analyses. Therefore, all observed discrepancies are solely concerned with the extent to which the ROMs' impedance values are optimal. In this analysis, we consider two configurations: (a) only one frequency a Ã 0 is considered for the lumped soil impedances and (b) two frequencies, one for the horizontal (a x 0 ) and one for the rotational (a yy 0 ) lumped soil impedances. The results of this analysis are shown in Figure 16. Figure 16 reveals two interesting results: the first confirms the hypothesis that the dimensionless frequency that controls the interaction is related to the rotational soil impedance. This can be seen since the dimensionless frequency a Ã 0 in Figure 16a and a u 0 in Figure 16b are very similar. The second result suggests that the optimal a Ã 0 is not only a function of the building geometry and stiffness itself (T), but also a function of the foundation embedment ratio and the signal frequency content. In fact, Conti et al. (2020) studied the seismic performance of bridge piers and concluded that the overall damping and natural frequency changes significantly with the layout of caisson and pile foundation. Therefore, improving the functional form of the dimensionless frequency a 0 ;B=VsT could potentially be achieved by including more parameters, such as a 0 ;(h=VsT ), (h=B), (D=B), (Bf 0 =Vs), n s . However, our set of results did not reveal a clear pattern on how this dependency could be taken into account. Our results did reveal, however, that evaluating the flexible-base period using the NIST recommended iterative method may result in less accurate results, especially when the soil is stiffer. Future work by the authors will focus on learning a parameterized functional form that can incorporate the characteristics of the structure, foundation and shaking frequency in the iterative Figure 16. Optimal dimensionless frequency a Ã 0 obtained using EnKI for (left) one dimensionless frequency and (right) for two dimensionless frequency.
approach of evaluating the global dynamic parameters of building structures accounting for SSI effects.

Summary and conclusions
The state-of-practice introduces a major simplification to transform the frequency domain analysis into a time domain analysis. It assumes that the frequency-dependent soil impedance functions can be replaced by single-valued functions evaluated at the flexible-base structure frequency. This simplification has been adopted to enable time history analyses using standard finite element packages. In this study, we found that despite the simplification outlined above, the NIST recommendations yield global dynamic behavior of building structures that compares well to direct finite element solutions. However, their performance deteriorates in capturing the response of the structure such as the floor displacement timeseries, 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. The latter can be easily seen in Figure 9 for the lumped rotational component when a 0 ! 0 which in particular occurs for flexible buildings and stiff soils. In this regard, results in Figure 12 show that the fixed-base period in general provides better results since a 0 (T ) is greater than a 0 (T ) and thus the lumped rotational dashpot is better estimated; 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 the phase of the foundation input motion was also incorporated; thus, the discrepancies are only attributed to the frequency at which ''true'' impedance functions for plane strain setting are read.
Moreover, we showed that the foundation input motion effects are negligible for cases considered here with embedment ratios smaller than D=B<0:5. The latter is due to the fact that foundation embedment are indeed small and vertically propagating shear waves do not create base slab averaging effects. Therefore, the free field motion can be directly applied without introducing substantial errors. Finally, we noticed that the selection of the dimensionless frequency a 0 that should be used to evaluate the impedance values depends on other dimensionless parameters, such as the frequency of the input signal and the embedment ratio. This suggests that determination of the representative a 0 is not a sole function of period elongation. In fact, Conti et al. (2020) showed that the effect on the SSI response should be referred to both the structure and the foundation system. Future work by the authors is aiming to learn a functional form that account for these parameters to estimate the most representative frequency of dynamic building-foundation-soil systems.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/ or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/ or publication of this article: The authors would like to acknowledge CONICYT-Becas Chile for the financial support to the first author for conducting doctoral studies at The California Institute of Technology.
F k + F k yy À hk xy À Á k yy k x À k 2 xy + F k xy À hk x À Á k 2 xy À k x k yy h 2 = F 1 k + h 2 k x À 2 hk xy + k yy k x k yy À k 2 xy ! : Note that, for the fixed-base system to be equivalent to the flexible-base system, the displacements must be equal, that is,D = D. By comparison of Equations 23 and 24, we conclude that the following relation between the stiffness of the two systems must hold 1 k = 1 k + h 2 k x À 2 hk xy + k yy k x k yy À k 2 xy : ð25Þ Next, we now represent the impedance K 2 C of the system in complex format, this is where b j = vc j =2 k j and the index j represents the horizontal, vertical, or coupled degree of freedom to be considered. Employing Equation 25 and the impedance represented in complex format, one can write 1 k 1 + 2b f i À Á= 1 k 1 + 2b i i ð Þ + h 2 k x 1 + 2b x i ð ÞÀ2 hk xy 1 + 2b xy i À Á + k yy 1 + 2b yy i À Á k x k yy 1 + 2b x i ð Þ 1 + 2b yy i À Á À k 2 xy 1 + 2b xy i If each term in the previous expression is multiplied by its complex conjugated and assuming that the damping ratio for all degree of freedom are small so that (b j ) 2 '0, then Equation 27 becomes where the variablesk 2 R + andb 2 R + in Equation 28 are defined as followŝ