2.5-dimension soil seismic response to oblique incident waves based on exact free-field solution

Abstract

With the three dimensional (3D) oblique incident waves exactly determined for the free field, the soil seismic responses in both frequency and time domains are studied by the 2.5 dimension (2.5D) finite/infinite element method. First, the free-field responses in frequency domain are solved exactly for 3D arbitrary incident P and SV waves, which requires no coordinate conversion or extra effort for SV waves with super-critical incident angles. Next, the earthquake spectra are incorporated by the concept of equivalent seismic forces on the near-field boundary, based only on the displacements input derived for unit ground accelerations of each frequency using the 2.5D approach. For the asymmetric 2.5D finite/infinite element model adopted, the procedure for soil seismic analysis is presented. The solutions computed by the proposed method are verified against those of Wolf’s and de Barros and Luco’s and for inversely calculated ground motions. Of interest is that abrupt variation in soil response occurs around the critical angle on the wave propagation plane for SV waves. In addition, the horizontal displacements attenuate with increasing horizontal incident angle, while the longitudinal ones increase inversely for 3D incident P and SV waves.

Graphical abstract

Keywords

3D oblique incident wave / seismic response analysis / 2.5D approach / infinite element / half-space

1 Introduction

To accommodate the needs of urban development, there is a global trend toward the utilization of underground space to offer various functions, such as subways, stations, and commercial complex, etc. However, the damages posed by earthquakes to underground structures should not be overlooked due to their inherent nature in after-disaster rescue and restoration. Unfortunately, the analysis techniques available for underground structures seem not comparable to those for above-ground structures, such as bridges and buildings. Concerning the seismic analysis of underground structures, there still exist some fundamental problems to overcome. As part of the effort, this paper is aimed at developing a rational procedure for analyzing the dynamic response of a half-space to three-dimensional (3D) arbitrary incident seismic waves using the 2.5D method, focusing on accurate representation of the incident waves using the exact solution for the free field. The following literature review will start from the theoretical researches and then move to the numerical methods, including those based mainly on the boundary element method (BEM) and the finite elements. Finally, the focus is placed on the 2.5D finite/infinite element method (FIEM), which will be adopted in this paper for its versatility and effectiveness.

Previously, scholars have studied analytically the 3D dynamic responses of a site under incident waves considering the scattered waves of underground structures or local topography by various mathematical means. For local topography, the diffraction and scattering of elastic waves were studied under various 3D surface irregularities [1–3]. For underground structures, Wolf [4] formulated the problem of dynamic soil−structure interaction with arbitrary incident waves theoretically. As for the 3D seismic responses, underground structures, such as buried pipes, embedded foundations, and lined tunnels, were also studied [5–7]. Inspired by Wolf’s work [4], the exact dynamic stiffness matrices of the 3D layered half-space were adopted by Liang and Ba [8] to obtain the analytical solution of the free-field response more precisely. Recently, the dynamic behavior of the half-space with transversely isotropic and poroelastic properties was analyzed by Zhang and Pan [9] and Zhou et al. [10]. Focusing on the surface waves, closed-form solutions were derived for an unsupported shallow circular opening under Rayleigh waves by Bobet et al. [11].

Considering the complexity and limitation of analytical solutions, attempts have been made to develop various numerical methods for soil simulation and wave analysis. In this regard, the half-space was usually divided into two parts as the near field for the part concern with underground structures or local topography, and the far field for the remaining infinite domain with radiation damping. For underground seismic analysis, the BEM of various forms has been adopted, such as the indirect boundary integral equation method [12,13], indirect boundary method (IBM) [14–17], direct boundary element method [18–20], and indirect boundary element method [21], etc. Although methods based on the BEM have been adopted widely, a proper Green’s function is essential to the problem of concern in the formulation of boundary elements, which is not always readily available. In addition, other types of absorbing boundaries were also developed for simulating the infinite domain in numerical analysis, such as the viscous-spring artificial boundary [22,23], scaled boundary finite elements [24], viscous absorbing boundaries [25], improved 3D time-domain artificial boundary [26], the substructure of artificial boundaries [27], the fast multipole IBM [28], etc. In some recent works [29,30], the viscous-spring artificial boundary was used to study the effect of local topography for complex soil models or structures. Overall, each of the methods mentioned above has its limitations, since not all the boundary parameters are well determined. Hence, it is necessary to develop methods that can simulate the far field more accurately, considering the complicated underground incident waves.

Meanwhile, infinite elements have been proposed by Ungless [31] and Bettess [32] to simulate the far field, which can be easily implemented, since infinite elements are a degenerated form of finite elements. Zhao et al. [33–35] used the finite/infinite element method (FIEM) to study the wave scattering problems in infinite media with local topography. This method allows the variations in material and geometry of the half-space to be easily simulated. Yang et al. [36] presented a dynamic condensation scheme for generating the impedance matrices for all frequencies automatically, together with the key parameters (e.g., required mesh range, element size, attenuation coefficients, wavenumbers,) specified for the finite and infinite elements used. Later, infinite elements have been applied to various problems including the soil−structure interaction [37,38], seismic analysis of cavities for vertical P- and SV-waves [39], and of double tunnels [40]. Recently, Yang et al. [41] improved the FIEM for 2D arbitrary incident P and SV waves in the soil to deal with the problem of water-filled valleys [42]. However, there is still a lack of a rational procedure for treating the 3D oblique incident waves using the FIEM.

In practice, the 2D approach has often been adopted to analyze long underground structures using their cross-section profiles. However, this suffers from the drawback that the results obtained are generally qualitative and unable to account for the wave propagation effect along the tunnel axis. Though a full 3D modeling can offer all the waves messages for long underground structures, the computational cost involved is prohibitive.

As a compromise between the 2D and 3D approaches, Yang and Hung [43] proposed the 2.5D FIEM, which used only the 2D profile, for analyzing the 3D wave propagation behavior in the half-space under the moving loads, assuming the soil−structure system to be uniform along the load-moving direction. This method requires the plane-strain element to be added a third out-of-plane degree of freedom (DOF) at each node, in addition to the two in-plane DOFs. Subsequently, this approach was applied to studying various problems, e.g., layered medium vibrations [44], train-induced vibration isolation [45], seismic response of underground tunnels for vertical incident waves [46], vibration mitigation of railway-side buildings [47], stress wave propagation [48], internal loads excitation [49], coupling effect of the seismic and high-speed train loads [50]. Nowadays, the 2.5D concept has been extended to various wave propagation problems using methods not restricted to the FIEM [51–57]. But, there is still a lack of a proper 2.5D method to deal with the 3D oblique incidence of seismic waves.

To take advantage of both the infinite elements and the 2.5D method, it is necessary to reformulate the theoretical framework of the 2.5D FIEM so that it can accommodate 3D arbitrary incident seismic waves. While the theoretical framework adopted is similar to the previous ones [41,42], the present paper advances in two aspects: 1) the original 2D formulation is extended to the 2.5D for 3D soil seismic analysis; 2) the exact free-field solutions for inputting the 3D arbitrary incident seismic waves based only on the boundary displacements are presented for the first time.

The contents of the present paper are as follows. Section 2: Based on the wave equations given in Electronic Supplementary Materials, the exact solutions of the free-field responses in frequency-space domain are derived for 3D arbitrary incident P or SV waves. Section 3: For a unit acceleration set on the origin (for P or SV waves), the equivalent seismic forces are obtained from the exact displacement for the near-field boundary in frequency-wavenumber (k_{\mathit{\text{z}}}) domain, laying the foundation for treating each frequency of the earthquake spectra. Section 4: The asymmetric 2.5D finite/infinite element model and approach are briefed, with key parameters specified, along with the procedure of seismic analysis summarized. Section 5: Through verification against Wolf’s and de Barros and Luco’s solutions and inverse computation, the reliability of 2.5D FIEM for 3D oblique incident waves is ensured. Section 6: A parametric study is conducted in both frequency and time domains, focusing on the effect of the horizontal incident angles on the 3D half-space seismic response. Section 7: Conclusions.

2 Free field under 3D arbitrary incident seismic waves

A 3D half-space subjected to 3D arbitrary incident P or SV waves (see Fig.1), i.e., with seismic waves propagating both along the x- and z-axes, is considered for seismic analysis. The x–y plane is usually designated for the profile of a long structure (such as a tunnel) directed along the z-axis, which will be discretized in a later section. As the first step of the half-space seismic analysis, the closed-form solution of the free-field response to be used for the equivalent seismic forces calculation will be elaborated in this section.

2.1 Closed-form solution for P waves

Similar to the 2D problem for 2D incident waves [41], the 3D oblique incident P waves (P_{\mathrm{I}}) can also induce waveform conversion and wave reflection on the free surface, meaning that both reflected P and SV waves (P_{\mathrm{R}} and {SV}_{\mathrm{R}}) exist on the wave propagation plane shown in Fig.2(a). As was mentioned, the x–y plane is designated as the profile of a long structure (such as a tunnel) directed along the z-axis. Let {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}} denote the incident angle between the incident P waves and the x-axis on the x–z plane, and {\theta }_{\mathrm{P}\mathrm{i}\mathrm{v}} the incident angle between the incident P waves and the y-axis on the wave propagation plane.

Since the wave propagation plane (shown in Fig.2(a)) does not coincide with either of the two vertical planes (x–y and y–z planes) in the 3D Cartesian coordinates, the wave numbers k_{\mathrm{P}\mathrm{i},\mathrm{h}},k_{\mathrm{P}\mathrm{P}\mathrm{r},\mathrm{h}},k_{\mathrm{P}\mathrm{S}\mathrm{r},\mathrm{h}} (see Eq. (A9) in Electronic Supplementary Materials) need to be decomposed along the x- and z-axes with respect to incident angle on the x–z plane {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}}. Therefore, one can also infer that all the incident and reflected waves have the same wave numbers along both the x- and z-axes for P waves. In this regard, two symbols, i.e., k_{x\mathrm{P}} and k_{\mathit{\text{z}}\mathrm{P}}, are used to denote the wave numbers along the x- and z-axes, respectively, for P waves. Then, the wave numbers k_{x\mathrm{P}} and k_{\mathit{\text{z}}\mathrm{P}} can be expressed as:

\begin{array}{}\text{(1a)}& k_{x\mathrm{P}}=k_{\mathrm{P}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{P}\mathrm{i}\mathrm{v}}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}},\end{array}

\begin{array}{}\text{(1b)}& k_{\mathit{\text{z}}\mathrm{P}}=k_{\mathrm{P}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{P}\mathrm{i}\mathrm{v}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}}.\end{array}

By letting r_{\mathrm{P}}^2=k_{x\mathrm{P}}^2+k_{\mathit{\text{z}}\mathrm{P}}^2-k_{\mathrm{P}}^2, r_{\mathrm{S}}^2=k_{x\mathrm{P}}^2+k_{\mathit{\text{z}}\mathrm{P}}^2-k_{\mathrm{S}}^2, the solutions of the transformed potential functions to the differential equations in Eq. (A8) in Electronic Supplementary Materials can be given as:

\begin{array}{}\text{(2a)}& {\hat{\hat{\Phi }}}_{\mathrm{P}}=I_{\mathrm{P}}\mathrm{e}\mathrm{x}\mathrm{p}\left(r_{\mathrm{P}}y\right)+R_{\mathrm{P}}\mathrm{e}\mathrm{x}\mathrm{p}\left({-r}_{\mathrm{P}}y\right),\end{array}

\begin{array}{}\text{(2b)}& {\hat{\hat{\Psi }}}_{xy,\mathrm{P}}^{\mathrm{S}\mathrm{V}}=R_{xy,\mathrm{S}}\mathrm{e}\mathrm{x}\mathrm{p}\left({-r}_{\mathrm{S}}y\right),\end{array}

\begin{array}{}\text{(2c)}& {\hat{\hat{\Psi }}}_{y\mathit{\text{z}},\mathrm{P}}^{\mathrm{S}\mathrm{V}}=R_{y\mathit{\text{z}},\mathrm{S}}\mathrm{e}\mathrm{x}\mathrm{p}\left({-r}_{\mathrm{S}}y\right),\end{array}

where I_{\mathrm{P}} and R_{\mathrm{P}} represent the amplitudes of the incident and reflected P waves, respectively, while R_{xy,\mathrm{S}} and R_{y\mathit{\text{z}},\mathrm{S}} are the amplitudes of the reflected SV waves on the x–y and y–z planes, respectively.

Substituting Eq. (2) into Eqs. (A5) and (A6) in Electronic Supplementary Materials yields the transformed displacements and stresses related to I_{\mathrm{P}}, R_{\mathrm{P}}, R_{xy,\mathrm{S}}, and R_{y\mathit{\text{z}},\mathrm{S}} as follows:

\begin{array}{}\text{(3a)}& {\left\{\hat{\hat{u}},\hat{\hat{v}},\hat{\hat{w}}\right\}}_{\mathrm{P}}^{\mathrm{T}}=\left[{\mathrm{\Gamma }}_{\mathrm{P}}\right]\left[E_{\mathrm{P}}\right]{\left\{I_{\mathrm{P}},{R_{\mathrm{P}},R}_{xy,\mathrm{S}},R_{y\mathit{\text{z}},\mathrm{S}}\right\}}^{\mathrm{T}},\end{array}

\begin{array}{}\text{(3b)}& {\left\{{\hat{\hat{\sigma }}}_{xx},{\hat{\hat{\sigma }}}_{yy},{\hat{\hat{\sigma }}}_{\mathit{\text{zz}}},{\hat{\hat{\tau }}}_{xy},{\hat{\hat{\tau }}}_{y\mathit{\text{z}}},{\hat{\hat{\tau }}}_{x\mathit{\text{z}}}\right\}}_{\mathrm{P}}^{\mathrm{T}}=\left[{\mathrm{\Lambda }}_{\mathrm{P}}\right]\left[E_{\mathrm{P}}\right]{\left\{I_{\mathrm{P}},{R_{\mathrm{P}},R}_{xy,\mathrm{S}},R_{y\mathit{\text{z}},\mathrm{S}}\right\}}^{\mathrm{T}},\end{array}

where \left[{\mathrm{\Gamma }}_{\mathrm{P}}\right], \left[{\mathrm{\Lambda }}_{\mathrm{P}}\right], and \left[E_{\mathrm{P}}\right] are the coefficient matrices that can be obtained by substituting Eq. (2) into Eqs. (A5) and (A6), which are not shown here for brevity. For the 3D problem considered, the stress boundary conditions of the free field are:

\begin{array}{}\text{(4a)}& {\hat{\hat{\sigma }}}_{yy}\left(y=0\right)=0,\end{array}

\begin{array}{}\text{(4b)}& {\hat{\hat{\tau }}}_{xy}\left(y=0\right)=0,\end{array}

\begin{array}{}\text{(4c)}& {\hat{\hat{\tau }}}_{y\mathit{\text{z}}}\left(y=0\right)=0.\end{array}

By combining Eq. (4) and Eq. (3b), and using the wave velocity ratio V=c_{\mathrm{P}}/c_{\mathrm{S}}=\sqrt{\left(\lambda +2\mu \right)/\mu }, the three amplitude ratios for 3D arbitrary incident P waves can be obtained as:

\begin{array}{}\text{(5a)}& f_{\mathrm{P}\mathrm{P}}=\frac{R_{\mathrm{P}}}{I_{\mathrm{P}}}=1-\frac{2K_{\mathrm{P}}{Q}_{\mathrm{P}}}{K_{\mathrm{P}}{Q}_{\mathrm{P}}-4r_{\mathrm{P}}{r}_{\mathrm{S}}\left(k_{x\mathrm{P}}^2+k_{\mathit{\text{z}}\mathrm{P}}^2\right)},\end{array}

\begin{array}{}\text{(5b)}& f_{\mathrm{P}\mathrm{S},xy}=\frac{R_{xy,\mathrm{S}}}{I_{\mathrm{P}}}=\frac{4\mathrm{i}{k}_{x\mathrm{P}}{r}_{\mathrm{P}}{Q}_{\mathrm{P}}}{K_{\mathrm{P}}{Q}_{\mathrm{P}}-4r_{\mathrm{P}}{r}_{\mathrm{S}}\left(k_{x\mathrm{P}}^2+k_{\mathit{\text{z}}\mathrm{P}}^2\right)},\end{array}

\begin{array}{}\text{(5c)}& f_{\mathrm{P}\mathrm{S},y\mathit{\text{z}}}=\frac{R_{y\mathit{\text{z}},S}}{I_{\mathrm{P}}}=\frac{4\mathrm{i}{k}_{\mathit{\text{z}}\mathrm{P}}{r}_{\mathrm{P}}{Q}_{\mathrm{P}}}{K_{\mathrm{P}}{Q}_{\mathrm{P}}-4r_{\mathrm{P}}{r}_{\mathrm{S}}\left(k_{x\mathrm{P}}^2+k_{\mathit{\text{z}}\mathrm{P}}^2\right)},\end{array}

where f_{\mathrm{P}\mathrm{P}}, f_{\mathrm{P}\mathrm{S},xy}, and f_{\mathrm{P}\mathrm{S},y\mathit{\text{z}}} denote the amplitude ratios of the reflected P and SV waves on the x–y and y–z planes to incident P waves, respectively, with K_{\mathrm{P}}=\left[2\left(k_{x\mathrm{P}}^2+k_{\mathit{\text{z}}\mathrm{P}}^2\right)-k_{\mathrm{S}}^2\right] and Q_{\mathrm{P}}=\left[2\left(k_{x\mathrm{P}}^2+k_{\mathit{\text{z}}\mathrm{P}}^2\right)-V^2{k}_{\mathrm{P}}^2\right].

Lastly, by substituting the Eq. (5) into Eq. (3) and employing the inverse Fourier transform with respect to k_{x\mathrm{P}} and k_{\mathit{\text{z}}\mathrm{P}}, one can express the frequency-space-domain displacements and stresses of the free field (in single Fourier transform) for P waves with fixed horizontal and vertical incident angles, i.e., {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}} and {\theta }_{\mathrm{P}\mathrm{i}\mathrm{v}}, as follows:

\begin{array}{rl}{\left\{\tilde{u},\tilde{v},\tilde{w}\right\}}_{\mathrm{P}}^{\mathrm{T}}=& \left[{\mathrm{\Gamma }}_{\mathrm{P}}\right]\left[E_{\mathrm{P}}\right]{\left\{1,f_{\mathrm{P}\mathrm{P}},f_{\mathrm{P}\mathrm{S},xy},f_{\mathrm{P}\mathrm{S},y\mathit{\text{z}}}\right\}}^{\mathrm{T}}\\ & \cdot I_{\mathrm{P}}\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{i}{k}_{x\mathrm{P}}x)\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{i}{k}_{\mathit{\text{z}}\mathrm{P}}\mathit{\text{z}})\end{array}

\begin{array}{}\text{(6a)}& {\left\{{\tilde{\sigma }}_{xx},{\tilde{\sigma }}_{yy},{\tilde{\sigma }}_{\mathit{\text{zz}}},{\tilde{\tau }}_{xy},{\tilde{\tau }}_{y\mathit{\text{z}}},{\tilde{\tau }}_{x\mathit{\text{z}}}\right\}}_{\mathrm{P}}^{\mathrm{T}},\end{array}

\begin{array}{}\text{(6b)}& =\left[{\mathrm{\Lambda }}_{\mathrm{P}}\right]\left[E_{\mathrm{P}}\right]{\left\{1,f_{\mathrm{P}\mathrm{P}},f_{\mathrm{P}\mathrm{S},xy},f_{\mathrm{P}\mathrm{S},yz}\right\}}^{\mathrm{T}}{I}_{\mathrm{P}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{i}{k}_{x\mathrm{P}}x\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{i}{k}_{\mathit{\text{z}}\mathrm{P}}\mathit{\text{z}}\right),\end{array}

where a single Fourier transform is denoted with “~”. The preceding expressions in Eqs. (6a) and (6b) are called the closed-form or exact solutions for the frequency-space-domain displacements and stresses, respectively, of the entire free field under 3D arbitrary incident P waves.

2.2 Closed-form solution for SV waves

Similar to the case for P waves, the 3D wave propagation behavior in the free field for SV waves was depicted in Fig.2(b), in which {\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}} denotes the incident angle between the incident SV waves and the x-axis on the x–z plane, and {\theta }_{\mathrm{S}\mathrm{i}\mathrm{v}} the incident angle between the incident SV waves and the y-axis on the wave propagation plane.

Then, by projecting all the wave numbers k_{\mathrm{S}\mathrm{i},\mathrm{h}},k_{\mathrm{S}\mathrm{P}\mathrm{r},\mathrm{h}}, k_{\mathrm{S}\mathrm{S}\mathrm{r},\mathrm{h}} (Eq. (A10) in Electronic Supplementary Materials) on the x- and z-axes with respect to the incident angle {\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}} on the x–z plane, it can be also concluded that the wave numbers along both the x- and z-axes (i.e., k_{x\mathrm{S}} and k_{\mathit{\text{z}}\mathrm{S}}, respectively) are the same for all the SV waves. They can be expressed as:

\begin{array}{}\text{(7a)}& k_{x\mathrm{S}}=k_{\mathrm{S}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{S}\mathrm{i}\mathrm{v}}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}},\end{array}

\begin{array}{}\text{(7b)}& k_{\mathit{\text{z}}\mathrm{S}}=k_{\mathrm{S}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{S}\mathrm{i}\mathrm{v}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}}.\end{array}

By letting r_{\mathrm{P}}^2=k_{x\mathrm{S}}^2+k_{\mathit{\text{z}}\mathrm{S}}^2-k_{\mathrm{P}}^2, r_{\mathrm{S}}^2=k_{x\mathrm{S}}^2+k_{\mathit{\text{z}}\mathrm{S}}^2-k_{\mathrm{S}}^2, the solutions to the differential equations in transformed form in Eq. (A8) in Electronic Supplementary Materials can be obtained as:

\begin{array}{}\text{(8a)}& {\hat{\hat{\Phi }}}_{\mathrm{S}}=R_{\mathrm{P}}\mathrm{e}\mathrm{x}\mathrm{p}\left({-r}_{\mathrm{P}}y\right),\end{array}

\begin{array}{}\text{(8b)}& {\hat{\hat{\Psi }}}_{xy,\mathrm{S}}^{\mathrm{S}\mathrm{V}}=I_{\mathrm{S},xy}\mathrm{e}\mathrm{x}\mathrm{p}\left(r_{\mathrm{S}}y\right)+R_{\mathrm{S},xy}\mathrm{e}\mathrm{x}\mathrm{p}\left({-r}_{\mathrm{S}}y\right),\end{array}

\begin{array}{}\text{(8c)}& {\hat{\hat{\Psi }}}_{y\mathit{\text{z}},\mathrm{S}}^{\mathrm{S}\mathrm{V}}=I_{\mathrm{S},y\mathit{\text{z}}}\mathrm{e}\mathrm{x}\mathrm{p}\left(r_{\mathrm{S}}y\right)+R_{\mathrm{S},y\mathit{\text{z}}}\mathrm{e}\mathrm{x}\mathrm{p}\left({-r}_{\mathrm{S}}y\right),\end{array}

where I_{\mathrm{S},xy} and R_{\mathrm{S},xy} represent the amplitudes of incident and reflected SV waves on the x–y plane, respectively, while I_{\mathrm{S},y\mathit{\text{z}}} and R_{\mathrm{S},y\mathit{\text{z}}} on the z–y plane, and R_{\mathrm{P}} is the amplitude of reflected P waves.

Substituting Eq. (8) into Eqs. (A5) and (A6) in Electronic Supplementary Materials yields the transformed displacements and stresses with respect to R_{\mathrm{P}}, I_{\mathrm{S},xy}, I_{\mathrm{S},y\mathit{\text{z}}}, R_{\mathrm{S},xy}, and R_{\mathrm{S},y\mathit{\text{z}}} as:

\begin{array}{}\text{(9a)}& {\left\{\hat{\hat{u}},\hat{\hat{v}},\hat{\hat{w}}\right\}}_{\mathrm{S}}^{\mathrm{T}}=\left[{\mathrm{\Gamma }}_{\mathrm{S}}\right]\left[E_{\mathrm{S}}\right]{\left\{{R_{\mathrm{P}},I_{\mathrm{S},xy},I_{\mathrm{S},y\mathit{\text{z}}},R}_{\mathrm{S},xy},R_{\mathrm{S},y\mathit{\text{z}}}\right\}}^{\mathrm{T}},\end{array}

\begin{array}{rl}& {\left\{{\hat{\hat{\sigma }}}_{xx},{\hat{\hat{\sigma }}}_{yy},{\hat{\hat{\sigma }}}_{\mathit{\text{zz}}},{\hat{\hat{\tau }}}_{xy},{\hat{\hat{\tau }}}_{y\mathit{\text{z}}},{\hat{\hat{\tau }}}_{x\mathit{\text{z}}}\right\}}_{\mathrm{S}}^{\mathrm{T}}\\ =& \left[{\mathrm{\Lambda }}_{\mathrm{S}}\right]\left[E_{\mathrm{S}}\right]{\left\{{R_{\mathrm{P}},I_{\mathrm{S},xy},I_{\mathrm{S},y\mathit{\text{z}}},R}_{\mathrm{S},xy},R_{\mathrm{S},y\mathit{\text{z}}}\right\}}^{\mathrm{T}},\end{array}

where \left[{\mathrm{\Gamma }}_{\mathrm{S}}\right], \left[{\mathrm{\Lambda }}_{\mathrm{S}}\right], and \left[E_{\mathrm{S}}\right] are the coefficient matrices for SV waves that can be easily obtained by substituting Eq. (8) into Eqs. (A5) and (A6).

Here, for unified operation, the amplitude ratios are specified as those of the waves components \left(R_{\mathrm{P}},I_{\mathrm{S},xy},I_{\mathrm{S},y\mathit{\text{z}}},R_{\mathrm{S},xy},R_{\mathrm{S},y\mathit{\text{z}}}\right) to the total incident SV waves (I_{\mathrm{S}}) on the wave propagation plane shown in Fig.2(b). First, the relationships of I_{\mathrm{S},xy}, I_{\mathrm{S},y\mathit{\text{z}}}, and I_{\mathrm{S}} can be given as:

\begin{array}{}\text{(10a)}& f_{\mathrm{S},xy}=\frac{I_{\mathrm{S},xy}}{I_{\mathrm{S}}}=\frac{I_{\mathrm{S},xy}}{\sqrt{I_{\mathrm{S},xy}^2+I_{\mathrm{S},y\mathit{\text{z}}}^2}}=\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}},\end{array}

\begin{array}{}\text{(10b)}& f_{\mathrm{S},y\mathit{\text{z}}}=\frac{I_{\mathrm{S},y\mathit{\text{z}}}}{I_{\mathrm{S}}}=\frac{I_{\mathrm{S},xy}}{\sqrt{I_{\mathrm{S},xy}^2+I_{\mathrm{S},y\mathit{\text{z}}}^2}}=\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}}.\end{array}

Then, by substituting the stress boundary conditions in Eq. (4) into Eq. (9b), one obtains the other amplitude ratios under 3D oblique incident SV waves as:

\begin{array}{}\text{(11a)}& f_{\mathrm{S}\mathrm{P}}=\frac{R_{\mathrm{P}}}{I_{\mathrm{S}}}=\frac{-4\mathrm{i}{r}_{\mathrm{S}}{K}_{\mathrm{S}}\left(k_{x\mathrm{S}}^2+k_{\mathit{\text{z}}\mathrm{S}}^2\right)}{\left[K_{\mathrm{S}}{Q}_{\mathrm{S}}-4r_{\mathrm{P}}{r}_{\mathrm{S}}\left(k_{x\mathrm{S}}^2+k_{\mathit{\text{z}}\mathrm{S}}^2\right)\right]k_{\mathrm{S}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{S}\mathrm{i}\mathrm{v}}},\end{array}

\begin{array}{}\text{(11b)}& f_{\mathrm{S}\mathrm{S},xy}=\frac{R_{x\mathrm{S}}}{I_{\mathrm{S}}}=\left[\frac{-8r_{\mathrm{S}}{r}_{\mathrm{P}}\left(k_{x\mathrm{S}}^2+k_{\mathit{\text{z}}\mathrm{S}}^2\right)}{\left[K_{\mathrm{S}}{Q}_{\mathrm{S}}-4r_{\mathrm{P}}{r}_{\mathrm{S}}\left(k_{x\mathrm{S}}^2+k_{\mathit{\text{z}}\mathrm{S}}^2\right)\right]}-1\right]\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}},\end{array}

\begin{array}{}\text{(11c)}& f_{\mathrm{S}\mathrm{S},y\mathit{\text{z}}}=\frac{R_{\mathit{\text{z}}\mathrm{S}}}{I_{\mathrm{S}}}=\left[\frac{-8r_{\mathrm{S}}{r}_{\mathrm{P}}\left(k_{x\mathrm{S}}^2+k_{\mathit{\text{z}}\mathrm{S}}^2\right)}{\left[K_{\mathrm{S}}{Q}_{\mathrm{S}}-4r_{\mathrm{P}}{r}_{\mathrm{S}}\left(k_{x\mathrm{S}}^2+k_{\mathit{\text{z}}\mathrm{S}}^2\right)\right]}-1\right]\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}},\end{array}

where f_{\mathrm{S}\mathrm{P}}, f_{\mathrm{S}\mathrm{S},xy}, and f_{\mathrm{S}\mathrm{S},y\mathit{\text{z}}} denote the amplitude ratios of the reflected P and SV waves on the x–y and z–y planes to the total incident SV waves, respectively, and Q_{\mathrm{S}}=\left[2\left(k_{x\mathrm{S}}^2+k_{\mathit{\text{z}}\mathrm{S}}^2\right)-V^2{k}_{\mathrm{P}}^2\right] and K_{\mathrm{S}}=\left[2\left(k_{x\mathrm{S}}^2+k_{\mathit{\text{z}}\mathrm{S}}^2\right)-k_{\mathrm{S}}^2\right].

Finally, substituting Eqs. (10) and (11) into Eq. (9) and performing the inverse Fourier transform with respect to k_{x\mathrm{S}} and k_{\mathit{\text{z}}\mathrm{S}}, one can obtain the frequency-space-domain displacements and stresses of the free field (in single Fourier transform) as follows:

\begin{array}{rl}{\left\{\tilde{u},\tilde{v},\tilde{w}\right\}}_{\mathrm{S}}^{\mathrm{T}}=& \left[{\mathrm{\Gamma }}_{\mathrm{S}}\right]\left[E_{\mathrm{S}}\right]{\left\{f_{\mathrm{S}\mathrm{P}},f_{\mathrm{S},xy},f_{\mathrm{S},yz},f_{\mathrm{S}\mathrm{S},xy},f_{\mathrm{S}\mathrm{S},yz}\right\}}^{\mathrm{T}}\\ & \cdot I_{\mathrm{S}}\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{i}{k}_{x\mathrm{S}}x)\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{i}{k}_{\mathit{\text{z}}\mathrm{S}}\mathit{\text{z}}),\end{array}

\begin{array}{rl}& {\left\{{\tilde{\sigma }}_{xx},{\tilde{\sigma }}_{yy},{\tilde{\sigma }}_{zz},{\tilde{\tau }}_{xy},{\tilde{\tau }}_{y\mathit{\text{z}}},{\tilde{\tau }}_{x\mathit{\text{z}}}\right\}}_{\mathrm{S}}^{\mathrm{T}}\\ =& \left[{\mathrm{\Lambda }}_{\mathrm{S}}\right]\left[E_{\mathrm{S}}\right]{\left\{f_{\mathrm{S}\mathrm{P}},f_{\mathrm{S},xy},f_{\mathrm{S},y\mathit{\text{z}}},f_{\mathrm{S}\mathrm{S},xy},f_{\mathrm{S}\mathrm{S},y\mathit{\text{z}}}\right\}}^{\mathrm{T}}\\ & \cdot I_{\mathrm{S}}\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{i}{k}_{x\mathrm{S}}x)\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{i}{k}_{\mathit{\text{z}}\mathrm{S}}\mathit{\text{z}}).\end{array}

Similarly, the preceding expressions in Eqs. (12a) and (12b), are the frequency-space-domain closed-form or exact solutions for the displacements and stresses, respectively, of the entire free field under 3D arbitrary incident SV waves.

Besides, for the 2D problems [41], one knows that for oblique incident SV waves with angle {\theta }_{\mathrm{S}\mathrm{i}\mathrm{v}} on the wave propagation plane, a critical angle {\theta }_{\mathrm{s}\mathrm{i}\mathrm{v}}^{\mathrm{c}\mathrm{r}} exists such that the reflected P waves disappear and transform into the surface waves. Moreover, the closed-form solutions for the response calculated from the transformed potential function {\mathrm{\Psi }}_{xy,\mathrm{S}}^{\mathrm{S}\mathrm{V}} in terms of \mathrm{e}\mathrm{x}\mathrm{p}\left({-r}_{\mathrm{S}}y\right) in frequency domain have been proved to be valid for SV waves with 2D sub- and super-critical incident angles [41]. Since the critical angle {\theta }_{\mathrm{s}\mathrm{i}\mathrm{v}}^{\mathrm{c}\mathrm{r}} is independent of the incident angle {\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}} on the x–z plane (see Fig.2(b)), it follows that the displacements and stresses derived in Eq. (12) from the potential functions {\mathrm{\Psi }}_{xy,\mathrm{S}}^{\mathrm{S}\mathrm{V}} and {\mathrm{\Psi }}_{y\mathit{\text{z}},\mathrm{S}}^{\mathrm{S}\mathrm{V}} (in terms of \mathrm{e}\mathrm{x}\mathrm{p}\left({-r}_{\mathrm{S}}y\right) in Eq. (8)) remain equally valid for 3D arbitrary incident SV waves with no limitation on the critical angle.

The incident-wave amplitudes \left(I_{\mathrm{P}},I_{\mathrm{S}}\right) contained in Eqs. (6) and (12) will be determined by additional boundary conditions for the ground motion input in Subsection 3.1. One advantage with the above procedure for solving the 3D free-field responses under the 3D arbitrary incident P and SV waves is that it can directly reveal the wave propagation behavior in the 3D space. Moreover, it differs from the existing ones, e.g., Zhang and Chopra [18], in that no coordinate conversion is required, as the solutions have been directly given in the 3D global coordinates, as in Eqs. (6) and (12). Given the closed-form solutions of the free field, the equivalent seismic forces that are necessary to the frequency-based finite/infinite element analysis will be formulated in the following section. Besides, the above process for dealing with the homogenous and isotropic half-space can be extended to the layered soils once the continuity conditions between soil layers are taken into account, which is part of the future work.

3 Equivalent seismic forces calculation for 3D arbitrary incident seismic waves

With reference to Fig.3 for the numerical model, the half-space on the x–y plane in the 3D Cartesian coordinates is divided into two regions as the near field that contains the underground structures or local topography and part of the surrounding soil of interest, and the remaining far field that extends to infinity. The near field is simulated by finite elements and the far field by infinite elements. To ensure the accuracy of the solution, the criteria for selecting the mesh range of the near field were given by Yang et al. for the moving train loads [36] and by Lin et al. [40] for seismic analysis, which will be elaborated in the section to follow. By assuming the soil and embedded structure of interest to be uniform along the z-axis, only the 2D profile of the soil and embedded structure on the x–y plane of the 3D soil model needs to considered in simulation, which has been referred to as the 2.5D approach [43]. In this study, the near field will be discretized into finite elements and the far field into infinite elements, as will be elaborated in Subsection 4.1.

For 2D problems [41], it was known that both the stresses and displacements on the boundaries \left(B_{\mathrm{L}},B_{\mathrm{R}},B_{\mathrm{B}}\subset B\right) of the near field (see Fig.3) enter into the equivalent seismic force calculation, according to Zhao and Valliappan [58]. However, for 3D problems, it would be difficult to include all the six stress components \left({\sigma }_{xx}{,\sigma }_{yy},{\sigma }_{\mathit{\text{zz}}},{\tau }_{xy},{\tau }_{y\mathit{\text{z}}},{\tau }_{x\mathit{\text{z}}}\right) on the boundary of the near field to compose the equivalent seismic forces. Instead, an alternative method that utilizes only the displacement responses of the single layer elements enclosing the near-field boundary B was presented by Liu et al. [27]. Based on the free-field solution presented in Section 2, the single layer elements adjacent to the near-field boundary will be expanded for the 2.5D formulation for calculating the equivalent seismic forces in frequency domain, given the ground seismic records in Subsection 3.2.

3.1 Free-field displacement responses for equivalent seismic forces calculation used

For earthquake input, the nodal responses of the single layer (see Subsection 3.2) adjacent to the near-field boundary to a unit ground acceleration will first be generated using the theory developed for oblique incident waves in Section 2, prior to computation of the equivalent seismic forces. The boundary conditions of a unit horizontal or vertical acceleration for SV or P waves, respectively, at the origin on the x–y plane (Fig.3) for each frequency \omega can be specified as follows:

\begin{array}{}\text{(13a)}& -{\omega }^2\hat{\hat{u}}\left(y=0\right)=1,\end{array}

\begin{array}{}\text{(13b)}& -{\omega }^2\hat{\hat{v}}\left(y=0\right)=1.\end{array}

They can be converted into the displacement boundary conditions in frequency-wavenumbers (k_x and k_{\mathit{\text{z}}}) domain as follows:

\begin{array}{}\text{(14a)}& \hat{\hat{u}}\left(y=0\right)=-1/{\omega }^2,\end{array}

\begin{array}{}\text{(14b)}& \hat{\hat{v}}\left(y=0\right)=-1/{\omega }^2.\end{array}

Combining Eq. (14) with the free-field displacements in Eqs. (3a) and (9a) for the P and SV waves, respectively, one can first derive the two amplitudes \left(I_{\mathrm{P}},I_{\mathrm{S}}\right) in frequency and wavenumbers (k_x and k_{\mathit{\text{z}}}) domain as:

\begin{array}{}\text{(15a)}& I_{\mathrm{P}}={-\left[{\omega }^2\left({r_{\mathrm{P}}-r}_{\mathrm{P}}{f}_{\mathrm{P}\mathrm{P}}+\mathrm{i}{k}_{x\mathrm{P}}{f}_{\mathrm{P}\mathrm{S},xy}+\mathrm{i}{k}_{\mathit{\text{z}}\mathrm{P}}{f}_{\mathrm{P}\mathrm{S},y\mathit{\text{z}}}\right)\right]}^{-1},\end{array}

\begin{array}{}\text{(15b)}& I_{\mathrm{S}}={-\left[{\omega }^2\left(-\mathrm{i}{k}_{x\mathrm{S}}{f}_{\mathrm{S}\mathrm{P}}+r_{\mathrm{S}}{f}_{\mathrm{S},xy}-r_{\mathrm{S}}{f}_{\mathrm{S}\mathrm{S},xy}\right)\right]}^{-1}.\end{array}

For the 2D problems [41], since the responses calculated from the frequency-based finite/infinite element model have been given in the space \left(x,y\right) domain, the equivalent seismic forces need also be specified in the same frequency-space \left(x,y\right) domain. But for the 2.5D FIEM to be presented in Subsection 4.1, the responses to be computed will be given in the wavenumber (k_{\mathit{\text{z}}}) domain on the x–y plane. Therefore, the equivalent seismic forces should also be specified in the frequency-wavenumber (k_{\mathit{\text{z}}}) domain on the near-field boundary for given P and SV waves in the free field. Accordingly, one needs to remove the transformed component containing the wave number (k_{\mathit{\text{z}}}) in Eqs. (6a) and (12a), for the P and SV waves, respectively. Accordingly, one can obtain the free-field displacements {\hat{\mathit{U}}}^{\mathrm{f}} (denoted with “^”, and superscript f for the free-field response) in frequency-wavenumber (k_{\mathit{\text{z}}}) domain for P and SV waves, respectively, as follows:

\begin{array}{}\text{(16a)}& {\left\{\hat{u},\hat{v},\hat{w}\right\}}_{\mathrm{P}}^{\mathrm{T}}=\left[{\Gamma }_{\mathrm{P}}\right]\left[E_{\mathrm{P}}\right]{\left\{1,f_{\mathrm{P}\mathrm{P}},f_{\mathrm{P}\mathrm{S},xy},f_{\mathrm{P}\mathrm{S},y\mathit{\text{z}}}\right\}}^{\mathrm{T}}{I}_{\mathrm{P}}\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{i}{k}_{x\mathrm{P}}x),\end{array}

\begin{array}{}\text{(16b)}& {\left\{\hat{u},\hat{v},\hat{w}\right\}}_{\mathrm{S}}^{\mathrm{T}}=\left[{\Gamma }_{\mathrm{S}}\right]\left[E_{\mathrm{S}}\right]{\left\{f_{\mathrm{S}\mathrm{P}},f_{\mathrm{S},xy},f_{\mathrm{S},y\mathit{\text{z}}},f_{\mathrm{S}\mathrm{S},xy},f_{\mathrm{S}\mathrm{S},y\mathit{\text{z}}}\right\}}^{\mathrm{T}}{I}_{\mathrm{S}}\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{i}{k}_{x\mathrm{S}}x).\end{array}

Then, substituting the incident-wave amplitudes calculated in Eqs. (15a) and (15b) into Eqs. (16a) and (16b) for P and SV waves, respectively, one can obtain the frequency-wavenumber (k_{\mathit{\text{z}}}) domain displacements {\hat{\mathit{U}}}_{\mathrm{s}\mathrm{l}}^{\mathrm{f}} ({\hat{\mathit{U}}}_{\mathrm{s}\mathrm{l}}^{\mathrm{f}}\subset {\hat{\mathit{U}}}^{\mathrm{f}}, and subscript sl for the DOFs of the single layer) of each point \left(x,y\right) on the single layer, to be given in Subsection 3.2, for each unit ground acceleration of each frequency \omega. Such a result can be looped over each frequency within the range of the earthquake spectrum considered, to calculate the equivalent seismic forces. The single layer will be further explained below.

3.2 Equivalent seismic forces on near-field boundary

To calculate the equivalent seismic forces for the 2.5D method, it is required that the displacements {\hat{\mathit{U}}}_{\mathrm{s}\mathrm{l}}^{\mathrm{f}} of the single layer (Eq. (16)) be made available for a given seismic record in the free field. Here, the earthquake data adopted in Fig.4 are the acceleration history of the 1999 Chi-Chi Earthquake recorded with a sampling interval of \mathrm{\Delta }t=0.005\mathrm{s} at TCU068 station in Taichung. The horizontal acceleration in Fig.4(a) is used for SV waves, and the vertical acceleration in Fig.4(b) for P waves. Moreover, by using the fast Fourier transform (FFT), the acceleration frequency spectra for the two directions were obtained and plotted in Fig.5, where the horizontal axis denotes the circular frequency in Hz ( = 2 rad/s).

Based on Nyquist’s theorem, the frequency range for sampling is selected as 0 to 1/\left(2\mathrm{\Delta }t\right), which is equal to 0–100 Hz in this paper. To ensure proper resolution in the frequency domain, the frequency increment adopted in the FFT for generating the spectrum is \mathrm{\Delta }\omega =\text{π}/[\mathrm{\Delta }t(N_{\mathrm{h}\mathrm{a}\mathrm{l}\mathrm{f}}-1)] [41] or approximately 0.069 rad/s, where N_{\mathrm{h}\mathrm{a}\mathrm{l}\mathrm{f}} is half of the total number of frequency increments of the FFT, being set equal to the number of sampling points (Fig.4). Based on the spectra for the Chi-Chi Earthquake in Fig.5, the actual frequency range included in seismic analysis, including calculation of the equivalent seismic forces, is 0–25 Hz \left({\omega }_{\mathrm{a}\mathrm{c}\mathrm{t}}\right).

On the other hand, in the FFT and inverse FFT, the same frequency increment \mathrm{\Delta }\omega should be used. Also, as indicated by the earthquake spectra in Fig.5, the range up to 25 Hz is good enough, which is a quarter of the valid range of 0–100 Hz, meaning that the ratio of the valid range 1/\left(2\mathrm{\Delta }t\right) to the actual range {\omega }_{\mathrm{a}\mathrm{c}\mathrm{t}} used is \gamma =4. It follows that the actual time increment \mathrm{\Delta }{t}_{\mathrm{a}\mathrm{c}\mathrm{t}} for the time-history calculation by the IFFT can be set as \mathrm{\Delta }{t}_{\mathrm{a}\mathrm{c}\mathrm{t}}=\gamma \mathrm{\Delta }t. Since the increment of \mathrm{\Delta }t= 0.005 s has been used by the earthquake data in Fig.4, the time increment adopted in the IFFT calculation is simply \mathrm{\Delta }{t}_{\mathrm{a}\mathrm{c}\mathrm{t}}= 0.02 s.

Now, one can proceed to calculate the equivalent seismic forces. For the 2D problems, the equivalent seismic forces calculation {\mathit{F}}_{\mathrm{e}\mathrm{q}} requires both the boundary displacements {\mathit{U}}_{\mathrm{b}}^{\mathrm{f}} and stresses {\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{f}} of the free field to be known according to Zhao and Valliappan’s [58], namely,

(17){\mathit{F}}_{\mathrm{e}\mathrm{q}}={\mathit{S}}_{\mathrm{r}\mathrm{b}}{\mathit{U}}_{\mathrm{b}}^{\mathrm{f}}+\mathit{A}{\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{f}},

where subscript ‘r’ and ‘b’ denote the DOFs related to the far field and the interfacing boundary of B in Fig.3, respectively, {\mathit{S}}_{\mathrm{r}\mathrm{b}} is the dynamic impedance matrix of the far field on the boundary B and A the area matrix for converting the surface tractions into nodal forces. To remove the stress-related term \mathit{A}{\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{f}} in the equivalent seismic forces calculation, one needs to start from the original 2D formulation. First, with the use of the equivalent seismic forces {\mathit{F}}_{\mathrm{e}\mathrm{q}} from Eq. (17), the frequency-domain equation of motion of the half-space for the 2D seismic analysis, say, using the FIEM, can be expressed as [39,40]:

(18)\begin{array}{rl}& \left[\begin{array}{cc}{\mathit{S}}_{\mathrm{n}\mathrm{n}}\left(\omega \right)& {\mathit{S}}_{\mathrm{n}\mathrm{b}}\left(\omega \right)\\ {\mathit{S}}_{\mathrm{b}\mathrm{n}}\left(\omega \right)& {\mathit{S}}_{\mathrm{b}\mathrm{b}}\left(\omega \right)+{\mathit{S}}_{\mathrm{r}\mathrm{b}}\left(\omega \right)\end{array}\right]\left[\begin{array}{c}{\mathit{U}}_{\mathrm{n}}\left(\omega \right)\\ {\mathit{U}}_{\mathrm{b}}\left(\omega \right)\end{array}\right]\\ =& \left[\begin{array}{c}0\\ {\mathit{S}}_{\mathrm{r}\mathrm{b}}\left(\omega \right){\mathit{U}}_{\mathrm{b}}^{\mathrm{f}}\left(\omega \right)+\mathit{A}{\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{f}}\left(\omega \right)\end{array}\right],\end{array}

where subscript ‘n’ denotes the DOFs related to the near field, {\mathit{S}}_{\mathrm{n}\mathrm{n}}\left(\omega \right),{\mathit{S}}_{\mathrm{n}\mathrm{b}}\left(\omega \right),{\mathit{S}}_{\mathrm{b}\mathrm{n}}\left(\omega \right),{\mathit{S}}_{\mathrm{b}\mathrm{b}}\left(\omega \right) are the impedance matrices of the near field and interfacing boundary, and {\mathit{U}}_{\mathrm{n}}\left(\omega \right),{\mathit{U}}_{\mathrm{b}}\left(\omega \right) the frequency-domain displacements of the near field and interfacing boundary, respectively.

Then, the stress-related term \mathit{A}{\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{f}} in Eq. (18) will be replaced by the single-layered displacement {\hat{\mathit{U}}}_{\mathrm{s}\mathrm{l}}^{\mathrm{f}}, as an extension of the procedure by Liu et al. [27] to the 2.5D case. Consider the single layer of near-field finite elements that are adjacent to the boundary (denoted by the solid circles in Fig.6, with the DOFs of boundary b circled), as they are related to the equivalent seismic forces {\mathit{F}}_{\mathrm{e}\mathrm{q}}. To distinguish from the DOFs of the near field n, the DOFs of the single layer mentioned here (other than the DOFs of boundary b) will be denoted by \overline{\mathrm{n}}, as marked in Fig.6, with the understanding that the DOFs of \overline{\mathrm{n}} are contained in those of n, i.e., U_{\overline{\mathrm{n}}}\left(\omega \right)\subset U_{\mathrm{n}}\left(\omega \right),. It follows that the equations of motion for the DOFs related to boundary b can be separated from Eq. (18) and written for the free field as follows:

(19)\begin{array}{rl}& \left[\begin{array}{cc}{\mathit{S}}_{\mathrm{b}\overline{\mathrm{n}}}\left(\omega \right)& {\mathit{S}}_{\mathrm{b}\mathrm{b}}\left(\omega \right)+{\mathit{S}}_{\mathrm{r}\mathrm{b}}\left(\omega \right)\end{array}\right]{\left[\begin{array}{cc}{\mathit{U}}_{\overline{\mathrm{n}}}^{\mathrm{f}}\left(\omega \right)& {\mathit{U}}_{\mathrm{b}}^{\mathrm{f}}\left(\omega \right)\end{array}\right]}^{\mathrm{T}}\\ =& {\mathit{S}}_{\mathrm{r}\mathrm{b}}\left(\omega \right){\mathit{U}}_{\mathrm{b}}^{\mathrm{f}}\left(\omega \right)+\mathit{A}{\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{f}}\left(\omega \right),\end{array}

where {\mathit{U}}_{\overline{\mathrm{n}}}^{\mathrm{f}}\left(\omega \right) are the nodal free-field displacements for the single-layer of the finite elements adjacent to the boundary, and {\mathit{U}}_{\mathrm{b}}^{\mathrm{f}}\left(\omega \right) are those of the boundary, both in frequency domain. Here, both {\hat{\mathit{U}}}_{\overline{\mathrm{n}}}^{\mathrm{f}}\left(\omega \right) and {\hat{\mathit{U}}}_{\mathrm{b}}^{\mathrm{f}}\left(\omega \right) contain {\hat{\mathit{U}}}_{\mathrm{s}\mathrm{l}}^{\mathrm{f}}\left(\omega \right). In arriving at Eq. (19), it is recognized that the boundary stiffness {\mathit{S}}_{\mathrm{b}\mathrm{n}}\left(\omega \right) relates only to the DOFs of the single layer \overline{n}, but not to the remaining DOFs of the near field n. Therefore, the product {\mathit{S}}_{\mathrm{b}\mathrm{n}}\left(\omega \right){\mathit{U}}_{\mathrm{n}}\left(\omega \right) in Eq. (18) can be replaced by {\mathit{S}}_{\mathrm{b}\overline{\mathrm{n}}}\left(\omega \right){\mathit{U}}_{\overline{\mathrm{n}}}^{\mathrm{f}}\left(\omega \right). Since the term {\mathit{S}}_{\mathrm{r}\mathrm{b}}\left(\omega \right){\mathit{U}}_{\mathrm{b}}^{\mathrm{f}}\left(\omega \right) appears on both sides of Eq. (19), it can just be eliminated. As a result, the stress-related term \mathit{A}{\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{f}}\left(\omega \right) is:

(20)\mathit{A}{\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{f}}\left(\omega \right)={\mathit{S}}_{\mathrm{b}\stackrel{-}{\mathrm{n}}}\left(\omega \right){\mathit{U}}_{\stackrel{-}{\mathrm{n}}}^{\mathrm{f}}\left(\omega \right)+{\mathit{S}}_{\mathrm{b}\mathrm{b}}\left(\omega \right){\mathit{U}}_{\mathrm{b}}^{\mathrm{f}}\left(\omega \right).

By substituting Eq. (20) into Eq. (17), the equivalent seismic forces {\mathit{F}}_{\mathrm{e}\mathrm{q}} in frequency domain for the 2D problems based only on the displacements input can be given as:

(21){\mathit{F}}_{\mathrm{e}\mathrm{q}}={\mathit{S}}_{\mathrm{b}\stackrel{-}{\mathrm{n}}}{\mathit{U}}_{\stackrel{-}{\mathrm{n}}}^{\mathrm{f}}+\left({\mathit{S}}_{\mathrm{b}\mathrm{b}}+{\mathit{S}}_{\mathrm{r}\mathrm{b}}\right){\mathit{U}}_{\mathrm{b}}^{\mathrm{f}},

For different parts of the near-field boundaries B_{\mathrm{L}},B_{\mathrm{R}},B_{\mathrm{B}} in Fig.3, the equivalent seismic forces {\mathit{F}}_{\mathrm{e}\mathrm{q}} given in Eq. (21) can be easily applied due to absence of the stress-related term. However, it should be noted that the term {\mathit{U}}_{\overline{\mathrm{n}}}^{\mathrm{f}} should be interpreted as the DOFs related to the single-layered finite elements adjacent to the boundary, which differs from the term {\mathit{U}}_{\mathrm{b}}^{\mathrm{f}} for the boundary DOFs.

Finally, the equivalent seismic forces {\mathit{F}}_{\mathrm{e}\mathrm{q}} to be used for 3D oblique incident seismic waves by the 2.5D approach will be computed. As for use by the 2.5D method, the displacements of the free field have been given in frequency-wavenumber (k_{\mathit{\text{z}}}) domain {\hat{\mathit{U}}}_{\mathrm{s}\mathrm{l}}^{\mathrm{f}} (including {\hat{\mathit{U}}}_{\overline{\mathrm{n}}}^{\mathrm{f}} and {\hat{\mathit{U}}}_{\mathrm{b}}^{\mathrm{f}}) in Eqs. (16a) and (16b) for P and SV waves, respectively. Now, one can simply substitute the displacements obtained for the single-layered nodal displacements ({\hat{\mathit{U}}}_{\overline{\mathrm{n}}}^{\mathrm{f}} and {\hat{\mathit{U}}}_{\mathrm{b}}^{\mathrm{f}}) of the near field in Eqs. (16a) and (16b) for P and SV waves, respectively, into the equation in Eq. (21) to obtain:

(22){\mathit{F}}_{\mathrm{e}\mathrm{q}}={\mathit{S}}_{\mathrm{b}\overline{\mathrm{n}}}^{2.5\mathrm{D}}{\hat{\mathit{U}}}_{\overline{\mathrm{n}}}^{\mathrm{f}}+\left({\mathit{S}}_{\mathrm{b}\mathrm{b}}^{2.5\mathrm{D}}+{\mathit{S}}_{\mathrm{r}\mathrm{b}}^{2.5\mathrm{D}}\right){\hat{\mathit{U}}}_{\mathrm{b}}^{\mathrm{f}},

where superscript 2.5D signifies the nautre of the associated dynamic stiffess. The finite and infinite elements to be used to construct the impedance matrices for the 2.5D method will be briefed below.

4 3D seismic analysis by 2.5D finite/infinite element method

In this study, the 2.5D FIEM will be extended to the 3D seismic analysis with oblique incident waves for the advantage of improved computational efficiency. For the 2.5D method, the assumption is that the soil and embedded structure (such as tunnels) are uniform along the z-direction, by which only the 2D profile of the soil and embedded structure on the x–y plane is needed to calculate the 3D wave propagation solutions for the half-space. Nevertheless, the 2.5D formulation differs from the 2D one in that each node of the finite elements contains three displacement DOFs, including the one normal to the profile \left(x,y\right) in Fig.7, while for the 2D formulation, only the two in-plane displacement DOFs are needed.

Since the near-field seismic response to the oblique incident waves is asymmetric with respect to the vertical central line, the origin of the 2.5D element mesh is placed at the left-top corner on the x–y plane in Fig.7. The near field will be modeled by finite elements and far field by infinite elements. The depth and width of the near field on the x–y plane (i.e., finite element mesh range) are set to be 30 m and 60 m along the y- and x-axes, respectively. The longitudinal length is also set to be 60 m along the z-axis (normal to the x–y plane) for the numerical analysis to follow.

4.1 2.5D equation of motion by finite/infinite element approach

A dynamic force acting on the free surface in time domain can be expressed in a general form as [43]

(23)f\left(x,y,\mathit{\text{z}},t\right)=\phi \left(x,y\right)\varphi \left(\mathit{\text{z}}\right)Q\left(t\right),

where \phi \left(x,y\right) is the distribution function on the x–y plane, and \varphi \left(\mathit{\text{z}}\right) the distribution along the z-axis, and Q\left(t\right) the magnitude. By applying the Fourier transform with respect to z and t, the force in Eq. (23) can be converted to the frequency-wavenumber (k_{\mathit{\text{z}}}) domain as:

(24)\hat{f}\left(x,y,k_{\mathit{\text{z}}},\omega \right)=\phi \left(x,y\right)\tilde{\varphi }\left(k_{\mathit{\text{z}}}\right)\tilde{Q}\left(\omega \right),

where \tilde{\varphi }\left(k_{\mathit{\text{z}}}\right) and \tilde{Q}\left(\omega \right) are the Fourier transform of \varphi \left(\mathit{\text{z}}\right) and Q\left(t\right), respectively. Next, the time-space-domain force can be recovered by using the inverse Fourier transform, as:

(25)\begin{array}{rl}f\left(x,y,\mathit{\text{z}},t\right)=& \phi \left(x,y\right){\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\tilde{\varphi }\left(k_{\mathit{\text{z}}}\right)\tilde{Q}\left(\omega \right)\\ & \cdot \mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{i}{k}_{\mathit{\text{z}}}\mathit{\text{z}})\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathrm{i}\omega t\right)\mathrm{d}{k}_{\mathit{\text{z}}}\mathrm{d}\omega .\end{array}

Accordingly, the dynamic force acting on the 3D half-space can be regarded as an integration of the harmonic functions \mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{i}{k}_{\mathit{\text{z}}}\mathit{\text{z}})\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathrm{i}\omega t\right). Let \hat{\mathit{U}}\left(\mathrm{i}\omega \right) denote the frequency response function (FRF) of the displacement to the unit harmonic load \phi \left(x,y\right)\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{i}{k}_{\mathit{\text{z}}}\mathit{\text{z}})\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathrm{i}\omega t\right). The total displacements \mathit{U} to the external force in time domain is

(26)\begin{array}{rl}\mathit{U}\left(x,y,\mathit{\text{z}},t\right)=& {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\tilde{\varphi }\left(k_{\mathit{\text{z}}}\right)\tilde{Q}\left(\omega \right)\hat{\mathbf{U}}\left(\mathrm{i}\omega \right)\\ & \cdot \mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{i}{k}_{\mathit{\text{z}}}\mathit{\text{z}})\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathrm{i}\omega t\right)\mathrm{d}{k}_{\mathit{\text{z}}}\mathrm{d}\omega .\end{array}

Using the 2.5D method, the displacement response \hat{\mathit{U}} of the near field and boundary in frequency-wavenumber (k_{\mathit{\text{z}}}) domain can be computed using the plane elements consisting of three DOFs at each node on the x–y plane, which makes it possible to separate variable z from the 2D profile \left(x,y\right).

The elements chosen herein are the quadratic 8-node (Q8) finite element and degenerated 5-node (Q5) infinite element in frequency domain. For each element, the global coordinates \left(x,y\right) and displacements \hat{\mathit{U}} of each point on the x–y plane can be interpolated by the nodal coordinates \left(x_i,y_i\right) and displacements {\hat{\mathit{U}}}_i, respectively, via the respective shape functions \left(M_i,N_i\right). As the Q8 finite element is available elsewhere, only the shape functions for the coordinates and displacements \left(M_i,N_i\right) of the Q5 infinite element are given below:

\begin{array}{}\text{(27a)}& M_1=-\frac{1}{2}\left(\xi -1\right)\left(\eta -1\right)\eta ,\end{array}

\begin{array}{}\text{(27b)}& M_2=\left(\xi -1\right)\left(\eta -1\right)\left(\eta +1\right),\end{array}

\begin{array}{}\text{(27c)}& M_3=-\frac{1}{2}\left(\xi -1\right)\left(\eta +1\right)\eta ,\end{array}

\begin{array}{}\text{(27d)}& M_4=\frac{1}{2}\xi \left(\eta +1\right),\end{array}

\begin{array}{}\text{(27e)}& M_5=-\frac{1}{2}\xi \left(\eta -1\right).\end{array}

\begin{array}{}\text{(28a)}& N_1=-\frac{\left(\eta -1\right)\eta }{2}P\left(\xi \right),\end{array}

\begin{array}{}\text{(28b)}& N_2=-\left(\eta -1\right)\left(\eta +1\right)P\left(\xi \right),\end{array}

\begin{array}{}\text{(28c)}& N_3=\frac{\left(\eta +1\right)\eta }{2}P\left(\xi \right),\end{array}

(29)P\left(\xi \right)={\mathrm{e}}^{-{\alpha }_{\mathrm{l}}\xi }{\mathrm{e}}^{-{\mathrm{i}k}_{\mathrm{l}}\xi },

where \left(\xi ,\eta \right) denotes the local coordinates, P\left(\xi \right) is the wave propagation function, the term with {\alpha }_{\mathrm{l}} denotes the decay factor for geometric attenuation to infinity and the term with k_{\mathrm{l}} the propagation of waves.

For the 2.5D approach, due to addition of the out-of-plane DOF of each node, the stress-strain matrices {\mathit{B}}^{2.5\mathrm{D}} in the local coordinates of both the finite and infinite elements no longer possess the property of symmetry, i.e.,

\begin{array}{}\text{(30a)}& {\mathit{B}}^{2.5\mathrm{D}}=\mathit{L}\mathit{N},\end{array}

\begin{array}{}\text{(30b)}& \mathit{L}={\left[\begin{array}{cccccc}{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}& 0& 0& 0& -\mathrm{i}{k}_{\mathit{\text{z}}}& {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }y}}\\ 0& {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }y}}& 0& -\mathrm{i}{k}_{\mathit{\text{z}}}& 0& {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}\\ 0& 0& -\mathrm{i}{k}_{\mathit{\text{z}}}& {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }y}}& {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}& 0\end{array}\right]}^{\mathrm{T}},\end{array}

\begin{array}{}\text{(30c)}& \mathit{N}=\left[\begin{array}{ccccccc}{N}_1& 0& 0& & N_j& 0& 0\\ 0& N_1& 0& \dots & 0& N_j& 0\\ 0& 0& N_1& & 0& 0& N_j\end{array}\right],\end{array}

where L and N are the matrixes for the deferential operator and the shape function for the displacements, respectively, and subscript j=8\mathrm{o}\mathrm{r}3 for the finite or infinite elements, respectively.

In assembly, the impedance matrices {\mathit{S}}^{2.5\mathrm{D}} (= {\mathit{K}}^{2.5\mathrm{D}}-{{\omega }^2\mathit{M}}^{2.5\mathrm{D}}) should be given in the global coordinates, where {\mathit{K}}^{2.5\mathrm{D}} and {\mathit{M}}^{2.5\mathrm{D}} denote the global stiffness and the mass matrices, respectively. By substituting the equivalent seismic forces {\mathit{F}}_{\mathrm{e}\mathrm{q}}\left({\omega ,k}_{\mathit{\text{z}}}\right) derived in Eq. (22) for a unit ground acceleration (for given \omega), the 2.5D equation of motion for the 3D soil seismic analysis can be presented in frequency-wavenumber (k_{\mathit{\text{z}}}) domain as:

(31){\mathit{S}}^{2.5\mathrm{D}}\left({\omega ,k}_{\mathit{\text{z}}}\right){\hat{\mathit{U}}}_j\left({\omega ,k}_{\mathit{\text{z}}}\right)={\mathit{F}}_{\mathrm{e}\mathrm{q}}\left({\omega ,k}_{\mathit{\text{z}}}\right),

where subscript j denotes P or SV waves. Accordingly, the nodal displacements {\hat{\mathit{U}}}_j\left({\omega ,k}_{\mathit{\text{z}}}\right) of the near field and boundary for each finite element in frequency-wavenumber (k_{\mathit{\text{z}}}) domain for a unit ground acceleration (for given \omega) can be obtained, which can be interpolated to yield the displacements \hat{\mathit{U}}\left({\omega ,k}_{\mathit{\text{z}}}\right) of any point within the near field of concern.

4.2 Relevant parameters for the finite/infinite elements

For 3D seismic analysis, some crucial parameters of the numerical elements need to be modified and selected when using the 2.5D method. First, the wave number k_{\mathit{\text{z}}} used in Eq. (30b) should be set equal to the wave numbers k_{\mathit{\text{z}}\mathrm{P}}\mathrm{a}\mathrm{n}\mathrm{d}{k}_{\mathit{\text{z}}\mathrm{S}} along the z-axis as given by the exact solutions in Eqs. (1b) and (7b) for the P and SV waves, respectively, to be consistent with the equivalent seismic forces calculation {\mathit{F}}_{\mathrm{e}\mathrm{q}}. Then, the decay factor {\alpha }_{\mathrm{l}} in Eq. (29) can be modified using the analogy to the case for the moving loads in Ref. [43], that is:

\begin{array}{}\text{(32a)}& {\alpha }_{\mathrm{P}}=\frac{1}{2R}\left(1+\frac{k_{\mathit{\text{z}}}^2}{k_{\mathit{\text{z}}}^2+k_{\mathrm{P}}^2}\right),\end{array}

\begin{array}{}\text{(32b)}& {\alpha }_{\mathrm{S}}=\frac{1}{2R}\left(1+\frac{k_{\mathit{\text{z}}}^2}{k_{\mathit{\text{z}}}^2+k_{\mathrm{S}}^2}\right),\end{array}

\begin{array}{}\text{(32c)}& {\alpha }_{\mathrm{R}}=\frac{1}{2R}\frac{k_{\mathit{\text{z}}}^2}{k_{\mathit{\text{z}}}^2+k_{\mathrm{R}}^2},\end{array}

where R is the distance between the center point of the free surface and any point on the near-field boundary (Fig.7), k_{\mathit{\text{z}}} is set to k_{\mathit{\text{z}}\mathrm{P}} for P waves, and to k_{\mathit{\text{z}}\mathrm{S}} for SV waves. The modified decay factors \left({\alpha }_{\mathrm{P}},{\alpha }_{\mathrm{S}},{\alpha }_{\mathrm{R}}\right) correspond to the region where P, SV, and Rayleigh (R) waves are dominant, respectively. Besides, the wave number k_j used in Eq. (29) should be modified as:

(33)k_j^{\prime }=\sqrt{k_j^2-k_{\mathit{\text{z}}}^2},

where the subscript j represents the R, P, and SV waves, respectively.

Similar to analysis of the 2D problems in Ref. [41], the R-, SV-, and P-wave numbers, i.e., k_{\mathrm{R}}^{\prime }, k_{\mathrm{S}}^{\prime }, and k_{\mathrm{P}} are chosen for the infinite elements set on the near- and below-surface boundaries on the two sides and deep bottom boundary, for incident P waves and for incident SV waves with sub-critical angles on the wave propagation plane. For incident SV waves with super-critical angles on the wave propagation plane, the SV-wave number k_{\mathrm{S}}^{\prime } is chosen for both the below-surface boundaries on two sides and the deep bottom boundary.

4.3 Seismic analysis procedure for 3D arbitrary incident waves

First, for the range of frequencies considered (Fig.5), i.e., {\omega }_{\mathrm{a}\mathrm{c}\mathrm{t}}= 25 Hz, a frequency increment, i.e., \mathrm{\Delta }\omega = 0.1 Hz ( = 0.63 rad/s) [46] is used to divide the earthquake spectrum into a number of intervals. Then, for each frequency \omega within the range, one can first solve the near-field and boundary displacements \hat{\mathit{U}} in frequency-wavenumber (k_{\mathit{\text{z}}}) domain by Eq. (31) for the P or SV waves, and obtain the displacements in frequency-space \left(x,y,\mathit{\text{z}}\right) domain by the inverse Fourier transform to Eq. (A7b) in Electronic Supplementary Materials with respect to k_{\mathit{\text{z}}\mathrm{P}} or k_{\mathit{\text{z}}\mathrm{S}} for P or SV waves, respectively, as:

(34)\tilde{\mathit{U}}\left(x,y,\mathit{\text{z}},\omega \right)=\hat{\mathit{U}}\left(x,y,k_{\mathit{\text{z}}},\omega \right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{i}{k}_{\mathit{\text{z}}}\mathit{\text{z}}\right),

in a way similar to derivation of Eq. (16). Then, one can calculate the acceleration responses \ddot{\tilde{\mathit{U}}}\left(\omega \right) in frequency-space \left(x,y,\mathit{\text{z}}\right) domain by multiplying \tilde{\mathit{U}}\left(\omega \right) by {-\omega }^2. Repeating the above procedure for the next frequency \omega +\mathrm{\Delta }\omega, one can obtain the component responses for each frequency \omega.

The above solutions have been given in frequency domain. The inverse transformation from frequency to time domain based on the formula in Eq. (A7b) in Electronic Supplementary Materials can be briefly summarized as follows [41]: 1) Solving the spectral responses by multiplying the component responses of the near field and boundary by the corresponding spectral amplitude of the earthquake spectrum (see Fig.5). 2) Employing the IFFT with respect to \omega by adding up all the spectral responses for the frequency \omega within the actual frequency range, i.e., 0–25 Hz, with increment \mathrm{\Delta }\omega (rad/s). By so doing, the integration with respect to \omega has been replaced by the summation of the spectral responses for each frequency \omega.

As a side note, the time-history responses for the region of concern were calculated with the actual time increment \mathrm{\Delta }{t}_{\mathrm{a}\mathrm{c}\mathrm{t}}=\gamma \mathrm{\Delta }t, where \mathrm{\Delta }t is the time increment of the Chi-Chi Earthquake data and \gamma = 4 is the ratio of the theoretical and actual frequency ranges used (see Subsection 3.2 for the explanation).

One benefit with the 2.5D FIEM for seismic analysis is that not only the reflective and refractive waves on the x–y plane, but the waves transmitting along the tunnel axis can be duly considered, similar to the moving train problems [43–45,47–49]. Since the 2.5D model consists only of the 2D plane strain elements, the overall computational cost is very small compared with the 3D modeling, due to the drastic reduction in elements, nodes and number of DOFs. However, the accuracy of the 2.5D FIEM in simulating 3D wave propagation can still be guaranteed. The other advantage is that infinite elements are topographically compatible with the main-stream finite elements. Therefore, the 2.5D FIEM is likely to be easily adopted by practicing engineers. It should be added that the present 2.5D FIEM for the seismic analysis is also applicable to the layered half-space, once the exact responses of the layered free-field are solved, which is part of the future work.

As a side remark, the total procedure for the 2.5D soil seismic analysis has been implemented in a self-designed program. The infinite elements used herein can be easily included in the existing finite element software, since they are a degenerated of the finite elements. The only new feature for the 2.5D formulation is that each node of the finite/infinite elements contains three displacement DOFs, including the one normal to the profile \left(x,y\right) in Fig.7, unlike the 2D formulation that has only two in-plane displacement DOFs.

5 Verification of the 2.5D method for 3D seismic analysis

5.1 Verification of free-field solution in frequency domain with Wolf’s

First, the reliability of the procedure for ground motion input by the 2.5D FIEM will be verified against Wolf’s solutions [4]. The range of frequency included in calculation is up to 25 Hz, identical to the actual frequency range of the earthquake spectra in Fig.5. The soil properties adopted are: density \rho = 2000 kg/m³, Poisson’s ratio \nu = 0.33, elastic modulus E_{\mathrm{s}}= 100 MPa, and corresponding shear waves velocity c_{\mathrm{S}}= 137.1 m/s, with the hysteresis damping effect ignored.

For selection of the finite element size L, concerning the criterion of element size L\le {\lambda }_{\mathrm{S}}/5 [40], where {\lambda }_{\mathrm{S}} denotes the shear wavelength, the required element size should be determined by the minimum shear wavelength {\left({\lambda }_{\mathrm{S}}\right)}_{\mathrm{m}\mathrm{i}\mathrm{n}}, corresponding to the upper limit of the calculated frequency range (i.e., {\omega }_{\mathrm{a}\mathrm{c}\mathrm{t}}=25\mathrm{H}\mathrm{z} and {\left({\lambda }_{\mathrm{S}}\right)}_{\mathrm{m}\mathrm{i}\mathrm{n}}= 5.5 m). Based on this, the finite element size L can be taken as 1 m. As for selection of the mesh range R, no requirement is needed for the narrow frequency range of the spectra in seismic analysis [40]. Thus, the mesh range of R = 30 m in Fig.7 is adopted merely for better observation of the near-field behavior of the soil. See the Ref. [40] for explanation of the criteria for selecting the mesh range R and finite element size L.

Two horizontal incident angles, i.e., {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}}\left(\mathrm{o}\mathrm{r}{\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}}\right)= 0° or {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}}\left(\mathrm{o}\mathrm{r}{\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}}\right)= 90°, are considered for P or SV waves. For a unit vertical acceleration of P waves (or horizontal of SV waves) applied at the origin of the x–y plane (Fig.7) for each frequency \omega (Section 3.1), the ground displacements of the points OP-1, OP-2, OP-3, can be computed using the present 2.5D procedure. To compare with Wolf’s solutions, the ground displacement amplitudes computed have been normalized with respect to the displacement magnitude of the input excitation at the origin, i.e., U_{\mathrm{P}} for the incident P waves, and U_{\mathrm{S}} for the SV waves. The normalized displacements for points OP-1, OP-2, OP-3, i.e., \left|\tilde{u}/U_{\mathrm{P}}\right|, \left|\tilde{v}/U_{\mathrm{P}}\right|, \left|\tilde{w}/U_{\mathrm{P}}\right|, \left|\tilde{u}/U_{\mathrm{S}}\right|, \left|\tilde{v}/U_{\mathrm{S}}\right|, \left|\tilde{w}/U_{\mathrm{S}}\right|, have been compared with Wolf’s against the vertical incident angles (0°−90°) in Fig.8 and Fig.9 for P waves for the horizontal incident angles of {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}} = 0° and 90°, respectively, and similar for the SV waves in Fig.10 and Fig.11.

As revealed by these figures, all the solutions calculated by the proposed 2.5D method for the horizontal incident angles of 0° and 90° match excellently with the Wolf’s solutions for the vertical incident angles (0°−90°) on the wave propagation plane for both P and SV waves. This has verified the reliability of the present 2.5D method for the 3D seismic response analysis. Of interest is that abrupt variation occurs on the soil response around the critical angle on the wave propagation plane of 30° (calculated by {\theta }_{\mathrm{S}\mathrm{i}}^{\mathrm{c}\mathrm{r}}={\mathrm{s}\mathrm{i}\mathrm{n}}^{-1}\left(1/V\right) [41]) in Fig.10 and Fig.11, for incident SV waves, similar to the phenomenon observed for 2D problems [41].

It should be noted that the longitudinal displacements \left(\left|\tilde{w}/U_{\mathrm{P}}\right|\mathrm{o}\mathrm{r}\left|\tilde{w}/U_{\mathrm{S}}\right|\right) vanish for {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}} \left(\mathrm{o}\mathrm{r}{\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}}\right)= 0°, and the horizontal displacements \left(\left|\tilde{u}/U_{\mathrm{P}}\right|\mathrm{o}\mathrm{r}\left|\tilde{u}/U_{\mathrm{S}}\right|\right) vanish for {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}} \left(\mathrm{o}\mathrm{r}{\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}}\right)= 90°, for the P or SV waves, which thus were not shown herein. From this analysis, it is confirmed that the mesh range and element size used (R = 30 m, L = 1 m) are adequate for the soil properties adopted. Meanwhile, the criterion of L\le {\lambda }_{\mathrm{S}}/5 for element size and the mesh range adopted are deemed adequate for seismic analysis, as pointed out in Ref. [41]. Unless noted otherwise, both L = 1 m and R = 30 m will be adopted in the analyses to follow.

In the above comparison with existing solutions, the vertical incident angle {\theta }_{\mathrm{P}\mathrm{i}\mathrm{v}} for P waves (or {\theta }_{\mathrm{S}\mathrm{i}\mathrm{v}} for SV waves) is allowed to vary from 0° to 90°. In fact, for {\theta }_{\mathrm{P}\mathrm{i}\mathrm{v}} (or {\theta }_{\mathrm{S}\mathrm{i}\mathrm{v}}) approaching 90° or for the body waves close to the x–z plane (Fig.2), the surface waves or Rayleigh waves are going to play a role and may change the trend of response on the soil surface. Such an effect has not been addressed in this paper.

5.2 Verification of scattered-field solution in frequency domain with de Barros and Luco’s

In this section, the present method will be used to study the 3D response of the underground space enclosing a cylindrical cavity (Fig.12), and compared with the solution by de Barros and Luco [15]. To this end, the buried depth ratio chosen is H/r=2, the buried depth is H=10 m and the radius r=5 m. The material properties of the soil selected are: density \rho = 2000 kg/m³, Poisson’s ratio \nu = 0.33, elastic modulus E_{\mathrm{s}}= 100 MPa, hysteresis damping \beta = 0.01 and corresponding shear waves velocity c_{\mathrm{S}}= 137.1 m/s. Also, the dimensionless frequency \eta =2rf/c_{\mathrm{S}} is 0.5, corresponding to the frequency of f= 7.5 Hz.

Based on Ref. [15], the unit vertical or horizontal acceleration is applied at the point (30,0,0) in Fig.12 for P or SV waves, and the horizontal and vertical incident angles, i.e., {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}} and {\theta }_{\mathrm{P}\mathrm{i}\mathrm{v}} for P waves and {\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}} and {\theta }_{\mathrm{S}\mathrm{i}\mathrm{v}} for SV waves in Fig.2, are taken as 30° and 45°, respectively. Then, the ground displacements on the x–y plane, i.e., \left|\tilde{u}/U_{\mathrm{P}}\right|,\left|\tilde{v}/U_{\mathrm{P}}\right|,\left|\tilde{w}/U_{\mathrm{P}}\right| for P waves and \left|\tilde{u}/U_{\mathrm{S}}\right|,\left|\tilde{v}/U_{\mathrm{S}}\right|,\left|\tilde{w}/U_{\mathrm{S}}\right| for SV waves, calculated by present 2.5D approach, along with de Barros and Luco’s [15], were plotted in Fig.13 and Fig.14 for P and SV waves, respectively. Evidently, for both P and SV waves, the present 2.5D solutions are in excellent agreement with de Barros and Luco’s, indicating that the scattered waves generated by the embedded cavity can be absorbed well by the infinite elements using the present 2.5D method. It should be added that for high-frequency scattered waves, which may be generated by sources like the moving train loads (close to 100 Hz), they can be also absorbed well by the infinite elements, as long as the element size criteria for finite elements are satisfied [47–49], due to the fact that high-frequency waves decay faster than the low ones.

5.3 Verification of free-field solution in time domain

The properties adopted for the soil are: density \rho = 2000 kg/m³, elastic Poisson’s ratio \nu = 0.25, modulus E_{\mathrm{s}}= 100 MPa, and hysteretic damping ratio \beta = 0.02. For the soil properties adopted herein, the mesh range and element size selected in Subsection 5.1 (R= 30 m, L= 1 m) remain adequate, based on the criteria for the element size and the mesh range proposed in Ref. [40]. and elaborated in Subsection 5.1 in this paper. The accuracy of the entire 2.5D seismic analysis procedure, including the transformation of response from frequency to time domain for 3D arbitrary incident seismic waves, as outlined in Subsection 4.3, will be verified herein. The horizontal and vertical incident angles {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}} and {\theta }_{\mathrm{P}\mathrm{i}\mathrm{v}}, respectively, for P waves, and {\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}} and {\theta }_{\mathrm{S}\mathrm{i}\mathrm{v}} for SV waves (see Fig.2(a) and Fig.2(b)) are all chosen to be 50°. By specifying the Chi-Chi Earthquake motions at the ground origin and using the theoretical framework given in Section 3 for incident waves, the equilibrium earthquake forces on the near-field boundary can be computed. With this, the ground motions can be inversely computed for point OP-1 on the ground of the x–y plane in Fig.7, where the original Chi-Chi Earthquake motions were imposed.

Both the original and inversely computed time-histories of the Chi-Chi Earthquake (see Fig.4 for original) at point OP-1 on the x–y plane were plotted in Fig.15(a) and Fig.15(b), for SV and P waves, respectively. Evidently, the inversely computed solutions and the original excitation at point OP-1 on the x–y plane (ground motion input point) for the Chi-Chi Earthquake match excellently. It is indicated that the ground motion input by the equivalent seismic forces and the IFFT for calculating the time-history responses are correct and reliable, which will be adopted in the time-domain analysis in the parametric study to follow.

6 Parametric study of 3D soil seismic responses

6.1 3D soil seismic responses in frequency domain for P waves

Using the same soil model in Fig.7, the ground input in Subsection 3.1, and the soil properties in Subsection 5.3, the effect of the horizontal incident angle {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}} (Fig.2(a)) on the soil response will be first assessed. Specifically, the element size and the mesh range selected, i.e., L = 1 m and R = 30 m, will be adopted herein for the same soil properties as in Subsection 5.3 used since their reliability has been ensured in Section 5 for both the frequency- and time-domain analyses. For fixed vertical incident angle, i.e., {\theta }_{\mathrm{P}\mathrm{i}\mathrm{v}}= 50°, the normalized vertical displacement amplitudes \left|\tilde{v}/U_{\mathrm{P}}\right| is also fixed. The horizontal and longitudinal displacement amplitudes in normalized form (\left|\tilde{u}/U_{\mathrm{P}}\right|, \left|\tilde{w}/U_{\mathrm{P}}\right|) vs. the horizontal incident angle {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}} (0° to 90°) of the observation points along the soil depth (Fig.7) on the x–y plane under the frequency f of 25 Hz were plotted in Fig.16(a) and Fig.16(b), respectively.

From Fig.16, one observes that the horizontal displacements \left|\tilde{u}/U_{\mathrm{P}}\right| on both the ground and underground attenuate with increasing horizontal incident angle {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}}, while the longitudinal ones \left|\tilde{w}/U_{\mathrm{P}}\right| increase inversely. Next, since the wave propagation plane (see Fig.2(a)) moves away from the x–y plane with increasing incident angle {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}}, the soil response along the x-axis decreases, while that along the z-axis increases gradually.

To further study the distribution of the 3D soil response, the frequency f= 25 Hz and incident angles of {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}}, {\theta }_{\mathrm{P}\mathrm{i}\mathrm{v}}= 50° were considered. The three displacements (\left|\tilde{u}/U_{\mathrm{P}}\right|,\left|\tilde{v}/U_{\mathrm{P}}\right|,\left|\tilde{w}/U_{\mathrm{P}}\right|) computed for the free surface along the x- and y-axes and normal to the x–y plane were plotted in Fig.17(a). Besides, to enhance the visual impression on the responses along the depth, Fig.17(b) has been devised such that the horizontal axis represents the displacements, and the vertical axis the depth. From these figures, all the displacements (\left|\tilde{u}/U_{\mathrm{P}}\right|,\left|\tilde{v}/U_{\mathrm{P}}\right|,\left|\tilde{w}/U_{\mathrm{P}}\right|) show a trend of monotonic decay along the horizontal x-axis, while they all fluctuate synchronously along the soil depth (y-axis) with relatively stable amplitudes. Besides, the trend of displacements along the longitudinal z-axis is similar to that along the horizontal x-axis, due to the fact that the attenuation term \mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{i}{k}_{\mathit{\text{z}}}\mathit{\text{z}}\right) about variable z is accompanied by the term \mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{i}{k}_{x}x\right) about variable x in Eq. (6).

6.2 3D soil seismic responses in frequency domain for SV waves

The properties adopted of the soil are identical those for incident P waves, together with the same ground input in Subsection 3.1. For the SV waves, the frequency considered is f= 25 Hz and the vertical incident angle is {\theta }_{\mathrm{S}\mathrm{i}\mathrm{v}}= 50°. The horizontal and longitudinal displacements (\left|\tilde{u}/U_{\mathrm{S}}\right|, \left|\tilde{w}/U_{\mathrm{S}}\right|) vs. the horizontal incident angle {\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}} (0° to 90°) of the observation points along the soil depth (Fig.7) on the x–y plane were plotted in Fig.18(a) and Fig.18(b), respectively. For the SV waves, the critical angle on the wave propagation plane calculated by {\theta }_{\mathrm{S}\mathrm{i}}^{\mathrm{c}\mathrm{r}}={\mathrm{s}\mathrm{i}\mathrm{n}}^{-1}\left(1/V\right) is 35.26° [41]. As can be seen from Fig.18, the trends for both the horizontal and longitudinal displacements (\left|\tilde{u}/U_{\mathrm{S}}\right|, \left|\tilde{w}/U_{\mathrm{S}}\right|) with increasing horizontal incident angle {\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}} are similar to the case for P waves presented above. Namely, the horizontal response attenuates with increasing horizontal incident angle {\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}}, while the longitudinal one \left|\tilde{w}/U_{\mathrm{S}}\right| increases inversely.

Similarly, to investigate the spatial distribution of displacements under incident SV waves, the frequency f= 25 Hz and incident angles of {\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}},{\theta }_{\mathrm{S}\mathrm{i}\mathrm{v}}= 50° were adopted, which means a super-critical condition. The displacements (\left|\tilde{u}/U_{\mathrm{S}}\right|,\left|\tilde{v}/U_{\mathrm{S}}\right|,\left|\tilde{w}/U_{\mathrm{S}}\right|) computed for the free surface and along the soil depth (vertical central line) on the x–y plane were plotted in Fig.19(a) and Fig.19(b), respectively. Obviously, except for the trends similar to the ones for P waves, the vertical displacement \left|\tilde{v}/U_{\mathrm{S}}\right| appears to be significantly larger than both the horizontal and longitudinal ones (\left|\tilde{u}/U_{\mathrm{S}}\right|, \left|\tilde{w}/U_{\mathrm{S}}\right|) on the ground. Moreover, for increasing depth, both the horizontal and longitudinal displacements (\left|\tilde{u}/U_{\mathrm{S}}\right|,\left|\tilde{w}/U_{\mathrm{S}}\right|) fluctuate synchronously, but they differ by half a period of fluctuation with the vertical ones \left|\tilde{v}/U_{\mathrm{S}}\right|, as revealed by Fig.19(b). Both the phenomena observed above are consistent with those of the 2D problems [41] for SV waves with vertical incident angles {\theta }_{\mathrm{S}\mathrm{i}\mathrm{v}} larger than the critical value, due to the fact that the reflected P waves (Fig.2(b)) have been converted into the surface waves along the ground [41]. This means that the original restraining effect of the reflected P waves on the vertical ground responses has been changed into that of transformed surface waves on the horizontal ones, due to the conversion of waveforms for the reflected P and SV waves. Consequently, nonsynchronous fluctuation was observed along the depth (Fig.19(b)) in the vertical direction (y-axis) and other two directions (x- and z-axes) due to the reversal of the restraining effect for occurrence of the surface waves.

6.3 3D soil seismic responses in time domain for P waves

The vertical motion of the Chi-Chi Earthquake in Fig.4(b) is chosen as the ground input for P waves at point (0,0,0) in Fig.3, with other soil parameters kept identical to those in Subsection 5.3. To assess the trend of wave propagation on the free surface in time domain, the analysis will be conducted for some specific instants prior to the full development of the earthquake to avoid the mingling of complicated reflected and refracted waves within the soil body. For instance, all the instants selected below, 0.45, 0.47, 0.49 s, belong to the initial stage of the seismic excitation span of 0–90 s, when the PGA is not developed yet (Fig.4). This is the reason why the acceleration amplitudes shown in the figures below are all much smaller than the PGA of the earthquake.

First, for vertical incident angle {\theta }_{\mathrm{P}\mathrm{i}\mathrm{v}} = 50°, the transient wave forms for vertical ground acceleration a_y at 0.45 s were plotted for three horizontal incident angles {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}} = 15°, 50°, 85° in Fig.20. Clearly, with increasing horizontal incident angle {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}}, the wave propagation plane (Fig.2(a)) gradually deviates from the x–y plane and simultaneously becomes close to the y–z plane. In other words, the direction of wave propagation on the x–z plane changes from the one parallel to the x-axis to the z-axis, which conforms to the normal direction of the transient wave forms on the x–z plane as indicated by red arrows in Fig.20. Since the incident waves imposed on the wave propagation plane “rotate” with the plane, the direction of incident waves on the x–z plane (free surface of soil) conforms to the normal direction of the transient wave forms.

To obtain the trend of ground acceleration over time for horizontal and vertical incident angles {\theta }_{\mathrm{P}\mathrm{i}\mathrm{h}},{\theta }_{\mathrm{P}\mathrm{i}\mathrm{v}} of 50°, the transient vertical ground acceleration a_y on the x–z plane at three moments t = 0.45, 0.47, 0.49 s were plotted in Fig.21. One observes intuitively that the transient wave form advances along its normal direction (red arrows) or the direction of incident waves on the x–z plane as time goes by. In addition, from the coordinate data of the marked points for the same wave crest in Fig.21, one also observes that the acceleration amplitude attenuates as well along the wave propagation direction. Consequently, for 3D incident P waves, it is inferred that the direction of wave propagation, direction of response attenuation, and direction of incident waves on the x–z plane coincide simultaneously.

6.4 3D soil seismic responses in time domain for SV waves

In this section, the horizontal motion of the Chi-Chi Earthquake in Fig.4(a) is chosen for SV waves as the horizontal ground input at point (0,0,0) in Fig.3. Similarly, three horizontal incident angles are considered, {\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}} = 15°, 50°, 85°. For vertical incident angle {\theta }_{\mathrm{S}\mathrm{i}\mathrm{v}}= 50°, the transient horizontal ground acceleration a_x on the x–z plane at 0.45 s for the three cases were plotted in Fig.22. Similar observations exist as those for P waves. Namely, the direction of incident waves on the x–z plane conforms to the normal direction of transient wave forms or wave propagation on the x–z plane. Moreover, the shape of transient wave forms on the x–z plane is reflective of the specific earthquake used as the ground input via the origin.

Finally, for horizontal and vertical incident angles {\theta }_{\mathrm{S}\mathrm{i}\mathrm{h}},{\theta }_{\mathrm{S}\mathrm{i}\mathrm{v}}= 50°, the transient horizontal accelerations on the x–z plane at t = 0.45, 0.47, 0.49 s were plotted in Fig.23. Except for the same characteristics of acceleration on the ground as those for P waves, the ground waves propagate more slowly than that for P waves, as can be intuitively appreciated by comparing the marked positions for t = 0.45, 0.47, 0.49 s on the z-axis with those in Fig.23. It is indicated that stronger ground responses will persist under the incident SV waves within the same propagation range on the ground.

7 Concluding remarks

For 3D arbitrary incident P or SV waves, the theory and procedure are presented for the 3D seismic analysis of a viscoelastic half-space. Using the exact solutions derived of the free-field displacements, equivalent seismic forces acting on the near-field boundary are calculated for earthquake spectrum for the 2.5D FIEM. Subsequently, with the asymmetric model used to account for oblique incident seismic waves, the 3D soil seismic analysis procedure is established. The reliability of the proposed method has been verified against the solutions of Wolf’s and de Barros and Luco’s, as well as those for 3D inversely calculated ground motions. Based on the theory, earthquake, and soil properties adopted in analysis, the following conclusions are drawn for the soil responses to 3D arbitrary incident seismic waves.

1) Effect of incident angles: a) For fixed vertical incident angle, the ground or underground horizontal displacements decrease, while longitudinal displacements increase monotonously with increasing horizontal incident angle, for both P and SV waves; b) Abrupt variation in soil response occurs for SV waves around the critical incident angle on the wave propagation plane.

2) Spatial irregularity: a) All the displacement components decay monotonously along the horizontal and longitudinal directions, for both P and SV waves; b) Along the soil depth, all the displacement components fluctuate synchronously for P waves, while the fluctuation of half a period exists between the vertical displacement amplitude and the other two displacement components for SV waves with super-critical incident angles, both with little attenuation.

3) Time domain: a) The propagating direction of surface waves, the attenuation direction of response amplitude, and the direction of incident waves on the x–z plane are all coincident for both P and SV waves; b) The ground waves generated for SV waves propagate more slowly than for P waves; c) The persistence of ground responses is stronger for SV waves than for P waves along the propagation direction of surface waves.

Research is undergoing on the seismic analysis of problems closer to engineering practices, such as a half-space enclosing underground structures (i.e., tunnels) or overlain with sea water under arbitrary incident waves, aimed at evaluating the catastrophic impact of earthquakes.

The contribution of the paper is to establish the theoretical framework for the seismic analysis involving 3D arbitrary incident waves, which is not available in the literature. The half-space is assumed to be homogeneous viscoelastic at this moment, but it can be easily extended for nonelastic materials, e.g., anisotropic, poroelastic, or inelastic soils, to meet the practical needs, if proper constitutive laws and/or yield criteria are employed, subject to the need of iteration to obtain the final solution.

References

[1] Sanchez-Sesma F J. Diffraction of elastic waves by three-dimensional surface irregularities. Bulletin of the Seismological Society of America, 1983, 73(6A): 1621–1636

[2] Lee J J, Langston C A. Wave propagation in a three-dimensional circular basin. Bulletin of the Seismological Society of America, 1983, 73(6A): 1637–1653

[3] Lee V W. Three-dimensional diffraction of plane P, SV & SH waves by a hemispherical alluvial valley. International Journal of Soil Dynamics and Earthquake Engineering, 1984, 3(3): 133–144

[4] WolfJ P. Dynamic Soil−structure Interaction. New Jersey: Prentice-Hall, 1985

[5] Wong K C, Datta S K, Shah A H. Three-dimensional motion of buried pipeline. I: Analysis. Journal of Engineering Mechanics, 1986, 112(12): 1319–1337

[6] Luco J E, Wong H L. Seismic response of foundations embedded in a layered half-space. Earthquake Engineering & Structural Dynamics, 1987, 15(2): 233–247

[7] Moore I D, Guan F. Three-dimensional dynamic response of lined tunnels due to incident seismic waves. Earthquake Engineering & Structural Dynamics, 1996, 25(4): 357–369

[8] Liang J W, Ba Z. Exact dynamic stiffness matrices of 3-D layered site and its Green’s functions. Earthquake Engineering and Engineering Vibration, 2007, 5: 7–17

[9] Zhang Z Q, Pan E. Time-harmonic response of transversely isotropic and layered poroelastic half-spaces under general buried loads. Applied Mathematical Modelling, 2020, 80: 426–453

[10] Zhou J C, Pan E, Lin C P. A novel method for calculating dislocation Green’s functions and deformation in a transversely isotropic and layered elastic half-space. Engineering Analysis with Boundary Elements, 2023, 152: 22–44

[11] Bobet A, Yu H, Tiwari N. Seismic response of shallow circular openings to Rayleigh waves. Tunnelling and Underground Space Technology, 2023, 135: 105036

[12] Mossessian T K, Dravinski M. Amplification of elastic waves by a three dimensional valley. Part 1: Steady state response. Earthquake Engineering & Structural Dynamics, 1990, 19(5): 667–680

[13] Mossessian T K, Dravinski M. Amplification of elastic waves by a three dimensional valley. Part 2: Transient response. Earthquake Engineering & Structural Dynamics, 1990, 19(5): 681–691

[14] Luco J E, Wong H L, de Barros F C P. Three-dimensional response of a cylindrical canyon in a layered half-space. Earthquake Engineering & Structural Dynamics, 1990, 19(6): 799–817

[15] de Barros F C P, Luco J E. Diffraction of obliquely incident waves by a cylindrical cavity embedded in a layered viscoelastic half-space. Soil Dynamics and Earthquake Engineering, 1993, 12(3): 159–171

[16] Luco J E, De Barros F C P. Seismic response of a cylindrical shell embedded in a layered viscoelastic half-space. I: Formulation. Earthquake Engineering & Structural Dynamics, 1994, 23(5): 553–567

[17] Luco J E, De Barros F C P. Three-dimensional response of a layered cylindrical valley embedded in a layered half-space. Earthquake Engineering & Structural Dynamics, 1995, 24(1): 109–125

[18] Zhang L, Chopra A K. Three-dimensional analysis of spatially varying ground motions around a uniform canyon in a homogeneous half-space. Earthquake Engineering & Structural Dynamics, 1991, 20(10): 911–926

[19] Stamos A A, Beskos D E. Dynamic analysis of large 3-D underground structures by the BEM. Earthquake Engineering & Structural Dynamics, 1995, 24(6): 917–934

[20] Stamos A A, Beskos D E. 3-D seismic response analysis of long lined tunnels in half-space. Soil Dynamics and Earthquake Engineering, 1996, 15(2): 111–118

[21] Ba Z, Fu J, Liu Y, Lee V W, Wang Y. Scattering of elastic spherical P, SV, and SH waves by three-dimensional hill in a layered half-space. Soil Dynamics and Earthquake Engineering, 2021, 147: 106545

[22] Liu J, Du Y, Du X, Wang Z, Wu J. 3D viscous-spring artificial boundary in time domain. Earthquake Engineering and Engineering Vibration, 2006, 5(1): 93–102

[23] Liu J, Gu Y, Li B, Wang Y. An efficient method for the dynamic interaction of open structure-foundation systems. Frontiers of Structural and Civil Engineering, 2007, 1(3): 340–345

[24] Du J, Lin G. Improved numerical method for time domain dynamic structure-foundation interaction analysis based on scaled boundary finite element method. Frontiers of Structural and Civil Engineering, 2008, 2(4): 336–342

[25] Hatzigeorgiou G D, Beskos D E. Soil−structure interaction effects on seismic inelastic analysis of 3-D tunnels. Soil Dynamics and Earthquake Engineering, 2010, 30(9): 851–861

[26] Li P, Song E. Three-dimensional numerical analysis for the longitudinal seismic response of tunnels under an asynchronous wave input. Computers and Geotechnics, 2015, 63: 229–243

[27] Liu J, Tan H, Bao X, Wang D, Li S. Seismic wave input method for three-dimensional soil−structure dynamic interaction analysis based on the substructure of artificial boundaries. Earthquake Engineering and Engineering Vibration, 2019, 18(4): 747–758

[28] Liu Z, Zhang M, Huang L, Zhang H. Broadband seismic isolation effect of three-dimensional seismic metamaterials in a semi-infinite foundation: Modelled by fast multipole indirect boundary element method. Engineering Analysis with Boundary Elements, 2023, 154: 7–20

[29] Yin C, Li W H, Wang W. Evaluation of ground motion amplification effects in slope topography induced by the arbitrary directions of seismic waves. Energies, 2021, 14(20): 6744

[30] Qiao H, Dai Z, Du X, Wang C, Long P, Jiao C. Parametric study of topographic effect on train-bridge interaction of a continuous rigid frame bridge during earthquakes. Bulletin of Earthquake Engineering, 2023, 21(1): 125–149

[31] UnglessR F. Infinite finite element. Dissertation for the Doctoral Degree. Vancouver: University of British Columbia, 1973

[32] Bettess P. Infinite elements. International Journal for Numerical Methods in Engineering, 1977, 11(1): 53–64

[33] Zhao C, Valliappan S, Wang Y C. A numerical model for wave scattering problems in infinite media due to P- and SV-wave incidences. International Journal for Numerical Methods in Engineering, 1992, 33(8): 1661–1682

[34] Zhao C, Valliappan S. Seismic wave scattering effects under different canyon topographic and geological conditions. Soil Dynamics and Earthquake Engineering, 1993, 12(3): 129–143

[35] Zhao C, Valliappan S. Incident P and SV wave scattering effects under different canyon topographic and geological conditions. International Journal for Numerical and Analytical Methods in Geomechanics, 1993, 17(2): 73–94

[36] Yang Y B, Kuo S R, Hung H H. Frequency-independent infinite element for analyzing semi-infinite problems. International Journal for Numerical Methods in Engineering, 1996, 39(20): 3553–3569

[37] Yun C B, Kim D K, Kim J M. Analytical frequency-dependent infinite elements for soil−structure interaction analysis in two-dimensional medium. Engineering Structures, 2000, 22(3): 258–271

[38] Kim D K, Yun C B. Time-domain soil−structure interaction analysis in two-dimensional medium based on analytical frequency-dependent infinite elements. International Journal for Numerical Methods in Engineering, 2000, 47(7): 1241–1261

[39] Yang Y B, Hung H H, Lin K C, Cheng K W. Dynamic response of elastic half-space with cavity subjected to P and SV waves by finite/infinite element approach. International Journal of Structural Stability and Dynamics, 2015, 15(7): 1540009

[40] Lin S Y, Hung H H, Yang J P, Yang Y B. Seismic analysis of twin tunnels by a finite/infinite element approach. International Journal of Geomechanics, 2017, 17(9): 04017060

[41] Yang Y B, Zhou Z, Zhang X, Wang X. Soil seismic analysis for 2D oblique incident waves using exact free-field responses by frequency-based finite/infinite element method. Frontiers of Structural and Civil Engineering, 2022, 16(12): 1530–1551

[42] Yang Y B, Zhang X, Zhou Z, Wang X, Li Z. Seismic analysis of a half-space containing a water- filled valley under 2D oblique incident waves by finite-infinite element method. Soil Dynamics and Earthquake Engineering, 2023, 169: 107872

[43] Yang Y B, Hung H H A. 2. 5D finite/infinite element approach for modelling visco-elastic bodies subjected to moving loads. International Journal for Numerical Methods in Engineering, 2001, 51(11): 1317–1336

[44] Yang Y B, Hung H H, Chang D W. Train-induced wave propagation in layered soils using finite/infinite element simulation. Soil Dynamics and Earthquake Engineering, 2003, 23(4): 263–278

[45] Hung H H, Yang Y B, Chang D W. Wave barriers for reduction of train-induced vibrations in soils. Journal of Geotechnical and Geoenvironmental Engineering, 2004, 130(12): 1283–1291

[46] Lin K C, Hung H H, Yang J P, Yang Y B. Seismic analysis of underground tunnels by the 2.5D finite/infinite element approach. Soil Dynamics and Earthquake Engineering, 2016, 85: 31–43

[47] Yang Y B, Ge P B, Li Q M, Liang X J, Wu Y T. 2.5D vibration of railway-side buildings mitigated by open or infilled trenches considering rail irregularity. Soil Dynamics and Earthquake Engineering, 2018, 106: 204–214

[48] Yang Y B, Liu S J, Li Q M, Ge P B. Stress waves in half-space due to moving train loads by 2.5D finite/infinite element approach. Soil Dynamics and Earthquake Engineering, 2019, 125: 105714

[49] Yang Y B, Li P L, Chen W, Li J, Wu Y T. 2.5D formulation and analysis of a half-space subjected to internal loads moving at sub-and super-critical speeds. Soil Dynamics and Earthquake Engineering, 2021, 142: 106550

[50] Xie W, Gao G, Song J, Wang Y. Ground vibration analysis under combined seismic and high-speed train loads. Underground Space, 2022, 7(3): 363–379

[51] Romero A, Galvín P, António J, Domínguez J, Tadeu A. Modelling of acoustic and elastic wave propagation from underground structures using a 2.5D BEM-FEM approach. Engineering Analysis with Boundary Elements, 2017, 76: 26–39

[52] He C, Zhou S, Guo P, Di H, Zhang X. Modelling of ground vibration from tunnels in a poroelastic half-space using a 2.5-D FE-BE formulation. Tunnelling and Underground Space Technology, 2018, 82: 211–221

[53] Noori B, Arcos R, Clot A, Romeu J. A method based on 3D stiffness matrices in Cartesian coordinates for computation of 2.5D elastodynamic Green’s functions of layered half-spaces. Soil Dynamics and Earthquake Engineering, 2018, 114: 154–158

[54] Ba Z, Fu Z, Liu Z, Sang Q A. 2.5D IBEM to investigate the 3D seismic response of 2D topographies in a multi-layered transversely isotropic half-space. Engineering Analysis with Boundary Elements, 2020, 113: 382–401

[55] Kausel E, de Oliveira Barbosa J M. Analysis of embankment underlain by elastic half-space: 2.5D model with paralongitudinal approximations to the half-space. Soil Dynamics and Earthquake Engineering, 2022, 155: 107090

[56] Zhu J, Li X J, Liang J W. 2.5D FE-BE modelling of dynamic responses of segmented tunnels subjected to obliquely incident seismic waves. Soil Dynamics and Earthquake Engineering, 2022, 163: 107564

[57] Zhang Q, Zhao M, Huang J Q, Du X L, Zhang G L A. 2.5D finite element method combined with zigzag-paraxial boundary for long tunnel under obliquely incident seismic wave. Applied Sciences, 2023, 13(9): 5743

[58] Zhao C, Valliappan S. An efficient wave input procedure for infinite media. Communications in Numerical Methods in Engineering, 1993, 9(5): 407–415

RIGHTS & PERMISSIONS

Higher Education Press