Effect of a large-scale three-dimensional sedimentary basin on Rayleigh wave propagation by using spectral element method combined with frequency-wavenumber method

Abstract

Site effects study has always been a key research topic in earthquake engineering. This study proposes a hybrid method to analyze large-scale three-dimensional sedimentary basin under Rayleigh (R) wave incidence. The proposed hybrid method includes two steps: 1) calculate the free field responses of layered sites subjected to R-wave using the frequency-wavenumber method; 2) Simulate the local site region using spectral element method with the equivalent forces input computed from the free field responses. A comprehensive verification study is conducted demonstrating the accuracy of this method. To investigate the effect of sedimentary basin on R-wave propagation, a parametric study is performed on the medium impedance contrast ratio of sedimentary basins and the incident seismic wave predominant frequency, revealing the scattering patterns of sedimentary basins under R-wave incidence. Finally, a practical case of the Wudu Basin in the Tibetan Plateau region of China is simulated. Results indicate significant amplification of R-wave by sedimentary basin, and the proposed hybrid method could serve as a reliable and efficient approach for large-scale R-wave propagation simulation.

Graphical abstract

Keywords

Rayleigh wave / hybrid method / basin effect / wave propagation / spectral element method

1 Introduction

Extensive seismic damage surveys, strong motion records, and theoretical studies indicate that sedimentary basins significantly impact seismic wave propagation, leading to local amplification or attenuation effects [1–5]. A typical example, the Mexico City, located over 400 km away from epicenter, suffered severe damage due to sedimentary basin effect in the 1985 Mexico City earthquake. Seismic waves repeatedly reflected within the basin, causing acceleration amplification, prolonging the duration to about 180 s, with a peak acceleration of 0.18g [6]. In the 1995 Kobe earthquake, severe damage occurred in a strip area at the Osaka basin edge, which suffered evidently heavier damage than other nearby places [7]. Surface waves are a crucial component of long-period ground motion, especially in sedimentary basins where Rayleigh-wave (R-wave)would be induced by basin [8]. Hence, it is important to study the scattering characteristics of R-wave for reliable seismic damage assessment and proper engineering site selection in complex surface topography.

The impact of surface motion on sedimentary basin is essentially due to its role as a heterogeneous scatterer, causing complex scattering (diffraction) and wave mode conversion. As early as the 1970s, numerous scholars conducted studies on the effects of sedimentary basins. Trifunac [9] derived an accurate analytical solution for the scattering of shear horizontal (SH) waves by a semicircular sedimentary basin situated within a two-dimensional isotropic homogeneous half-space, employing the wave function expansion technique. An analytical solution regarding the scattering of R-wave by a shallow circular sedimentary valley was offered by Todorovska and Lee [10], who also employed the wave function expansion approach. Nonetheless, the limitations of analytical methods arise from their foundational assumptions, and only simplified terrains in half-space could be accounted. As a result, to accurately simulate more intricate and realistic scenarios, numerical methods have been developed over the last few decades [11].

The main numerical methods include: boundary element method (BEM), finite difference method (FDM), finite element method (FEM), and spectral element method (SEM). Wherein, BEM gains popularity for its merits of dimension reduction, accommodation of radiation boundary conditions, and accurate characterization of the wave field. Kawase and Aki [12] employed the discrete wave number BEM to examine the response of semicircular canyons and deep basins to incident R-wave in two-dimensional settings. Semblat et al. [13] mplemented the BEM to assess the impact of incident SH waves on alluvial basins in Caracas, Venezuela. Sánchez-Sesma et al. [14,15] utilized the indirect BEM to explore the interaction of R-wave with alluvial valleys across both two-dimensional and three-dimensional models. It is noteworthy that various BEMs may be difficult to solve in complex terrain problems and have limited guidance for practical engineering. To investigate site effects in complex medium, many researchers employed FDM and FEM to simulate the seismic wave field of local sites. Frankel and Vidale [16] used the FDM to simulate the seismic wave propagation in a three-dimensional (3D) model of the Santa Clara sedimentary basin located in San Jose, California. Opršal and Zahradnik [17] utilized the FDM to investigate the seismic response of the two-dimensional (2D) Volvi Lake sedimentary basin under the influence of compressive (P-) wave and shear vertical (SV-)waves. The seismic characteristics of near-surface topography influenced by R-wave was examined by Wang et al. [18] using the FDM. Zhang et al. [19] used the FEM to analyze the seismic response of a 3D Wudu Basin model under incident SV-wave. Qiang et al. [20] employed the FEM to perform a detailed study on the basin amplification effect of an ideal 2D basin model. While FDM and FEM can be utilized for simulating seismic wave fields in complex heterogeneities, it is important to note that a significant number of elements would be established when solving large-scale broadband 3D wave fields, which leads to a substantial consumption of computational resources.

Actually, since 1998 Komatitsch and Vilotte first introduced the SEM into seismology [21], the SEM has widely been utilized to address the issue of large-scale ground motion simulation [22–28]. The SEM, which combines the geometrical flexibility of the FEM with the accuracy of the pseudo spectral method, is a powerful tool to handle broadband wave propagation problem due to its relatively low computation demand. It is noteworthy that existing seismic simulations based on the SEM predominantly address source problems, typically considering the ground motion of local sites caused by fault rupture within the computational area. However, there are limitations in simulating plane wave propagation, which is important in the study of site effects for engineering seismic design. Tong et al. [29,30] used the propagator stiffness matrix method and SPECFEM3D to achieve plane P- and SV-waves field simulation in layered sites. Liang et al. [31] employed the frequency-wavenumber (FK) and SEM to simulate the plane P-, SV-, and SH-waves propagation in localized regions. To improve the computational efficiency of wavefield simulation in complex localized regions, a number of researchers have turned to the development of hybrid simulation techniques [32–34]. However, few studies are available on SEM simulations of 3D irregular sites in layered half-space under arbitrarily incident R-wave.

For shallow earthquakes, body waves have significantly attenuated in areas far from the epicenter, and surface waves may play a major role. Among them, R-wave propagate far with large amplitudes, which may result in soil liquefaction and severe structure damage. Miller et al. [35] shows the percentage of three types of elastic seismic waves in actual measured ground motions: surface waves account for 67.3%, shear waves for 25.8%, and compressional waves for 6.9%. In earthquake engineering, it is generally accepted that the destructive force of ground motion is primarily caused by shear waves and surface waves [36]. Additionally, the dispersive characteristics of R-wave in a layered half-space are an important issue in research. Dispersion is a phenomenon where the propagation speed of a wave varies with its frequency, causing spatial variations in ground motion. R-wave exhibit slow energy attenuation and long propagation distances, making them a major cause of seismic hazards in the mid-to-far field. Additionally, complex local topography alters the spatial distribution of dynamic stiffness, thereby exerting influence on R-wave propagation and seismic responses. At present, many studies focused on P- and SV-waves [37,38], and R waves are rarely studied. The equivalent nodal force method was employed by Zhang et al. [39] to model the R-wave propagation in layered locations via the FEM. Nguyen et al. [40] combined thin-layer method with domain reduction method to achieve FEM simulations in layered sites for optimizing meso-scale (sub-km) simulations. These hybrid methods are able to simulate R-wave propagation in layered sites, but there are certain limitations when dealing with real surfaces in large-scale regions.

In this study, we propose a robust method to analyze 3D local sites subjected to R-wave incidence in layered soils. The method consists of two steps: first, the FK method is utilized to solve the free field response. Subsequently, equivalent forces calculated from the free field responses are applied to the boundary of the local region, and the propagation of R-wave in the local region is simulated with a fine SEM model using open-source software SPECFEM3D. The article is organized as follows. Section 2 outlines the overall methodology and provides a brief overview of the formulations related to FK, SEM, and equivalent forces. Section 3 provides comprehensive verification, and Section 4 conducts a parametric study of the sedimentary basin model under incident R-wave. In Section 5, the propagation of R-wave in the Wudu Basin area of the Qinghai-Tibet Plateau in China is simulated. Finally, Section 6 summarizes the main conclusions of this study.

2 Formulations of the proposed method

In this study, we integrated the FK method, capable of solving free-field R-wave responses, into SPECFEM3D. The FK-SE hybrid method combines the advantages of the FK method and SEM. Advantage of the proposed method mainly lies in its superior efficiency for large scale ground motion simulations due to its frequency-wavenumber domain solutions of the FK method. Although the FK method is limited to layered half-space studies, but it still serves well for certain study purposes, especially for local site effect study which are not concerned about topography in a large-scale domain.

Fig.1 illustrates the FK-SE hybrid method. Fig.1(a) illustrates the model of a layered site with the 3D local topography (such as canyon, mountain, basin) in our study. The model in Fig.1(a) can be divided into two models as shown in Fig.1(b) and Fig.1(c). Fig.1(b) shows a model of a layered half-space with the local region. The layered half-space consists of N layers with thickness h_i for ith layer. Note that N could be an arbitrary integer number. The layered medium is assumed to be homogeneous, isotropic, and viscoelastic, which is described by two Lamé’s constants of v and G, mass density ρ, and damping ratio ζ. R-wave, with an incident angle φ with the y-axis, is propagated from the far-field to the local regions. The SE model is shown in Fig.1(c). The SPECFEM3D is used to simulate the wave field of local region, and the Stacey absorbing boundary condition is adopted to absorb the outgoing waves at boundary.

Specific implementation steps of the proposed method are as follows.

1) Obtain the coordinates of all boundary Gauss-Lobatto-Legendre (GLL) nodes where the input force will be applied.

2) Determine the incident wave time step size dt, and identify the frequency points (w_i) to be used in FK simulation.

3) Establish the global stiffness matrix [K_P-SV] for R-wave in the layered site (free field) at each frequency point. Solve for the dispersion curves c_a(ω) of each mode to establish the relationship between phase velocity c_a and frequency ω.

4) Substitute the corresponding dispersion curves c_a(ω) into the global stiffness matrix [K_P-SV] to compute the corresponding eigenvectors. This step obtains the in-plane displacement amplitudes for each layer in the frequency-wavenumber domain.

5) Utilize the amplitude coefficients of upward and downward waves in the layers, and calculate the displacement and stress amplitudes at the boundary GLL nodes.

6) Transform the assigned displacement and stress amplitudes using Eq. (1) to obtain the in-plane and out-of-plane displacements and stresses in the frequency domain.

7) Use inverse FFT to obtain the displacement and stress time histories.

8) Calculate the equivalent forces {F} with simulated free field response at the boundary GLL nodes, and perform SEM analysis with the equivalent forces as input.

(1)\left\{\begin{array}{l}u\\ v\\ w\end{array}\right\}=\left\{\begin{array}{ccc}\sin \phi & \cos \phi & 0\\ \cos \phi & \sin \phi & 0\\ 0& 0& 1\end{array}\right\}\left\{\begin{array}{l}u\\ 0\\ w\end{array}\right\} .

2.1 The frequency-wavenumber method

The FK method is a semi-analytical technique designed for solving seismic wavefields in layered half-spaces. It directly solves wave equations, avoiding cumulative errors and ensuring high accuracy. The FK method uses frequency-wavenumber domain Green’s functions to resolve broadband frequency seismic wave propagation, making it suitable for engineering applications. It supports multi-scale calculations, accommodating simulations ranging from kilometer-scale crustal formations to meter-scale geotechnical layers. This special feature enables detailed modeling of near-surface velocity structure and its influence on seismic ground motion.

Several FK methods have been developed over the past few decades, including the discrete wave-number integration method [41], propagator matrix method [42,43], reflection transmission coefficient method [44], and stiffness matrix method [45,46], which has shown good potential for various applications [47].

Notably, the authors have developed a global exact stiffness matrix method based on Wolf’s theory [48] to evaluate ground motion characteristics [49–52]. These advancements offer a robust foundation for solving free field motions in layered sites subjected to incident plane waves, utilizing the FK method in this work.

The seismic motion in a 3D horizontal layer with thickness h is illustrated in Fig.2. The amplitude of the displacements and stresses in the frequency-wavenumber domain along the depth of z-axis can be expressed as follows:

(2)\{\begin{array}{rl}& U(\mathrm{z};\omega ,k)=\frac{V_{\mathrm{P}}^{\ast }k}{\omega }\left(A_{\mathrm{P}}{\mathrm{e}}^{\mathrm{i}k{\phi }_1\mathrm{z}}+B_{\mathrm{P}}{\mathrm{e}}^{-\mathrm{i}k{\phi }_1\mathrm{z}}\right)-\frac{V_{\mathrm{S}}^{\ast }k}{\omega }{\phi }_2\left(A_{\mathrm{S}\mathrm{V}}{\mathrm{e}}^{\mathrm{i}k{\phi }_2\mathrm{z}}-B_{\mathrm{S}\mathrm{V}}{\mathrm{e}}^{-\mathrm{i}k{\phi }_2\mathrm{z}}\right),\\ & W(\mathrm{z};\omega ,k)=-\frac{V_{\mathrm{P}}^{\ast }k}{\omega }{\phi }_1\left(A_{\mathrm{P}}{\mathrm{e}}^{\mathrm{i}k{\phi }_1\mathrm{z}}-B_{\mathrm{P}}{\mathrm{e}}^{-\mathrm{i}k{\phi }_1\mathrm{z}}\right)-\frac{V_{\mathrm{S}}^{\ast }k}{\omega }\left(A_{\mathrm{S}\mathrm{V}}{\mathrm{e}}^{\mathrm{i}k{\phi }_2\mathrm{z}}+B_{\mathrm{S}\mathrm{V}}{\mathrm{e}}^{-\mathrm{i}k{\phi }_2\mathrm{z}}\right),\end{array}

(3)\{\begin{array}{rl}& {\tau }_{\mathrm{z}x}(\mathrm{z};\omega ,k)=\frac{k^2{G}^{\ast }}{\omega }\left[2\mathrm{i}{V}_{\mathrm{P}}^{\ast }{\phi }_1\left(A_{\mathrm{P}}{\mathrm{e}}^{\mathrm{i}k{\phi }_1\mathrm{z}}-B_{\mathrm{P}}{\mathrm{e}}^{-\mathrm{i}k{\phi }_1\mathrm{z}}\right)+V_{\mathrm{S}}^{\ast }\left(1-{\phi }_2^2\right)\left(A_{\mathrm{S}\mathrm{V}}{\mathrm{e}}^{\mathrm{i}k{\phi }_2\mathrm{z}}+B_{\mathrm{S}\mathrm{V}}{\mathrm{e}}^{-\mathrm{i}k{\phi }_2\mathrm{z}}\right)\right],\\ & {\tau }_{\mathrm{z}\mathrm{z}}(\mathrm{z};\omega ,k)=\frac{k^2{G}^{\ast }}{\omega }\left[\mathrm{i}{V}_{\mathrm{P}}^{\ast }\left(1-{\phi }_2^2\right)\left(A_{\mathrm{P}}{\mathrm{e}}^{\mathrm{i}k{\phi }_1\mathrm{z}}+B_{\mathrm{P}}{\mathrm{e}}^{-\mathrm{i}k{\phi }_1\mathrm{z}}\right)-2\mathrm{i}{V}_{\mathrm{S}}^{\ast }\left(A_{\mathrm{S}\mathrm{V}}{\mathrm{e}}^{\mathrm{i}k{\phi }_2\mathrm{z}}+B_{\mathrm{S}\mathrm{V}}{\mathrm{e}}^{-\mathrm{i}k{\phi }_2\mathrm{z}}\right)\right],\end{array}

where U and W represent the displacement components in x- and z-directions, respectively, and τ denotes the stress components. The wave amplitudes of the up-going waves are denoted as A_P and A_SV, while B_P and B_SV correspond to the wave amplitudes of the down-going waves. The compressional (P-) wave velocity V_P and shear (S-) wave velocity V_S are defined as V_{\mathrm{P}}=\sqrt{G/\rho } and V_{\mathrm{S}}= \sqrt{(v+2G)/\rho }, respectively. With material damping ξ, the velocities are adjusted to V_{\mathrm{P}}^{\ast }= V_P (1 + 2iξ) and V_{\mathrm{S}}^{\ast }= V_S(1 + 2iξ), where G is replaced by G^* = G(1 + 2iξ). k denotes the radial wavenumber, and ω denotes the circular frequency while i = \sqrt{-1} is the imaginary unit. The two parameters, φ₁, φ₂, are defined as follows:

(4){\phi }_1=\sqrt{\left[\frac{{\omega }^2}{{\left(V_{\mathrm{P}}^{\ast }\right)}^2}-1\right]}, {\phi }_2=\sqrt{\left[\frac{{\omega }^2}{{\left(V_{\mathrm{S}}^{\ast }\right)}^2}-1\right]}.

Based on Eqs. (2) and (3), the amplitudes of displacement and stress at the top (subscript “1”) and bottom (subscript “2”) surfaces of the layer with respect to the wave amplitudes can be expressed as

(5){\left\{U_1,\mathrm{i}{W}_1,U_2,\mathrm{i}{W}_2\right\}}^{\mathrm{T}}= \left[{\mathit{D}}_{\mathrm{P}-\mathrm{S}\mathrm{V}}\right]{\left\{A_{\mathrm{P}},B_{\mathrm{P}},A_{\mathrm{S}\mathrm{V}},B_{\mathrm{S}\mathrm{V}}\right\}}^{\mathrm{T}},

(6){\left\{-{\tau }_{\mathrm{z}x1},-\mathrm{i}{\tau }_{\mathrm{z}\mathrm{z}1},{\tau }_{\mathrm{z}x2},\mathrm{i}{\tau }_{\mathrm{z}\mathrm{z}2}\right\}}^{\mathrm{T}}=\left[{\mathit{T}}_{\mathrm{P}-\mathrm{s}\mathrm{V}}\right]{\left\{A_{\mathrm{P},}{B}_{\mathrm{P},}{A}_{\mathrm{S}\mathrm{V}},B_{\mathrm{S}\mathrm{V}}\right\}}^{\mathrm{T}},

where D and T denote the displacement and stress coefficient matrix, respectively. The expression of coefficient matrix D and T can be seen in Appendix A in Electronic Supplementary Materials. Then, using Eqs. (5) and (6), the dynamic stiffness matrix of a layer is established and the relationship of stress and displacement of a layer can be expressed as

(7)\begin{array}{rl}& \{-{\tau }_{\mathrm{z}x1},-\mathrm{i}{\tau }_{\mathrm{z}\mathrm{z}1},{\tau }_{\mathrm{z}x2},\mathrm{i}{\tau }_{\mathrm{z}\mathrm{z}2}{\}}^{\mathrm{T}}\\ =& [{\mathit{T}}_{\mathrm{P}-\mathrm{S}\mathrm{V}}{\left[{\mathit{D}}_{\mathrm{P}-\mathrm{S}\mathrm{V}}\right]}^{-1}]\{U_1,\mathrm{i}{W}_1,U_2,\mathrm{i}{W}_2{\}}^{\mathrm{T}}.\end{array}

It is worth noting that the matrices D_P-SV and T_P-SV involved in Eq. (7) contain divergent exponential terms in the upward wave amplitude as shown in Eq. (8), which could cause the matrix to be singular.

(8)\exp \left(\mathrm{i}h\sqrt{\frac{{\omega }^2}{(c_{\mathrm{P}}^{\ast }{)}^2}-k}\right),\exp \left(\mathrm{i}h\sqrt{\frac{{\omega }^2}{(c_{\mathrm{S}}^{\ast }{)}^2}-k}\right).

Therefore, the divergent exponential terms in the above matrices are uniformly extracted into the wave amplitude vector, eliminating all divergent exponential terms in matrices D and T. Thus, Eq. (7) is revised to Eq. (9):

(9)\begin{array}{rl}& \{-{\tau }_{\mathrm{z}x1},-\mathrm{i}{\tau }_{\mathrm{z}\mathrm{z}1},{\tau }_{\mathrm{z}x2},\mathrm{i}{\tau }_{\mathrm{z}\mathrm{z}2}{\}}^{\mathrm{T}}\\ =& [{\mathit{T}}_{\mathrm{P}-\mathrm{S}\mathrm{V}}^{\prime }\left[{\mathit{D}}_{\mathrm{P}-\mathrm{S}\mathrm{V}}^{\prime }{]}^{-1}\right]\{U_1,\mathrm{i}{W}_1,U_2,\mathrm{i}{W}_2{\}}^{\mathrm{T}},\end{array}

where the superscript ‘r’ indicates the corrections made to the matrix. The amplitudes of external loads are specified as P₁ = − τ_zx1, Q₁ = − τ_zy1, iR₁ = − iτ_zz1 (top surface) and P₂ = τ_zx2, Q₂ = τ_zy2, R₂ = iτ_zz2 (bottom surface). Let {\mathit{K}}_{\mathrm{P}-\mathrm{S}\mathrm{V}}^{l_i} represent the in-plane dynamic stiffness matrix of layer l_i. By reorganizing Eq. (9), the dynamic stiffness matrix of the layer element is obtained.

(10)\{P_1,\mathrm{i}{R}_1, P_2,\mathrm{i}{R}_2{\}}^{\mathrm{T}}=[{\mathit{K}}_{\mathrm{P}-\mathrm{S}\mathrm{V}}^{l_i}{]}_{4\times 4}\{U_1,\mathrm{i}{W}_1,U_2,\mathrm{i}{W}_2{\}}^{\mathrm{T}}.

Similarly, the radiation condition for the half-space contains the down-going wave only. The dynamic stiffness matrix of a half-space (subscript “0”) can be determined as Eq. (11):

(11)\{P_0,\mathrm{i}{R}_0{\}}^{\mathrm{T}}=[ {\mathit{K}}_{\mathrm{P}-\mathrm{S}\mathrm{V}0}{]}_{2\times 2}\{U_0 ,\mathrm{i}{W}_0{\}}^{\mathrm{T}}.

According to the continuity condition at the interfaces between adjacent layers, the global stiffness matrix for the entire layered model is assembled by combining the stiffness matrices from each layer as described in Eqs. (10) and (11). The assembly process is depicted in Fig.2. Then, the dynamic equilibrium equation of the layered model can be written as Eq. (12):

(12)\begin{array}{rl}& \{P_1,\mathrm{i}{R}_1,\dots ,P_i,\mathrm{i}{R}_i,\dots ,P_{n+1},\mathrm{i}{R}_{n+1}{\}}^{\mathrm{T}}\\ =& [{\mathit{K}}_{\mathrm{P}-\mathrm{S}\mathrm{V}}{]}_{2N+2}\{U_1,\mathrm{i}{W}_1,\dots ,U_i,\mathrm{i}{W}_i,\dots ,U_{n+1},\mathrm{i}{W}_{n+1}{\}}^{\mathrm{T}}.\end{array}

For R-wave, the load vector is zero, that is, in Eq. (12), {P₁,iR₁,...,P_i,iR_i,...,P_n+1,iR_n+1}^T = 0. By using Newton’s iterative method to set the determinant of the stiffness matrix |K_P-SV| = 0, the dispersion equation is obtained, which determines the relationship between phase velocity c_a and frequency ω for different modes. The corresponding c_a(ω) is obtained, which is the dispersion curve for that mode. After obtaining the dispersion curves of each mode, substitute the pairs of (ω, c_a(ω)) into the global stiffness matrix [K_P-SV] to obtain the corresponding eigenvectors, yielding the displacement amplitude vectors for each layer {U₁,iW₁,...,U_i,iW_i,...,U_n+1,iW_n+1}^T. Then, for any layer, the amplitude coefficients of the upward and downward waves in that soil layer can be obtained using Eq. (5), and the in-plane displacement and stress amplitudes at any location within the layer can be determined. Through the coordinate transformation in Eq. (1), the in-plane and out-of-plane displacement and stress amplitudes are obtained, representing the free wave field in the layered site under incident R-wave. Calculate the velocity amplitudes at each point using the displacement amplitudes. Using inverse FFT to obtain the time histories of the velocities in three directions {\dot{u}}_x^{\mathrm{F}\mathrm{K}}, {\dot{u}}_y^{\mathrm{F}\mathrm{K}}, {\dot{u}}_{\mathrm{z}}^{\mathrm{F}\mathrm{K}} and the six stresses {\sigma }_{xx}^{\mathrm{F}\mathrm{K}}, {\sigma }_{yy}^{\mathrm{F}\mathrm{K}}, {\sigma }_{\mathrm{z}\mathrm{z}}^{\mathrm{F}\mathrm{K}}, {\sigma }_{xy}^{\mathrm{F}\mathrm{K}}, {\sigma }_{x\mathrm{z}}^{\mathrm{F}\mathrm{K}}, {\sigma }_{y\mathrm{z}}^{\mathrm{F}\mathrm{K}}.

2.2 Spectral element method

SPECFEM3D, which is an open-source SEM code, is designed for 3D wave field simulations [53,54]. It is capable of simulating both global and local seismic wave propagation processes and perform full-waveform imaging. High-precision parallel computation is implemented using MPICH, and GPU computation is also available.

The SEM uses Lagrange polynomials as basic functions and GLL collocation points for interpolation within each element, facilitating diagonalization of the mass matrix and enabling efficient parallel processing. However, the grid size, integration method and element type all affect the stability and accuracy of wave field calculation by SEM. Therefore, the Courant−Friedrichsen−Lewy (CFL) condition should be satisfied to impose constraints on the time-step size Δt as shown in Eq. (13):

(13)\mathrm{\Delta }t⩽C^{\ast }{\left(\mathrm{\Delta }h/V_{\mathrm{P}}\right)}_{\mathrm{m}\mathrm{i}\mathrm{n}},

where Δh denotes the minimum of distances between adjacent GLL nodes and C is the CFL constant. Meanwhile, the element size should satisfy that a wavelength is represented by at least five GLL points, which can be expressed as:

(14)d⩽(N /5 ){\lambda }_{\mathrm{m}\mathrm{i}\mathrm{n}},

where N is the order of the Lagrange polynomial, and λ_min is the shortest wavelength in the propagating medium which corresponds to the maximum frequency of interest.

2.3 Seismic input from frequency-wavenumber method to spectral element method

The free field motions of the layered half-space as the form of the equivalent forces are applied on the boundary of SEM model. Diffracted waves are absorbed using a simple and approximate absorbing boundary condition based on Clayton and Engquist [55]. The scattered waves propagate outward from the SEM model, and the wave field U_total at the boundary consists of two components: the incident wave field U_FK and the scattered wave field U_SC within the region, as shown in Eq. (15). The incident wave field are calculated by the FK method.

(15)U_{\mathrm{t}\mathrm{o}\mathrm{t}}=U_{\mathrm{F}\mathrm{K}}+U_{\mathrm{S}\mathrm{C}}.

At the SEM model boundary, the analytical expression of the incident wave field should be subtracted from the total wave field at the outer edge of the computational domain. And the absorbing boundary condition is then applied to the remaining scattered wave field, which must satisfy Eq. (16).

(16)\left(\mathit{\sigma }-{\mathit{\sigma }}_{\mathrm{F}\mathrm{K}}\right)\cdot \hat{\mathit{n}}=\rho V_{\mathrm{P}}\left[\hat{\mathit{n}}\cdot \left(\dot{\mathit{u}}-{\dot{\mathit{u}}}_{\mathrm{F}\mathrm{K}}\right)\right]\hat{\mathit{n}}+\rho V_{\mathrm{S}}\left[\hat{\mathit{s}}\cdot \left(\dot{\mathit{u}}-{\dot{\mathit{u}}}_{\mathrm{F}\mathrm{K}}\right)\right]\hat{\mathit{s}},

where \hat{\mathit{n}} is the unit outward normal of the boundary, and \hat{\mathit{s}} is the unit vector tangential to the boundary. The equivalent force required at the boundary of the SEM model is obtained by Eq. (17).

(17)\mathit{F}=\rho V_{\mathrm{P}}\left(\hat{\mathit{n}}\cdot {\dot{\mathit{u}}}_{\mathrm{F}\mathrm{K}}\right)\hat{\mathit{n}}+\rho V_{\mathrm{S}}\left(\hat{\mathit{s}}\cdot {\dot{\mathit{u}}}_{\mathrm{F}\mathrm{K}}\right)\hat{\mathit{s}}+{\mathit{\sigma }}_{\mathrm{F}\mathrm{K}}\cdot \hat{\mathit{n}},

where {\mathit{\sigma }}_{\mathrm{F}\mathrm{K}} and {\dot{\mathit{u}}}_{\mathrm{F}\mathrm{K}} are obtained by the FK method in Subsection 2.1.

We integrate the FK method for solving the free field under R-wave incidence discussed in Subsection 2.1 into the SPECFEM3D program. This method notably reduces the time spent on data retrieval from external files stored on the computer’s hard disk. Moreover, given the smaller time step size requirement by SEM, the FK calculation results are interpolated during computation, thereby reducing computer memory usage. This is more conducive to solving large-scale site effect problems, providing a powerful tool for local site effect analysis in earthquake engineering. The proposed method is versatile and applicable to any propagating mode of R-wave, including the fundamental mode and other modes. In this work, we focus on the fundamental mode of R-wave propagation, which is commonly predominant.

3 Method verification

Three numerical models are employed to verify the accuracy of the proposed method through time-domain and frequency-domain results. Model 1 is a horizontally layered site, Model 2 is a hemispherical sedimentary valley in 3D half-space, and Model 3 is an irregular alluvial valley in 3D half-space.

3.1 Case 1: A horizontally layered site

Fig.3 illustrates a three-layer horizontally layered site used to validate the accuracy of the time-domain results. Thickness of the first two layers is 200 and 300 m, respectively. Shear wave velocities for the three layers are 600, 1000, 1500 m/s from top to bottom. Density and Poisson’s ratio are 2000 kg/m³ and 0.3, respectively, for all layers. Size of the SEM model is 850 m × 850 m × 800 m and the R-wave incident angle φ is 30°. A Ricker wavelet is used as the incident wave form, which is defined as follows:

(18)S(t)=S_{\mathrm{a}}\left[2{\left(\pi f_0\left(t-t_0\right)\right)}^2-1\right]{\mathrm{e}}^{-{\left(\pi f_0\left(t-t_0\right)\right)}^2},

where S_a = 0.1 m is the amplitude, f₀ = 1 Hz is a central frequency, t₀ = 2 s is the time when peak value is reached.

The displacements at observation points P₁(0,0), P₂(0,−100), P₃(0,−200), P₄(0,−350), P₅(0,−500), and P₆(0,−600) are computed and compared with those from Korres et al. [34], as illustrated in Fig.4. Symbols u, v, and w denote the displacement time histories in x, y, and z directions. Clearly, the simulated results are in good agreement with results from Korres et al. [34].

3.2 Case 2: A hemispherical sedimentary basin in half space

In Fig.5, a hemispherical sedimentary basin in a homogeneous soil domain subjected to R-wave propagation is investigated. Basin radius a is 50 m. Shear wave velocity of the basin and half-space is 500 and 1000 m/s, respectively. Density of the basin and half-space is 2000 and 3000 kg/m³, respectively. Poisson’s ratio is 0.3, for basin and half-space. Size of the SEM model is 400 m × 400 m × 200 m and the R-wave incident angles φ is 90°. Within the range of [–2,2] for x/a, the observation points are evenly distributed on the ground. The dimensionless frequency η of the sedimentary basin (ωa/πV_S) is selected to be 0.5 and 0.75 for two cases of comparison. A Ricker wavelet with S_a = 1 m, f₀ = 2.5 Hz, and t₀ = 2 s is taken as the incident R-wave. Frequency-domain solutions are derived by applying the Fourier-Hankel transform to the time-domain results obtained using the proposed method. Fig.6 shows the comparisons between results of proposed approach and those from Dravinski and Mossessian [56] using indirect boundary integral equation method.

3.3 Case 3: An irregular alluvial valley in three-dimensional half space

To confirm the effectiveness of the proposed method within the time domain, we refer to study from Sánchez-Sesma and Luzón [15]. The model parameters are h = 0.4/a, b = 0.7a and a = 4000 m, as shown in Fig.7. Shear wave velocity of the basin and half-space is 1000 and 2000 m/s, respectively. Density of the basin and half-space is 1600 and 2000 kg/m³, respectively. Poisson’s ratio of the basin and half-space is 0.35 and 0.25, respectively. The R-wave incident angles φ is 90°.

Displacements u and w in x- and y-axis of the ground surface were collected at 50 evenly spaced observation points, ranging from −7280 to 7680 m. Fig.8 presents the synthetic seismograms generated by the proposed method, demonstrating a good agreement with the results from Sánchez-Sesma and Luzón [15]. Note that local extrema (including troughs and ridges) from results of Sánchez-Sesma and Luzón [15] are added as red dashed lines in Fig.8 for comparison.

4 Case study and analysis

In this section, the proposed method is employed to further reveal the propagation characteristics of R-wave in large-scale 3D sites. Taking a 3D trapezoidal sedimentary basin as an example, the study investigates the impact of different medium impedance contrast ratios (IC) and incident seismic wave predominant frequencies (PF) on R-wave propagation. As shown in Fig.9, the subscripts “b” and “l” denote sediment and soil layer, respectively. The top half-width of the basin is 2 km, the bottom half-width of the basin is 1.5 km, the depth of basin is 0.5 km and the edge angle is 45°. Density of the basin and half-space is 1700 and 2600 kg/m³, respectively. Shear wave velocity of the half-space is 1500 m/s while Shear wave velocity of the basin is varied for different study purposes. In addition, the compressive wave velocity is determined as two times of shear wave velocity.

Size of the SEM model is 10 km × 10 km × 2 km. The total SEM elements number is 355600 and the resolved highest frequency is approximately 10 Hz. Using a computer with an Intel Xeon Platinum 8124 M processor, the computation was performed with 32 processes, totaling 40000 time steps, with each simulation taking 3 h and 5 min.

4.1 Case 1: Effect of medium impedance contrast ratio

Based on the basin model described above, we assume the sedimentary layer in the basin is a homogeneous medium while the other parameters remain unchanged. The V_s,b are selected to be 400, 600, and 1000 m/s, respectively, while corresponding V_p,b are twice of V_s,b. Material density ρ remains unchanged. The corresponding IC (IC = ρ_bV_s,b/ρ_lV_s,b) are 0.175, 0.349, and 0.436, respectively. The Ricker wavelet, following Eq. (18), with S_a = 1 m, f₀ = 1 Hz, and t₀ = 2 s is taken for the incident R-wave. The influence of these IC on the propagation of R-wave will be analyzed.

The seismogram synthetics of horizontal and vertical components and the peak displacement distribution at observation points along the y-axis obtained from different IC models are shown in Fig.10. With the increase of the IC, the wave field becomes relatively simpler from complex, with a significant decrease in amplitude, and the propagation velocity of surface waves gradually increases. Damping of the soil layer causes R-wave to gradually attenuate as they propagate into the basin. This indicates that as IC increases, the transmission capability of energy strengthens and the decay rate of R-wave increases. In addition, the ability of sedimentary layers to capture the energy of R-wave entering the basin significantly decreases.

In Fig.10, the peak values of u_x and u_z at surface observation points outside the basin are approximately 1.0 and 1.5, respectively. Ground motion suppression occurs near the corners of the basin (y/a = −1 and 1), where the amplitude is even smaller than the peak ground motion that outside the basin on the bedrock site. Additionally, the smaller the IC between the interior and exterior medium of the basin, the more apparent the suppression effect on ground motion in this area near the corners of the basin. In the region y/a = [−1,−0.5] and [1,0.5], a significant edge effect of the basin is observed, where the peak ground motion noticeably increases. The maximum peak displacement of u_y is approximately 50%−76.0% of u_x, indicating that the horizontal y component cannot be ignored when considering the seismic response of structures within the basin. Notably, when IC = 0.175, notable amplification occurs in the middle region of the basin (y/a = [−0.5,0.5]). The amplification of u_y is due to the scattering of R-wave at the edge of the basin. As the IC decreases, the amplitude at the basin surface gradually increases. When IC = 0.175, the maximum amplitude reaches about three times of incident wave amplitude for all displacement components, which are evidently larger than other IC cases.

Fig.11 displays the peak displacement distribution of observation points within the surface range x = [−3000 m, 3000 m] and y = [−3000 m, 3000 m] obtained from models with different IC values. Note that the observation points at 50 m equally space intervals are selected for the plot. It is evident that smaller IC values lead to a more complex distribution of displacement peaks within the basin. The amplification effect is particularly pronounced for basins with smaller IC values, and amplification is stronger in the central region of the basin compared to the edges. Overall, both horizontal and vertical component peak values of R-wave decrease with increasing IC. These alternating changes in R-wave intensity may be attributed to interference between incident and reflected R-wave.

For models with smaller IC values, the peak value u_x appears in a wedge-shaped area, u_y appears at the far end of the incident wave, and u_z appears near the corner of the basin edge where the incident wave front is located. For models with larger IC values, the maximum R-wave values are generally located near the basin edge, with the distribution pattern closely related to the basin’s shape. While the peak distribution is more complex for low IC models, the general distribution of larger and smaller peaks remains relatively similar across different models. The position of the maximum peak value at the basin edge tends to shift to the right as the velocity of the soil layer decreases. It is worth noting that when IC = 0.175, the peak displacement in all three directions reaches its maximum value. Specifically, at IC = 0.175, the peak in the u_z direction is 1.46 times that of IC = 0.436, while u_x is 2.23 times and u_y is 2.86 times higher than at IC = 0.436. This significant amplification can lead to considerable damage to structures within the basin. However, seismic response within the basin decreases as IC increases, indicating that IC is not the sole factor influencing seismic response. Other factors affecting wave propagation and interference processes may play important roles in determining the seismic response.

Fig.12 illustrates the distribution of peak ground motion at observation points (−1750,−1750), (−1500,−1500), (−1000,−1000), (0,0), (1000,1000), (1500,1500), and (1750,1750) along depth at 50 m interval. With the decrease of the IC, the distribution of peak displacements along the depth direction becomes more complex, with gradual increases in peak displacements. Beyond a depth of 500 m (outside the basin), the peak displacements attenuate gradually and become less complex. The peak displacements along depth at observation points P₄, P₅, and P₆ (near the far end of the incident wave) exhibit more complexity compared to P₁, P₂, and P₃ (near the front end of the incident wave). This complexity arises due to the intricate interaction between incident and reflected R-wave within the basin, leading to varied propagation patterns of seismic waves. With larger IC values, the transmission capability of energy increases, leading to a higher attenuation rate of R-wave and a weakened ability of the basin to capture the energy entering R-wave. The peak value of u_z along depth near the surface exceeds that at the surface itself. At IC = 0.262, for instance, the peak value of u_z at P₂ is 4.2 times higher than that of the incident wave, highlighting potential risks to underground structures.

4.2 Case 2: Effect of incident seismic wave predominant frequencies

This section studies the effect of incident seismic wave PF. The V_s,b is set to 800 m/s, and the V_p,b is twice the V_s,b, and other material parameters remain unchanged. The incident wave is a Ricker wavelet following Eq. (18), with S_a = 1 m and t₀ = 2 s. Using Ricker wavelets with central frequencies of f_p1 = 0.5 Hz, f_p2 = 1 Hz, and f_p3 = 2 Hz as incident waves, the influence of PF on the seismic characteristics of basin ground motion are analyzed.

As shown in Fig.13, the wave response with different PFs exhibit significant differences in the same basin. PF strongly influences the amplitude of ground acceleration time history in the basin. As PF increases, the displacement time history becomes progressively more complex, with a corresponding increase in the complexity of peak value distributions. Clear reflected R-wave are observed in the time history of u_x within the basin at f_p1 = 0.5 Hz. Displacements in the y direction shows minor or negligible amplification since it is caused by scattered waves. The intense fluctuations in time histories of u_z can also be observed. At f_p2 = 1 Hz, the wave field generated within the basin becomes more complex compared to f_p1 = 0.5 Hz, with a more complex distribution of peak displacements and larger peaks. Further increasing to f_p3 = 2 Hz, the wave field complexity increases, yet the peak displacement decreases. This is because the energy becomes more concentrated with the increase of PF, and the interference between incident and reflected R-wave during propagation leads to significant differences in ground motion within the basin. The fundamental frequency of the basin, estimated by the formula f = V_s/4H, is approximately 0.5 Hz. Interestingly, although resonance is more likely to happen when PF approaches the fundamental frequency of the basin, peak displacement actually occurs at f_p2 = 1 Hz. Although resonance is more likely to occur when the PF is closer to the fundamental frequency of the basin, it may not be pronounced for R-wave propagating along the surface.

Fig.14 and Fig.15 present the response spectra of observation points along the x-axis at y = −1700 m and y = −500 m on the basin surface for three different PFs of incident Ricker wavelets. Results illustrate the significant influence of PF on short-period components (t < 2 s) of the ground response spectrum in the basin, while showing relatively minor effects on the long-period components (t > 2 s). Furthermore, despite incident R-wave having the same PF and amplitude, the response spectra differ at various locations within the basin model. This variability underscores the necessity for detailed studies focusing on different regions of the basin when investigating the effects of R-wave incidence.

For observation points along y = −1700 m, the u_x is larger at the edge of the basin, the u_y has a larger response at the far end of the incident wave, and the u_z is greatly affected by the PF of the incident wave. For observation points along y = −500 m, as the PF increases, the u_x, u_y, and u_z all have larger response spectra at the far end of the incident wave location. Smaller PFs correspond to a wider distribution range of response spectrum values with smaller amplitudes. The maximum response spectrum values consistently occur at f_p2 = 1 Hz, aligning with the patterns observed in the time-domain curves presented in Fig.13. Moreover, with increasing PF, the distribution of response spectrum values becomes increasingly complex. This complexity arises from significant interference effects between incident and reflected R-wave, leading to varied seismic responses within the basin.

5 Illustrative example: Wudu Basin simulation

The Wudu Basin is located in Jiangyou City, Sichuan Province, China. It extends approximately 4 km in the east–west direction and 6 km in the north–south direction, with a relatively thin sediment. Due to the basin effect, the Wudu area suffered severe earthquake damage in the 2008 Wenchuan earthquake compared to the surrounding areas [18]. Therefore, studying the seismic effects of such actual basins has significant theoretical and practical value.

Fig.16 shows the details of Wudu basin topography. The targeted area is from 104.73° to 104.84°E and from 31.82° to 31.93°N. The seismic wave field simulation of R-wave propagation in Wudu basin is performed using the hybrid method proposed in this paper. Size of the SEM model is 10 km × 10 km × 0.6 km. The sediment within the Wudu basin is divided into four layers, and the material parameters of each layer are obtained from borehole data [18], as shown in Tab.1. The digital elevation information for this area is sourced from the National Geospatial Data Cloud. Using meshing software Trelis, a 3D SEM model containing undulating terrain and uneven basin bedrock is established. The mesh size is 200 m outside the basin and 20 m inside the basin. This structured mesh consists of 337000 elements and 378000 nodes, as illustrated in Fig.16. Five observation points, shown as yellow triangles, are selected to analyze the ground motion time history under R-wave incidence. Observation points P₁ and P₅ are located outside the basin, while P₂, P₃, and P₄ are located inside the basin. The computation was performed using a computer with an Intel Xeon Platinum 8124M processor, utilizing 32 processes and totaling 200000 time steps, with the entire model computation taking 14 h and 57 min.

The Ricker wavelet following Eq. (18), with S_a = 1 m, f₀ = 1 Hz, and t₀ = 2 s is taken as the incident wave. Fig.17 shows the displacement time history of five observation points in Wudu Basin. Compared to the incident wave, the duration of ground motion in both the interior and exterior of the basin has increased, with durations outside the basin typically around 2–3 s, while those inside the basin are generally greater than 20 s. The peak displacement values of u_x at observation points P1 and P5, located outside the basin, are small, with u_y being nearly negligible. However, the displacements at observation points inside the basin are significantly amplified, with the maximum u_x reaching 6.25 m and u_y reaching 2.08 m, both occurring at observation point P3. The elevation difference within the basin is around 20 m, and the terrain fluctuates significantly in the eastern region of the basin outside. Therefore, the displacements of observation points P2 and P3 are larger than P4. The thickness of the sediment cover in the area where P2 is located is relatively thin, leading to reduced amplitude, with the maximum u_x reaching 4.29 m and u_y reaching 1.13 m. This indicates that the mountains significantly influence the amplification of R-wave incidence in the basin. Additionally, the thickness of the sediment cover significantly influences the amplification of R-wave.

Fig.18 shows snapshots of the R-wave propagation on the surface between 3.0 s and 6.4 s derived from the Wudu basin. At 3 s, the R-wave reach the edge of the basin. Around 3.8 s, the R-wave reach the interior of the basin, causing slight vibrations. Due to the amplification effect of soft soil, the R-wave amplitude inside the basin is much greater than in the surrounding areas. Since the propagation speed of R-wave inside the basin is much lower than outside, the surface displacement wave field snapshots from 4.6 to 5.2 s show that the arrival of the R-wave front inside the basin is significantly delayed compared to the outside region. Additionally, it is clearly observed that R-wave continue to propagate along the surface in the external area without significant focusing. Around 5.8 to 6.4 s, as the R-wave front passes through local region, the focusing of surface waves can still be clearly observed inside the basin.

Fig.19 shows the spatial distribution of the ground motion amplification factor in the Wudu Basin. The amplification factor is defined as the ratio of the maximum response at each point in the simulation area to the peak value of the input wave. Due to the influence of the 3D basin effect, there is also notable amplitude in the y-direction, i.e., the scattering direction. It is noteworthy that the distribution of the ground motion amplification factor outside the basin is closely related to the topography. Mountain peaks or ridges outside the basin generally have a larger amplification factor, with the maximum reaching 3 which is located in the middle of the basin. In marsh, the amplification factor is generally smaller (around 1.5). Inside the basin, the x-direction motion amplification factor is much larger than that outside the basin, generally greater than 5, with a maximum of 6.3. This demonstrates the amplification effect of the basin on R-wave. Additionally, regions within the basin with thick cover layers have higher amplification factors compared to regions with thin sediment cover.

6 Conclusions

This study presents a combined approach to analyze the effects of R wave incidence on 3D site problems in layered soil conditions. The method involves two steps: 1) solving the free field response of R-wave in layered soil using the FK method to obtain effective forces, and 2) perform SEM simulation with input effective forces using SPECFEM3D. This hybrid method is capable of handling site effect problems for large-scale areas under R-wave incidence since the FK method is highly efficient for large-scale free field problems. In addition, integrating FK method into SPECFEM3D program greatly saves the effort for writing and reading the large amount of equivalent force input files. The robustness and accuracy of the method are verified on free field and local sites subjected to R-wave. Additionally, a parametric study of R-wave propagation in a sedimentary basin model is conducted. Finally, a practical case study of the Wudu Basin, a 10 km × 10 km area in the Qinghai-Tibet Plateau region of China, is analyzed to investigate the topographic effects under R-wave incidence.

Results show that the IC and PF of the basin have significant impact on the ground motion under R-wave incidence. The surface displacement distribution in the basin differs significantly from P- and SV-wave incidences. Main findings of this study are summarized below.

1) As the IC decreases, the ground motion within the basin becomes more complex, with the maximum displacement amplified up to 6.34 times. Conversely, as IC increases, the energy transmission capability is enhanced, and the R-wave attenuation rate increases.

2) The PF also significantly affects the basin response. A higher PF results in a more complex response of the basin surface, although the displacement amplitude does not increase with the PF.

3) The existence of the Wudu Basin notably extends the duration of ground shaking. Due to the interplay between the uneven surface topography and the sedimentary basin, seismic motion amplification within the basin shows considerable variation. Within the basin, the amplification factor for motion in the x-direction, typically exceeding 5 with a peak value of 6.3, is significantly higher than outside. The irregular topography plays an important role in the transmission of R-wave.

4) The hybrid method proposed in this paper can more accurately and comprehensively describe the impact of large-scale sites on R-wave propagation. This confirms the importance of topographic effects in seismic wave field simulation and provides a basis for site effect assessment and seismic response analysis of structures.

Although the proposed hybrid method effectively addresses the site effect problems under R-wave incidence in large-scale 3D complex sites, this paper simplifies the geotechnical parameters of the site. Future research should integrate a nonlinear soil model into SPECFEM3D which allows nonlinear dynamic analysis of near-surface soil under strong ground motion, enhancing the method’s applicability to practical engineering problems.

References

[1] Sánchez-Sesma F J, Francisco J, Velázquez S A. On the seismic response of a dipping layer. Wave Motion, 1987, 9(5): 387–391

[2] Aki K. Local site effects on weak and strong ground motion. Tectonophysics, 1993, 218(1–3): 93–111

[3] Kawase H. The cause of the damage belt in Kobe: “The basin-edge effect,” constructive interference of the direct S-wave with the basin-induced diffracted/Rayleigh waves. Seismological Research Letters, 1996, 67(5): 25–34

[4] Kakoty P, Dyaga S M, Molina Hutt C. Impacts of simulated M9 Cascadia subduction zone earthquakes considering amplifications due to the Georgia sedimentary basin on reinforced concrete shear wall buildings. Earthquake Engineering & Structural Dynamics, 2021, 50(1): 237–256

[5] Marafi N A, Eberhard M O, Berman J W, Wirth E A, Frankel A D. Effects of deep basins on structural collapse during large subduction earthquakes. Earthquake Spectra, 2017, 33(3): 963–997

[6] Anderson J G, Bodin P, Brune J N, Prince J, Singh S K, Quaas R, Onate M. Strong ground motion from the Michoacan, Mexico, earthquake. Science, 1986, 233: 1043–1049

[7] Pitarka A, Irikura K, Iwata T, Sekiguchi H. Three-dimensional simulation of the near-fault ground motion for the 1995 Hyogo-ken Nanbu (Kobe), Japan, earthquake. Bulletin of the Seismological Society of America, 1998, 88(2): 428–440

[8] Brissaud Q, Bowden D C, Tsai V C. Extension of the basin Rayleigh-wave amplification theory to include basin-edge effects. Bulletin of the Seismological Society of America, 2020, 110(3): 1305–132

[9] Trifunac M D. Surface motion of a semi-cylindrical alluvial valley for incident plane SH waves. Bulletin of the Seismological Society of America, 1971, 61(6): 1755–1770

[10] Todorovska M I, Lee V W. A note on scattering of Rayleigh waves by shallow circular canyons: Analytical approach. Bulletin of the Indian Society of Earthquake Technology, 1991, 28(2): 1–16

[11] Poursartip B, Fathi A, Tassoulas J L. Large-scale simulation of seismic wave motion: A review. Soil Dynamics and Earthquake Engineering, 2020, 129: 105909

[12] Kawase H, Aki K. A study on the response of a soft basin for incident S, P, and Rayleigh waves with special reference to the long duration observed in Mexico City. Bulletin of the Seismological Society of America, 1989, 79(5): 1361–1382

[13] Semblat J F, Duval A M, Dangla P. Seismic site effects in a deep alluvial basin: Numerical analysis by the boundary element method. Computers and Geotechnics, 2002, 29(7): 573–585

[14] Sánchez-Sesma F J, Ramos-Martínez J, Campillo M. An indirect boundary element method applied to simulate the seismic response of alluvial valleys for incident P, S and Rayleigh waves. Earthquake Engineering & Structural Dynamics, 1993, 22(4): 279–295

[15] Sánchez-Sesma F J, Luzon F. Seismic response of three-dimensional alluvial valleys for incident P, S, and Rayleigh waves. Bulletin of the Seismological Society of America, 1995, 85(1): 269–284

[16] Frankel A, Vidale J. A three-dimensional simulation of seismic waves in the Santa Clara Valley, California, from a Loma Prieta aftershock. Bulletin of the Seismological Society of America, 1992, 82(5): 2045–2074

[17] Opršal I, Zahradnik J. Three-dimensional finite difference method and hybrid modeling of earthquake ground motion. Journal of Geophysical Research. Solid Earth, 2002, 107(B8): ESE–2

[18] Wang L, Xu Y, Xia J, Luo Y. Effect of near-surface topography on high-frequency Rayleigh-wave propagation. Journal of Applied Geophysics, 2015, 116: 93–103

[19] Zhang X, Peng X, Li X, Zhou Z, Mebarki A, Dou Z, Nie W. Seismic effects of a small sedimentary basin in the eastern Tibetan plateau based on numerical simulation and ground motion records from aftershocks of the 2008 Mw7.9 Wenchuan, China earthquake. Journal of Asian Earth Sciences, 2020, 192: 104257

[20] Qiang S, Wang H, Wen R, Liu Q, Zhou Y. Investigating the effects of structural parameters on seismic aggravation of two-dimensional sedimentary valleys. Soil Dynamics and Earthquake Engineering, 2023, 171: 107964

[21] Komatitsch D, Vilotte J P. The spectral element method: An efficient tool to simulate the seismic response of 2D and 3D geological structures. Bulletin of the Seismological Society of America, 1998, 88(2): 368–392

[22] Komatitsch D, Liu Q, Tromp J, Suss P, Stidham C, Shaw J H. Simulations of ground motion in the Los Angeles basin based upon the spectral-element method. Bulletin of the Seismological Society of America, 2004, 94(1): 187–206

[23] Shen W, Yang D, Xu X, Yang S, Liu S. 3D simulation of ground motion for the 2015 Mw7.8 Gorkha earthquake, Nepal, based on the spectral element method. Natural Hazards, 2022, 112(3): 2853–2871

[24] Touhami S, Gatti F, Lopez-Caballero F, Cottereau R, de Abreu Corrêa L, Aubry L, Clouteau D. SEM3D: A 3D high-fidelity numerical earthquake simulator for broadband (0–10 Hz) seismic response prediction at a regional scale. Geosciences, 2022, 12(3): 112

[25] Ba Z, Zhao J, Zhu Z, Zhou G. 3D physics-based ground motion simulation and topography effects of the 05 September 2022 Mw6.6 Luding earthquake, China. Soil Dynamics and Earthquake Engineering, 2023, 172: 108048

[26] Noureddine M A, de Martin F, Elmeouche R, Ababsa F, Beaufils M, Sammuneh M A. Source to site seismic ground motion prediction on an earth dam. Japanese Geotechnical Society Special Publication, 2024, 10(31): 1188–1193

[27] Castro-Cruz D, Gatti F, Lopez-Caballero F, Hollender F, El Haber E, Causse M. Blind broad-band (0–10 Hz) numerical prediction of the 3-D near field seismic response of an Mw6.0 extended fault scenario: application to the nuclear site of Cadarache. Geophysical Journal International, 2023, 232(1): 581–600

[28] De Martin F, Chaljub E, Thierry P, Sochala P, Dupros F, Maufroy E, Hollender F. Influential parameters on 3D synthetic ground motions in a sedimentary basin derived from global sensitivity analysis. Geophysical Journal International, 2021, 227(3): 1795–1817

[29] Tong P, Chen C W, Komatitsch D, Basini P, Liu Q. High-resolution seismic array imaging based on an SEM-FK hybrid method. Geophysical Journal International, 2014, 197(1): 369–395

[30] Tong P, Komatitsch D, Tseng T L, Hung S H, Chen C W, Basini P, Liu Q A. 3-D spectral-element and frequency-wave number hybrid method for high-resolution seismic array imaging. Geophysical Journal International, 2014, 41(20): 7025–7034

[31] Liang J, Wu M, Ba Z, Liu Y. A hybrid method for modeling broadband seismic wave propagation in 3D localized regions to incident P, SV, and SH waves. International Journal of Applied Mechanics, 2021, 13(10): 2150119

[32] MoncadaV SLopez-caballero F. Coupled approach for the assessment of basin effects on the seismic demand of nonlinear structures. In: Proceedings of 8th International Conference on Earthquake Geotechnical Engineering. Osaka: Japanese Geotechnical Society Special Publication, 2024

[33] SotoVLopez-Caballero F. Integrated assessment of basin effects on seismic damage: A coupled 3D SEM-FEM approach with domain reduction method. 2024 (available at the website of Research Square)

[34] Korres M, Lopez-Caballero F, Alves Fernandes V, Gatti F, Zentner I, Voldoire F, Castro-Cruz D. Enhanced seismic response prediction of critical structures via 3D regional scale physics-based earthquake simulation. Journal of Earthquake Engineering, 2023, 27(3): 546–574

[35] Miller G F, Pursey H, Bullard E C. On the partition of energy between elastic waves in a semi-infinite solid. Mathematical and Physical Sciences, 1955, 233(1192): 55–69

[36] Wong H L. Effect of surface topography on the diffraction of P, SV, and Rayleigh waves. Bulletin of the Seismological Society of America, 1982, 72(4): 1167–1183

[37] Yang Y, Yu H, Yuan Y, Lu D, Huang Q. A numerical framework for underground structures in layered ground under inclined P-SV waves using stiffness matrix and domain reduction methods. Frontiers of Structural and Civil Engineering, 2023, 17(1): 10–24

[38] YangYZhou ZWangXZhangXWangZ. 2.5-dimension soil seismic response to oblique incident waves based on exact free-field solution. Frontiers of Structural and Civil Engineering, 2024, 18(2): 216-235

[39] Zhang J, Lan J, Li H, Ba Z, Zhuang H, Xu Z. Seismic wave-field modeling of 3D irregularities in layered half-space under P-, S-, and Rayleigh waves with arbitrary incident azimuth. Soil Dynamics and Earthquake Engineering, 2024, 182: 108669

[40] Nguyen K T, Kusanovic D S, Asimaki D. Three-dimensional nonlinear soil–structure interaction for Rayleigh wave incidence in layered soils. Earthquake Engineering & Structural Dynamics, 2022, 51(11): 2752–2770

[41] Bouchon M, Aki K. Discrete wave-number representation of seismic-source wave fields. Bulletin of the Seismological Society of America, 1977, 67(2): 259–277

[42] HaskellN A. Vincit Veritas: A Portrait of the Life and Work of Norman Abraham Haskell. Washington, D.C.: American Geophysical Union, 1990

[43] Zhao X, Zhu J B, Zhao J, Cai J. Study of wave attenuation across parallel fractures using propagator matrix method. International Journal for Numerical and Analytical Methods in Geomechanics, 2012, 36(10): 1264–1279

[44] Singla V K, Gupta V K. Surface rotations due to kinematic shear dislocation point source in a multilayered elastic medium. Bulletin of the Seismological Society of America, 2019, 109(1): 433–447

[45] Kausel E. Generalized stiffness matrix method for layered soils. Soil Dynamics and Earthquake Engineering, 2018, 115: 663–672

[46] Kausel E, Roësset J M. Stiffness matrices for layered soils. Bulletin of the Seismological Society of America, 1981, 71(6): 1743–1761

[47] Cao Z, Tao X, Tao Z, Tang A. Kinematic source modeling for the synthesis of broadband ground motion using the F-K approach. Bulletin of the Seismological Society of America, 2019, 109(5): 1738–1757

[48] WolfJ P. Dynamic Soil–Structure-Interaction. Englewood Cliffs, NJ: Prentice-Hall, Inc., 1985

[49] Ba Z, Liang J, Lee V W. Wave propagation of buried spherical SH-, P1-, P2-and SV-waves in a layered poroelastic half-space. Soil Dynamics and Earthquake Engineering, 2016, 88: 237–255

[50] Ba Z, Niu J, Liu Y, Liang J. Dynamic response of a multi-scale layered saturated porous half-space due to seismic dislocation source by using a revised dynamic stiffness matrix method. Applied Mathematical Modelling, 2023, 120: 217–245

[51] Ba Z, Liu Y, Liang J, Sang Q, Wu M, Zhang Y. The dynamic stiffness matrix method for seismograms synthesis for layered transversely isotropic half-space. Applied Mathematical Modelling, 2022, 104: 205–227

[52] Ba Z, Sang Q, Wu M, Liang J. The revised direct stiffness matrix method for seismogram synthesis due to dislocations: from crustal to geotechnical scale. Geophysical Journal International, 2021, 227(1): 717–734

[53] Komatitsch D, Tromp J. Spectral-element simulations of global seismic wave propagation—I. Validation. Geophysical Journal International, 2002, 149(2): 390–412

[54] Komatitsch D, Tromp J. Spectral-element simulations of global seismic wave propagation—II. Three-dimensional models, oceans, rotation and self-gravitation. Geophysical Journal International, 2002, 150(1): 303–318

[55] Clayton R, Engquist B. Absorbing boundary conditions for acoustic and elastic wave equations. Bulletin of the Seismological Society of America, 1977, 67(6): 1529–1540

[56] 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

RIGHTS & PERMISSIONS

Higher Education Press