Abstract
In this paper, threedimensional (3D) scattering of obliquely incident Rayleigh waves by a saturated alluvial valley embedded in a layered halfspace is studied in the frequency domain by using the indirect boundary element method. The total responses are assumed as the sum of the freefield responses and scatteredfield responses. The freefield responses are calculated using the direct stiffness method and the scatterfield responses outside and inside of the valley are simulated by applying two sets of fictitious moving distributed loads and pore pressures on the interface of the valley. The amplitudes of the fictitious distributed loads and pore pressures are determined from the boundary conditions. The method is validated by comparing the results with the results published in the literature, and numerical results are obtained for a saturated valley embedded in a uniform saturated halfspace and in a single elastic layer overlying elastic halfspace. The effects of incident frequency, incident angle, porosity, drainage condition, and soil layers on the dynamic responses are discussed. It is found that 3D scattering is different from the 2D case. As the porosity decreased, the pore pressures around the valley became smaller but their oscillations became violent. The wave fields in a layered site are determined by the “dynamic interaction between valley and the layered halfspace” and the dispersion characteristics of Rayleigh waves for the given layeredhalfspace.
1. Introduction
In earthquake engineering and seismology fields, it is wellrecognized that local topographies can generate large amplifications and spatial variations in seismic ground motion. For example, during the 1985 earthquake in Mexico, structures built in alluvial basin suffered heavy damage, whereas structures built over bedrock experienced significantly less damage. Extensive theoretical and experimental studies have been conducted on this subject over the past 40 years.
In studies of alluvial valley, scattering of seismic waves can be 2D (incident wave field with incident angle perpendicular to the valley) or 3D (incident wave field at an arbitrary arrive angle to the valley). Followed by the pioneering work of Aki and Larner [1] and Trifunac [2], the 2D valley effects on the wave fields have been studied by many authors using analytical or numerical methods. Analytical solutions have been restricted to simple geometries in a uniform halfspace[35], and more general cases are solved using numerical methods such as the boundary integral equation method [6, 7], the boundary element methods [810], the finite difference method [11], and the AkiLarner method [12]. Based on the Biot’s theory of elastic waves propagating in a fluidsaturated poroelastic medium [1315], the 2D scattering by a valley in a uniform saturated halfspace has also been treated using the analytical and the boundary integral equation method [1618].
Compared to the 2D scattering, studies on the 3D responses of an alluvial valley are relatively less. Pedersen et al. [19] presented the indirect boundary element method with full space Green’s functions to obtain the 3D scattering of 2D topographies in a uniform halfspace. Barros and Luco [20] used an indirect boundary integral method based on moving Green’s functions for a layered halfspace to study the responses of a valley for obliquely incident P, SV, and SH waves. Liang et al. [21] studied the 3D responses of a circulararc alluvial valley in a uniform halfspace using the wave function expansion method. However, studies on the 3D scattering by a 2D saturated valley in a uniform or layered saturated halfspace have not been reported in the literature.
Based on the dynamicstiffness matrices and Green’s functions of distributed loads acting on inclined lines in layered halfspace, a special indirect boundary element method was first proposed by Wolf [22]. The advantage of this method is its efficiency in computation and flexibility for complex boundary configurations, and the method has been extended to the 2D poroelastic case by Liang and You [23] and Liang et al. [24] and then to the 3D seismic responses of 2D topographies in an elastic layered halfspace by Ba and Liang [25]. However, the above methods cannot be directly applied to a 3D wave scattering problem by 2D saturated alluvial valley in a fluidsaturated layered halfspace, and addressing this problem can be of great significance.
In this article, an indirect boundary element method (IBEM) is presented to study the 3D scattering of obliquely incident Rayleigh waves by a saturated alluvial valley. This was performed by deriving the 3D dynamicstiffness matrices, moving Green’s functions of distributed loads, and pore pressures acting on inclined lines in a fluidsaturated layered halfspace. The accuracy of the method is verified by comparison with related results, and the effects of incident frequency, incident angle, porosity, drainage condition, and soil layers on the dynamic responses are discussed in this paper.
2. Model
The model considered in this paper is an infinitely long saturated alluvial valley embedded in a layered halfspace with its axis along the $y$axis (Fig. 1(a)). The layered halfspace is composed of horizontal saturated layers and an underlying saturated halfspace. The seismic excitation is expressed by Rayleigh waves with an arbitrary incident angle $\psi $ with respect to the $y$ axis (Fig. 1(b)). Even though the model of the valley is 2D, the responses are 3D for obliquely incident wave fields.
Fig. 1Model of 3D scattering by twodimensional valley
a)
b)
The 3D responses of the 2D valley have a unique feature with the wave fields at two different cross sections of the valley being identical but with a shift in phase in the frequency domain. This allows only one cross section of the valley be used for calculation to obtain the whole 3D dynamic responses, while wave fields at other sections can be obtained by phase offset. In this paper, only the interface of the valley is discretized into inclined straight lines for halfspace Green’s functions are used in the proposed IBEM. The total responses are assumed to be the sum of the freefield and scatteredfield responses. Firstly, the freefield responses are calculated to determine the displacements and stresses at the interface of the valley, and then one set of fictitious distributed loads and pore pressures moving parallel to the $y$ axis are applied at the line elements (interface of the valley) to represent the scattered wave field outside of the valley (scattered wave field in the layered halfspace), and another set of distributed loads and pore pressures on the same line elements to represent the scattered wave field inside of the valley (Fig. 2). The amplitudes of the two sets of fictitious moving distributed loads and pore pressures can be determined by the boundary conditions at the interface of the valley, and finally, the dynamic responses arising from the waves in the freefield and from the fictitious moving distributed loads and pore pressures are summed to obtain the total responses.
Fig. 2Discretization of the interface of the alluvial valley
2.1. Biot’s theory
The equations of motion of a poroelastic medium [1315] can be expressed as:
$\alpha M\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\left(\mathrm{d}\mathrm{i}\mathrm{v}\mathbf{U}\right)+M\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\left(\mathrm{d}\mathrm{i}\mathrm{v}\mathbf{w}\right)={\rho}_{f}\ddot{\mathbf{U}}+m\ddot{\mathbf{w}}+b\dot{\mathbf{w}},$
where, $\mathbf{U}$ and $\mathbf{w}$ are displacement vectors of the solid frame and of the pore fluid with respect to the solid frame, respectively; $\mathbf{w}=n(\mathbf{W}\mathbf{U})$; $\mathbf{W}$ is the displacements vector of the pore fluid; $\rho =(1n){\rho}_{s}+n{\rho}_{f}$ represents the total density; ${\rho}_{s}$ is solid grain’s density; ${\rho}_{f}$ is pore fluid’s density; $n$ is porosity; $b$ is a parameter that accounts for internal friction due to relative motion between the solid frame and the pore fluid; $m={\rho}_{22}/{n}^{2}=(n{\rho}_{f}+{\rho}_{a})/{n}^{2}$ represents a densitylike parameter; ${\rho}_{a}$ is the couple density between the solid grain and the pore fluid; $G$ and $\lambda $ are Lame’s constants of the solid frame, ${\lambda}_{c}=\lambda +{\alpha}^{2}M$; $\alpha $ and $M$ are Biot’s parameters, which account for compressibility of the two phase material.
To solve the Biot’s governing equations in Cartesian coordinate system, two kinds of Fourier transformation are involved: one with respect to time and frequency, and the other with respect to the two horizontal coordinates and wave numbers. In this paper, the Fourier transformations for the two horizontal coordinates and time are defined as follows:
$f\left(x,y,z,t\right)\text{=}{\int}_{\infty}^{\infty}{\int}_{\infty}^{\infty}{\int}_{\infty}^{\infty}F\left({k}_{x},{k}_{y},z,\omega \right){e}^{i\omega ti{k}_{x}xi{k}_{y}y}d{k}_{x}d{k}_{y}d\omega ,$
where $x$, $y$ and $t$ denote the two horizontal coordinates and time, while ${k}_{x}$, ${k}_{y}$ and $\omega $ represent the two horizontal wave numbers and frequency. Performing the Fourier transform on Eq. (1) and solving the resulting ordinary differential equations, the displacement amplitudes of the solid frame and of the pore fluid with respect to the solid frame in the frequency and wave number domain can be written as:
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{k}_{Sz}{k}_{x}\beta \left({A}_{SV}{e}^{i{k}_{Sz}z}{B}_{SV}{e}^{i{k}_{Sz}z}\right){k}_{y}\beta \left({A}_{SH}{e}^{i{k}_{Sz}z}+{B}_{SH}{e}^{i{k}_{Sz}z}\right),$
${U}_{y}={k}_{y}\left[{a}_{1}\left({A}_{P1}{e}^{i{k}_{P1z}z}+{B}_{P1}{e}^{i{k}_{P1z}z}\right)+{a}_{2}\left({A}_{P2}{e}^{i{k}_{P2z}z}+{B}_{P2}{e}^{i{k}_{P2z}z}\right)\right]$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{k}_{Sz}{k}_{y}\beta \left({A}_{SV}{e}^{i{k}_{Sz}z}{B}_{SV}{e}^{i{k}_{Sz}z}\right)+{k}_{x}\beta \left({A}_{SH}{e}^{i{k}_{Sz}z}+{A}_{SH}{e}^{i{k}_{Sz}z}\right),$
${U}_{z}={k}_{P1z}{a}_{1}\left({A}_{P1}{e}^{i{k}_{P1z}z}{B}_{P1}{e}^{i{k}_{P1z}z}\right)+{k}_{P2z}{a}_{2}\left({A}_{P2}{e}^{i{k}_{P2z}z}{B}_{P2}{e}^{i{k}_{P2z}z}\right)$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\frac{\left({A}_{SV}{e}^{i{k}_{Sz}z}+{B}_{SV}{e}^{i{k}_{Sz}z}\right)}{\beta},$
${w}_{z}={k}_{P1z}{a}_{1}{D}_{1}\left({A}_{P1}{e}^{i{k}_{P1z}z}{B}_{P1}{e}^{i{k}_{P1z}z}\right)+{k}_{P2z}{a}_{2}{D}_{2}\left({A}_{P2}{e}^{i{k}_{P2z}z}{B}_{P2}{e}^{i{k}_{P2z}z}\right)$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\frac{{D}_{3}\left({A}_{SV}{e}^{i{k}_{Sz}z}+{B}_{SV}{e}^{i{k}_{Sz}z}\right)}{\beta}.$
where:
$\beta =\frac{1}{\sqrt{{k}_{x}^{2}+{k}_{y}^{2}}},{D}_{i}=\frac{\left({\lambda}_{c}+2\mu \rho {V}_{Pi}^{2}\right)}{\left(\rho {V}_{Pi}^{2}\alpha M\right)},i=1,\mathrm{}2,\mathrm{}\mathrm{}\mathrm{}{D}_{3}=\frac{{\rho}_{f}{\omega}^{2}}{\left(ib\omega m{\omega}^{2}\right)},$
and ${A}_{P1}$, ${B}_{P1}$, ${A}_{P2}$, ${B}_{P2}$, ${A}_{SV}$, ${B}_{SV}$, ${A}_{SH}$ and ${B}_{SH}$ are the amplitudes of the upgoing and downgoing $P1$, $P2$, $SV$ and $SH$ waves, respectively, and ${k}_{P1}$, ${k}_{P2}$ and ${k}_{S}$ are the wave numbers associated with the $P1$, $P2$ and $S$ waves, which are determined by Eq. (4):
where:
${\beta}_{2}=\frac{\left[\left(m{\omega}^{2}ib\omega \right)\rho {\omega}^{2}{\rho}_{f}^{2}{\omega}^{4}\right]}{\left[M\left(\lambda +2\mu \right)\right]}.$
2.2. 3D stiffness matrix of saturated, layered halfspace and solutions of the free field responses
The constitutive equations for a poroelastic medium have the following form [1315]:
where ${\delta}_{ij}$ is the Dirac delta function; $e={U}_{i,i}$ represents the dilation of the solid skeleton; $\theta ={w}_{i,i}$ represents the volume of fluid injected into a unit volume bulk material; ${\sigma}_{i,j}$ is the stresses in the bulk porous medium; and ${p}_{f}$ is the pore pressure.
The displacement amplitudes at the top and the bottom of a saturated layer can be expressed as Eq. (7a) developed based on Eqs. (3), and the external load amplitudes can be expressed as Eq. (7b) based on Eqs. (3) and Eqs. (6), by introducing ${Q}_{x1}={\tau}_{zx1}$, ${Q}_{y1}={\tau}_{zy1}$, ${Q}_{z1}=i{\sigma}_{z1}$, ${Q}_{f1}=i{p}_{f1}$, ${Q}_{x2}={\tau}_{zx2}$, ${Q}_{y2}={\tau}_{zy2}$, ${Q}_{z2}=i{\sigma}_{z2}$ and ${Q}_{f2}$=$i{p}_{f2}$:
where ${U}_{x1}$, ${U}_{y1}$, ${U}_{z1}$, ${U}_{x2}$, ${U}_{y2}$ and ${U}_{z2}$ are the solid frame’s displacement amplitudes; ${w}_{z1}$ and ${w}_{z2}$ are the displacement amplitudes of pore fluid with respect to the solid frame; ${\tau}_{zx1}$, ${\tau}_{zy1}$, ${\sigma}_{z1}$, ${\tau}_{zx2}$, ${\tau}_{zy2}$ and ${\sigma}_{z2}$ are the total stresses amplitudes; and ${p}_{f1}$ and ${p}_{f2}$ are the pore pressure amplitudes. To achieve symmetry, all vertical components in Eq. (7) are multiplied with $i$. Eliminating the amplitudes ${A}_{P1}$, ${B}_{P1}$, ${A}_{P2}$, ${B}_{P2}$, ${A}_{SV}$, ${B}_{SV}$, ${A}_{SH}$ and ${B}_{SH}$ in Eq. (7), leads to the 3D dynamic stiffness matrix of the saturated soil layer ${\mathbf{S}}_{P1P2SVSH}^{L}$. Applying loads at the free surface of the saturated halfspace, only downgoing waves with amplitudes ${B}_{P1}$, ${B}_{P2}$, ${B}_{SV}$, and ${B}_{SH}$ are developed, as the radiation condition states that no energy propagates from infinity, it excludes the upgoing waves with ${A}_{P1}\text{,}$${A}_{P2}\text{,}$${A}_{SV}$ and ${A}_{SH}$. Therefore, by setting ${A}_{p1}={A}_{p2}={A}_{SV}={A}_{SH}=$0 in Eq. (8) and eliminating ${B}_{P1}$, ${B}_{P2}$, ${B}_{SV}$ and ${B}_{SH}$, the saturated halfspace matrix ${\mathbf{S}}_{P1P2SVSH}^{R}$ can be obtained. Considering the continuity of displacements and equilibrium of stresses and pore pressure at the layer interfaces, the global 3D stiffness matrix of a multilayered site ${\mathbf{S}}_{P1P2SVSH}$ are obtained by assembling the layer matrices ${\mathbf{S}}_{P1P2SVSH}^{L}$ and the halfspace stiffness matrix ${\mathbf{S}}_{P1P2SVSH}^{R}$. Once ${\mathbf{S}}_{P1P2SVSH}$ is obtained, the discretized dynamicequilibrium equation of the layered saturated site can be expressed as:
where, $\mathbf{U}$ and $\mathbf{Q}$ are the vector of the displacement amplitudes and vector of external load amplitudes at the interfaces of the saturated layer.
To obtain freefield responses of the incident Rayleigh waves in a layered halfspace, the relationship between the apparent velocity $c$ and the frequency $\omega $ must be established. Using the 3D dynamic stiffness matrix of the layered site, the relationship between $c$ and $\omega $ can be established by setting $\left{\mathbf{S}}_{P1P2SVSH}\right=0$. Newton iteration method can be used to obtain the $c\left(\omega \right)$ for the given $\omega $. Substituting these values of $c\left(\omega \right)$ in the global dynamic matrix leads to nontrivial solutions for the dynamicstiffness matrix ${\mathbf{S}}_{P1P2SVSH}$, which corresponds to the socalled Rayleigh wave modes. Thus, the ratios of displacement amplitudes at the layers’ interfaces are obtained. Once the displacements at the interfaces of the layers are obtained, the amplitudes of the upgoing waves and downgoing waves in the layers can be determined using Eq. (8a). Then the freefield responses are obtained, which consist of the displacements of the solid frame in the $x$, $y$ and $z$ directions and the displacement of the pore fluid with respect to the solid phase ${\left\{{\mathbf{U}}_{xf}\left(S\right),{\mathbf{U}}_{yf}\left(S\right),{\mathbf{U}}_{zf}\left(S\right),\mathrm{}{\mathbf{w}}_{zf}\left(\mathrm{S}\right)\right\}}^{T}$, and the stresses in the $x$, $y$ and $z$ directions and pore pressure ${\left\{{\mathbf{t}}_{xf}\left(S\right),{\mathbf{t}}_{yf}\left(S\right),{\mathbf{t}}_{zf}\left(S\right),{\mathbf{t}}_{pf}\left(S\right)\right\}}^{T}$ (Fig. 2). The subscript ‘$f$’ indicates the variables that represent terms in the freefield responses.
2.3. Moving Green’s functions and solutions of the scattered wave field
The Green’s functions used in this paper denote the dynamic responses at any point in the freefield system (without the valley) when moving distributed loads in the $x$, $y$ and $z$ directions and the pore pressures acting on an inclined line embedded in the saturated layered halfspace (Fig. 3). These versions of green’s functions for a layered halfspace originate from Wolf [22], which were extended to the poroelastic case by Liang et al. [23] and to the 2.5 D form by Ba and Liang [25]. Taking the distributed load in the $x$ direction moving along the $y$ axis as an example, the load amplitude in the time and space domain, arising from the value of density ${p}_{x0}$ can be expressed as:
where, $\theta $ (0° $<\theta <$ 180°) is the angle of the line measured from the horizontal plane and $c$ is the moving velocity. $c={c}_{R}^{1}/\mathrm{c}\mathrm{o}\mathrm{s}\psi $ can be obtained from Fig. 1(b) for Rayleigh waves with the arriving angle $\psi $ with respect to the $y$ axis. ${c}_{R}^{1}$ is the Rayleigh wave’s velocity of the top soil layer. Performing Fourier transforms with respect to the two horizontal coordinates and also with respect to time on Eq. (9), the load amplitude in the frequency and wave number domain can be expressed as:
Fig. 3Distributed load acting on an inclined line in soil layers
As only part of the layer is loaded, two additional interfaces are introduced at the top and the bottom boundary of the load. First, the layer on which the distributed load acts is fixed at the two interfaces, and the corresponding reaction forces ${R}_{x1}$, ${R}_{y1}$, ${R}_{z1}$, ${R}_{f1}$ at the top interface and ${R}_{x2}$, ${R}_{y2}$, ${R}_{z2}$, ${R}_{f2}$ at the bottom interface) are calculated to achieve this condition. This way, the analysis is restricted to the loaded layer. The directions of these forces are then reversed, and are applied as external loads to the total system to determine the global response, and the dynamic responses of the total system induced by the external forces can be implemented using the direct stiffness method. The total response is finally obtained by adding the results of analysis of the fixed layer of the reaction forces (Fig. 3).
It should be noted that the above calculations are performed in the frequency and wave number domain, and the Green’s functions in the time and space domain can be obtained using inverse Fourier transformation given by Eq. (11):
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}={\int}_{\mathrm{\infty}}^{\mathrm{\infty}}F\left({k}_{x},{k\mathrm{\text{'}}}_{y},z,\omega \right)\mathrm{e}\mathrm{x}\mathrm{p}\left(i{k}_{x}x\right)d{k}_{x},$
where ${k}_{y}^{\mathrm{\text{'}}}=\omega /c$, $F\left({k}_{x},{k}_{y},z,\omega \right)$ is the dynamic response in the frequency and wave number domain, and $f(x,y,z,w)$ is the response in the frequency and space domain$.$
Using the Green’s functions described above, the scattered wave fields outside and inside of the valley are simulated by two sets of moving Green’s functions for the layered half space and for the valley, respectively. Therefore, when the free surface is drained (permeable), the displacements and stresses at the interface of the valley $S$ arising from the scattered wave fields outside and inside of the valley can be expressed as:
${\left\{{\mathbf{t}}_{xg}^{L}\left(S\right),{\mathbf{t}}_{yg}^{L}\left(S\right),{\mathbf{t}}_{zg}^{L}\left(S\right)\uff0c{\mathbf{t}}_{pg}^{L}\left(S\right)\right\}}^{T}=\left[{g}_{t}^{L}\left(S\right)\right]{\left\{\mathbf{p}{}_{x1}{}^{\mathrm{}}\mathrm{},\mathbf{p}{}_{y1}{}^{\mathrm{}}\mathrm{},\mathbf{p}{}_{z1}{}^{\mathrm{}}\mathrm{},{\mathbf{p}}_{f1}\right\}}^{T},$
${\left\{{\mathbf{U}}_{xg}^{V}\left(S\right),{\mathbf{U}}_{yg}^{V}\left(S\right),{\mathbf{U}}_{zg}^{V}\left(S\right)\right\}}^{T}=\left[{g}_{u}^{V}\left(S\right)\right]{\left\{\mathbf{p}{}_{x2}{}^{\mathrm{}}\mathrm{},\mathbf{p}{}_{y2}{}^{\mathrm{}}\mathrm{},\mathbf{p}{}_{z2}{}^{\mathrm{}}\mathrm{},{\mathbf{p}}_{f2}\right\}}^{T},$
${\left\{{\mathbf{t}}_{xg}^{V}\left(S\right),{\mathbf{t}}_{yg}^{V}\left(S\right),{\mathbf{t}}_{zg}^{V}\left(S\right),{\mathbf{t}}_{pg}^{V}\left(S\right)\right\}}^{T}=\left[{g}_{t}^{V}\left(S\right)\right]{\left\{\mathbf{p}{}_{x2}{}^{\mathrm{}}\mathrm{},\mathbf{p}{}_{y2}{}^{\mathrm{}}\mathrm{},\mathbf{p}{}_{z2}{}^{\mathrm{}}\mathrm{},{\mathbf{p}}_{f2}\right\}}^{T},$
where $\left[{g}_{u}^{L}\left(S\right)\right]$ and $\left[{g}_{t}^{L}\left(S\right)\right]$ are the displacements and stresses Green’s functions for the layered halfspace with a drained free surface, and $\left[{g}_{u}^{V}\left(S\right)\right]$ and $\left[{g}_{t}^{V}\left(S\right)\right]$ are the displacements and stresses Green’s functions for the valley. In Eq. (12), the subscript ‘$g$’ indicates results attributable to the moving distributed loads and pore pressures, while the superscripts ‘$L$’ and ‘$V$’ denote parameters associated with the layered halfspace and with the valley, respectively. ${\left\{\mathbf{p}{}_{x1}{}^{\mathrm{}}\mathrm{},\mathbf{p}{}_{y1}{}^{\mathrm{}}\mathrm{},\mathbf{p}{}_{z1}{}^{\mathrm{}}\mathrm{},{\mathbf{p}}_{f1}\right\}}^{T}$ and ${\left\{\mathbf{p}{}_{x2}{}^{\mathrm{}}\mathrm{},\mathbf{p}{}_{y2}{}^{\mathrm{}}\mathrm{},\mathbf{p}{}_{z2}{}^{\mathrm{}}\mathrm{},{\mathbf{p}}_{f2}\right\}}^{T}$ are the load and pore pressure amplitude vectors acting on the interface $S$ to calculate the green’s functions for layered halfspace and for the valley, respectively.
Similarly, when the free surface is undrained (impermeable), the displacements and stresses at the interface of the valley $S$ arising from the scattered wave fields outside and inside of the valley can be expressed as:
${\left\{{\mathbf{t}}_{xg}^{L}\left(S\right),{\mathbf{t}}_{yg}^{L}\left(S\right),{\mathbf{t}}_{zg}^{L}\left(S\right),{\mathbf{w}}_{zg}^{L}\left(S\right)\right\}}^{T}=\left[{g}_{ut}^{L}\left(S\right)\right]{\left\{\mathbf{p}{}_{x1}{}^{\mathrm{}}\mathrm{},\mathbf{p}{}_{y1}{}^{\mathrm{}}\mathrm{},\mathbf{p}{}_{z1}{}^{\mathrm{}}\mathrm{},{\mathbf{p}}_{f1}\right\}}^{T},$
${\left\{{\mathbf{U}}_{xg}^{V}\left(S\right),{\mathbf{U}}_{yg}^{V}\left(S\right),{\mathbf{U}}_{zg}^{V}\left(S\right)\right\}}^{T}=\left[{g}_{uu}^{V}\left(S\right)\right]{\left\{\mathbf{p}{}_{x2}{}^{\mathrm{}}\mathrm{},\mathbf{p}{}_{y2}{}^{\mathrm{}}\mathrm{},\mathbf{p}{}_{z2}{}^{\mathrm{}}\mathrm{},{\mathbf{p}}_{f2}\right\}}^{T},$
${\left\{{\mathbf{t}}_{xg}^{V}\left(S\right),{\mathbf{t}}_{yg}^{V}\left(S\right),{\mathbf{t}}_{zg}^{V}\left(S\right),{\mathbf{w}}_{zg}^{V}\left(S\right)\right\}}^{T}=\left[{g}_{ut}^{V}\left(S\right)\right]{\left\{\mathbf{p}{}_{x2}{}^{\mathrm{}}\mathrm{},\mathbf{p}{}_{y2}{}^{\mathrm{}}\mathrm{},\mathbf{p}{}_{z2}{}^{\mathrm{}}\mathrm{},{\mathbf{p}}_{f2}\right\}}^{T},$
where, $\left[{g}_{uu}^{L}\left(S\right)\right]$, $\left[{g}_{ut}^{L}\left(S\right)\right]$, $\left[{g}_{uu}^{V}\left(S\right)\right]$ and $\left[{g}_{ut}^{V}\left(S\right)\right]$ are the displacements and stresses Green’s functions for the layered halfspace and for the valley with an undrained free surface, respectively.
2.4. Boundary conditions
If the free surface of the layered halfspace and the interface of the valley have a drained boundary condition, the boundary conditions at the valley interface $S$ can be expressed as:
${{\int}_{s}\left[W\left(s\right)\right]}^{T}\left(\left[\begin{array}{c}{\mathbf{U}}_{xg}^{L}\left(S\right)\\ {\mathbf{U}}_{yg}^{L}\left(S\right)\\ {\mathbf{U}}_{zg}^{L}\left(S\right)\\ {\mathbf{w}}_{g}^{L}\left(S\right)\end{array}\right]+\left[\begin{array}{c}{\mathbf{U}}_{xf}\left(S\right)\\ {\mathbf{U}}_{yf}\left(S\right)\\ {\mathbf{U}}_{zf}\left(S\right)\\ {\mathbf{w}}_{f}\left(S\right)\end{array}\right]\right)ds={{\int}_{s}\left[W\left(s\right)\right]}^{T}\left[\begin{array}{c}{\mathbf{U}}_{xg}^{V}\left(S\right)\\ {\mathbf{U}}_{yg}^{V}\left(S\right)\\ {\mathbf{U}}_{zg}^{V}\left(S\right)\\ {\mathbf{w}}_{g}^{V}\left(S\right)\end{array}\right]ds,$
where, $\left[W\left(s\right)\right]$ is the weighting function. If $\left[W\left(s\right)\right]$ takes the value 1 in a single element and zero in all the others, the integral can be evaluated over each element separately. Substituting Eq. (12) into Eq. (14) results in:
$\left[{V}_{p}^{L}\right]{\left\{{\mathbf{p}}_{x1},{\mathbf{p}}_{y1},{\mathbf{p}}_{z1},{\mathbf{p}}_{f1}\right\}}^{T}+\left[{V}_{f}\right]=\left[{V}_{p}^{V}\right]{\left\{{\mathbf{p}}_{x2},{\mathbf{p}}_{y2},{\mathbf{p}}_{z2},{\mathbf{p}}_{f2}\right\}}^{T},$
where:
$\left[{T}_{f}\right]={\int}_{s}{\left[W\left(s\right)\right]}^{T}\left[{t}_{f}\left(s\right)\right]ds,\left[{V}_{p}^{L}\right]={\int}_{s}{\left[W\left(s\right)\right]}^{T}\left[{g}_{u}^{L}\left(s\right)\right]ds,$
$\left[{V}_{p}^{V}\right]={\int}_{s}{\left[W\left(s\right)\right]}^{T}\left[{g}_{u}^{V}\left(s\right)\right]ds,\left[{V}_{f}\right]={\int}_{s}{\left[W\left(s\right)\right]}^{T}\left[{v}_{f}\left(s\right)\right]ds.$
Combining the solutions of Eqs. (15) and (12), the displacements and stresses outside of the valley will be:
$\left\{\begin{array}{c}{t}_{x}\left(x,y,z\right)\\ {t}_{y}\left(x,y,z\right)\\ {t}_{z}\left(x,y,z\right)\\ {t}_{f}\left(x,y,z\right)\end{array}\right\}=\left\{\begin{array}{c}{t}_{x}\left(x,y,z\right)\\ {t}_{y}\left(x,y,z\right)\\ {t}_{z}\left(x,y,z\right)\\ {t}_{pf}\left(x,y,z\right)\end{array}\right\}+\left[{g}_{t}^{L}\left(x,y,z\right)\right]{\left[\begin{array}{l}{\mathbf{p}}_{x1}\\ {\mathbf{p}}_{y1}\\ {\mathbf{p}}_{z1}\\ {\mathbf{p}}_{f1}\end{array}\right]}^{T},$
and the displacements and the stresses inside of the valley will be:
If the free surface of the layered halfspace and the interface of the valley are both impermeable (undrained boundary), the boundary conditions are: continuity of the displacements of soil frame, equilibrium of the stresses, and zero flow through the interface of the valley. Similar to the calculation process described above for the drained boundary case, the responses outside and inside of the valley can be obtained for the undrained boundary case.
Fig. 4Surface displacement amplitudes compared with those of Dravinski and Mossessian [7]
3. Verification of the method
The present method can be reduced to a 2D problem in an elastic medium when the incident angle $\psi =\mathrm{}$90° and the saturated parameters ${\rho}_{f}$, $M$, $m$, $\alpha $ and $b$ are set as small values (i.e., 10^{3}). The 2D results of Rayleigh waves for a semicircular valley embedded in a uniform elastic halfspace provided by Dravinski and Mossessian [7] are used for comparison. The parameters used in the analysis are as follows: Poisson’s ratios of both the valley and the halfspace are $\nu =$1/3, damping ratio is $\zeta =$0.005, the shear wave velocity ratio of the valley to the halfspace is ${c}_{s}^{V}/{c}_{s}^{H}=$0.5, the mass density ratio is ${\rho}_{s}^{V}/{\rho}_{s}^{H}=$2/3 (note that the superscript ‘$V$’ represents the valley and ‘$H$’ represents the halfspace). The incident frequency $\eta =$0.5 and 1.0 are defined as $\eta =\omega a/\pi {c}_{s}^{H}$, where, $a$ is the radius of the valley. Fig. 4 shows that the results from the analysis are in good agreement with those provided by Dravinski and Mossessian [7].
4. Numerical results and discussion
4.1. Responses of a valley in a uniform halfspace
A semicircular saturated valley embedded in a uniform fluidsaturated halfspace is first analyzed for both drained and undrained boundary conditions. The material parameters of the halfspace and of the valley are presented in Table 1 for the case of $n=$0.3. The internal friction due to relative motion between the solid frame and the pore fluid is not considered in this analysis (i.e. $b=$0). The dimensionless incident frequency is defined as $\eta =\omega a/\left(\pi \sqrt{G/\rho}\right)$, where $G$ and $\rho $ are shear modulus and total density of the saturated halfspace, respectively, and $a$ is the radius of the valley. Fig. 5 presents the amplitudes of the surface displacement for dimensionless incident frequencies $\eta =$ 1.0 and 2.0, incident angles $\psi =$30°, 60° and 90°. The solid frame displacement amplitudes in the $x\text{,}$$y$ and $z$ directions are defined as $\left{\stackrel{}{U}}_{x}\right=\left{U}_{x}/({U}_{x0}+{U}_{y0})\right\text{,}$$\left{\stackrel{}{U}}_{y}\right=\left{U}_{y}/({U}_{x0}+{U}_{y0})\right$ and $\left{\stackrel{}{U}}_{z}\right=\left{U}_{z}/({U}_{x0}+{U}_{y0})\right$, where ${U}_{x0}$ and ${U}_{y0}$ are the two freefield horizontal displacements of the obliquely incident Rayleigh waves.
Table 1Parameters for the saturated valley and the saturated halfspace
$n$  $G$ (Pa)  $M$ (Pa)  $m$ (Pa)  ${\rho}_{s}$ (kg/m^{3})  ${\rho}_{f}$ (kg/m^{3})  $\nu $  $\alpha $  
Halfspace  0.1  156.33E8  182.16E8  55000  2650  1000  0.25  0.2762 
0.3  37E8  60.72E8  7222  2650  1000  0.25  0.8287  
Alluvial valley  0.1  39.08E8  143.71E8  55000  2650  1000  0.25  0.2762 
0.3  9.25E8  47.92E8  7222  2605  1000  0.25  0.8287 
Fig. 5 indicates that the presence of the saturated valley can cause very large amplification of the surface ground motion locally, which is strongly dependent upon the incident angles and the drainage condition. The repeated interference of the transmitted waves in the saturated valley leads to a large amplification on the surface of the valley. The relatively shorter wavelengths in the valley made the amplitude oscillations on the surface more violent. The displacement amplitudes in the left side (i.e., the incident side) of the valley are more complex than those on the right side because of the filtering effect of the valley.
The results in Fig. 5 also show that the 3D responses (i.e., for $\psi \ne $ 90°) are different from the 2D case (i.e., for $\psi =\mathrm{}$90°). These results indicate that the obliquely incident wave fields cannot be simply decomposed into inplane wave fields and antiplane wave fields, and then be calculated as the total inplane and antiplane responses. As the incident angle increased, the displacement amplitudes in the $x$ direction increased gradually, while the displacement amplitudes in the $y$ direction decreased gradually. For the 2D case ($\psi =$90°), only the inplane displacements (displacements in the $x$ and $z$ directions) are present. The drainage condition also has a significant effect on the surface displacement amplitudes, and the amplitudes are larger and more complex for the undrained boundary condition than those for the drained case. The reason for this phenomenon can be attributed to a stronger ‘encapsulation effect’ in the case of an undrained saturated valley.
Fig. 6 shows the contours of the pore pressure amplitudes in the range of $2a\le x\le 2a$ and $0\le z\le 2a$ for $\eta =$0.5 and 2.0. The incident angles are $\psi =$30° and 60° and the drainage condition is undrained. The calculation model, material parameters, and also the definition of the dimensionless frequency used are same as for the results presented in Fig. 5. The dimensionless pore pressure amplitudes are defined as $\left{\stackrel{}{p}}_{f}\right=\left{p}_{f}/\left(G\right({U}_{x0}+{U}_{y0}\left)\right)\right$.
It can be seen that larger amplitudes and more violent oscillations of pore pressures are found inside the valley compared to those in the outside saturated halfspace. The pore pressure amplitudes at the right side of the valley are much larger, in contrary to the distribution of the displacement amplitudes. As the incident frequency increased, the amplitudes of the pore pressure increased significantly and the oscillations became more violent. Fig. 6 also shows that the amplitudes of the pore pressure became larger with increasing incident angle, and increasing incident angle increased the differences in pore pressure amplitudes between the left and right side of the valley.
Fig. 7 shows the contours of the amplitudes of the pore pressure in the range of $2a\le x\le 2a$ and $0\le z\le 2a$ for $\eta =$0.5, with an undrained boundary condition and $\psi =$ 30°^{}and 60°. The material parameters of the saturated halfspace and of the saturated valley are presented in Table 1 for the case of $n=$0.1 and the internal friction is not considered ($b=$0). Comparisons between the results in Fig. 7 ($n=$0.3) and results in Fig. 6 ($n=$0.1) show that porosity has significant effect on wave scattering around the valley. As the porosity decreased, the pore pressure amplitudes decreased significantly, while the oscillations of the pore pressures became more violent. This indicates that less energy is taken up by the pore fluid than by the solid frame. The results in Fig. 7 also show that the pore pressure amplitudes become more concentrated inside the valley for smaller porosity (i.e. $n=$0.1).
Fig. 5Surface displacements amplitudes for different obliquely incident angels (n= 0.3)
Fig. 6Contours of pore pressure amplitudes around the valley (n= 0.3)
Fig. 7Contours of pore pressure amplitudes around the valley (n= 0.1)
4.2. Responses of a valley in a layered halfspace
A semicircular valley embedded in a single elastic layer overlying an elastic halfspace is analyzed considering both drained and undrained boundary conditions. The material parameters of the saturated valley are shown in Table 1 for the case of $n=$0.3 and material damping ratio of 5 % is considered in this analysis. The material parameters of the elastic layer are as follows: mass density ${\rho}^{L}=$1855 KN/m^{3}, shear modulus ${G}^{L}=$3.7 MPa, Poisson’s ratio ${v}^{L}=$0.25 and damping ratio ${\zeta}^{L}=$5 %. The material parameters of the underlying elastic halfspace are as follows: ${\rho}^{R}=$1855 KN/m^{3}, shear modulus ratio of the halfspace to the layer is $\sqrt{{G}^{R}/{G}^{L}}=$5.0, Poisson’s ratio ${v}^{R}=$0.25 and damping ratio ${\zeta}^{R}=$2 %. The thickness of the layer is $H/a=$2.0 and the incident angle is $\psi =$45°. The dimensionless incident frequency is defined as $\eta =\omega a/\left(\pi \sqrt{{G}^{L}/{\rho}^{L}}\right)$. Fig. 8 shows the amplitudes of the surface displacement for the first three modes of the Rayleigh waves for $\eta =$1.0 and 2.0. The displacement amplitudes in the $x$, $y$ and $z$ directions are defined as $\left{\stackrel{}{U}}_{x}\right=\left{\stackrel{}{U}}_{x}/({U}_{x0}+{U}_{y0})\right\text{,}$$\left{\stackrel{}{U}}_{y}\right=\left{\stackrel{}{U}}_{y}/({U}_{x0}+{U}_{y0})\right$ and $\left{\stackrel{}{U}}_{z}\right=\left{\stackrel{}{U}}_{z}/({U}_{x0}+{U}_{y0})\right$, where ${U}_{x0}$ and ${U}_{y0}$ are the two freefield horizontal displacements of the obliquely incident Rayleigh waves in the layered halfspace.
Results presented in Fig. 8 indicate that the surface displacement amplitudes can vary with different Rayleigh modes, depending on the dimensionless frequency and drainage conditions. Largest vertical displacement amplitudes are found in the third Rayleigh mode with maximum vertical displacement amplitudes at 2.2, 2.4 and 9.4, for Rayleigh modes 1, 2 and 3, respectively, for $\eta =$ 1.0. The largest vertical displacement amplitudes occurred at the second mode for $\eta =$2.0 with the maximum vertical displacement amplitudes at 2.2, 4.1 and 1.2 for the Rayleigh mode 1, 2 and 3, respectively, for $\eta =$2.0).
Fig. 8Surface displacement amplitudes for the first three modes (GR/GL= 5, H/a= 2.0)
Fig. 9 shows pore pressure amplitude contours inside of the valley for the first three Rayleigh modes, with dimensionless frequency $\eta =$ 1.0 and 2.0, and incident angle $\psi =$45°. Both the drained and undrained boundary conditions are considered in the analysis. The calculation model and material parameters are same as those in Fig. 8. The pore pressure amplitudes in Fig. 9 are defined as $\left{\stackrel{}{p}}_{f}\right=\left{p}_{f}/\left({G}^{L}\right({U}_{x0}+{U}_{y0}\left)\right)\right$.
Fig. 9 shows that both the values and the spatial distributions of the pore pressure amplitudes are different for different Rayleigh modes, depending on the incident frequency, incident angle, and the drainage condition. The pore pressure amplitude is largest at the third Rayleigh mode for $\eta =$1.0, while it is largest at the second Rayleigh mode for $\eta =$2.0. The oscillations of the pore pressure amplitudes are most violent for the first Rayleigh mode, less violent for the second Rayleigh mode, and the least violent for the third Rayleigh mode, which may be attributed to the faster apparent velocity (corresponding to longer wavelength) for the higher Rayleigh modes (Table 2). With increasing incident frequency, the pore pressure amplitudes become larger and the spatial distributions become more complex. It can also be seen for Fig. 9 that the drainage boundary condition has significant effects on the values and spatial distributions of the pore pressure amplitudes and the oscillations are more violent for the undrained boundary condition.
Fig. 10 shows the amplitudes of the surface displacement and Fig. 11 shows the contours of pore pressure amplitudes inside of the valley for the first two Rayleigh modes, illustrating the effects of the stiffness of the soil layer. In this results, the same calculation model and the same materials as in Fig. 9 are used, except the shear modulus ratio of the halfspace to the layer is ${G}^{R}/{G}^{L}=$2.0.
By comparing the results in Fig. 8 and Fig. 9 with Fig. 10 and Fig. 11 for $\eta =$1.0, it can be seen that the values and the spatial distributions of displacement (pore pressure) amplitudes are similar for the first Rayleigh mode, but are significantly different for the second mode. This is because for the scattered Rayleigh waves in a layered site, the resulting responses depend on two aspects: (1) the dynamic interaction between the valley and the layered site, and (2) the dispersion characteristics of the Rayleigh waves. The large difference in the apparent velocity of the second mode (Table 2) leads to the significant difference in the dynamic responses in Figs. 89 and in Figs. 1011, although the dynamic interaction keep invariant for the unchanged thickness of the layer (The dynamic interaction is determined by the thickness of the layer for invariant shape of valley).
Fig. 9Contours of pore pressure amplitudes for the first three modes (GR/GL= 5, H/a= 2.0)
Fig. 10Surface displacement amplitudes for the first two modes (GR/GL= 52, H/a= 2.0)
Fig. 11Contours of pore pressure amplitudes for the first two modes (GR/GL= 2, H/a= 2.0)
Fig. 12 shows the amplitudes of the surface displacement and Fig. 13 shows the contours of pore pressure amplitudes inside of the valley for the first three Rayleigh modes, illustrating the effects of the soil layer thickness. The same materials and calculation model are used as in Fig. 9, except that the thickness of the layer is $H/a=$4.0.
By comparing the results in Fig. 12 and Fig. 13 with results in Fig. 8 and Fig. 9 for $\eta =$2.0, it can be seen that the values and spatial distributions of the displacement (pore pressure) amplitudes are different for all three modes, which may be attributed to the changed dynamic interaction between valley and layered site (due to changes in the thickness of the layer) and also to the different values of the apparent velocity for the two layered halfspace (Table 2). It can also be seen from Fig. 12 and Fig. 13 that the values and spatial distributions of the displacement (pore pressure) amplitudes are very similar for the second and the third modes because of the closer apparent velocities of the two modes (Table 2).
Fig. 12Surface displacement amplitudes for the first three modes (GR/GL= 5, H/a= 4.0)
Table 2Apparent velocities for the first three Rayleigh modes
${G}^{R}/{G}^{L}=$5, $H/a=$2  ${G}^{R}/{G}^{L}=$2, $H/a=$ 2  ${G}^{R}/{G}^{L}=$ 5, $H/a=$ 4  
$\eta =$ 1.0  $\eta =$ 2.0  $\eta =$ 1.0  $\eta =$ 2.0  
Mode1  1309.88+67.76$i$  1301.69+64.94$i$  1308.56+67.17$i$  1301.69+64.94$i$ 
Mode2  2226.44+156.01$i$  1511.82+89.30$i$  2085.67+133.77$i$  1431.49+73.30$i$ 
Mode3  3524.38+388.14$i$  2413.09+149.86$i$  1517.60+99.14$i$  
For layered halfspace of ${G}^{R}/{G}^{L}=$2 and $H/a=$ 2, there are no higher modes than the second mode at the frequency $\eta =$1.0 
Fig. 13Contours of pore pressure amplitudes for the first three modes (GR/GL= 5, H/a= 4.0)
5. Conclusions
The 3D scattering of obliquely incident Rayleigh waves by a saturated alluvial valley embedded in a layered halfspace is studied using the indirect element method in the frequency domain. Results from the analysis indicated that as the porosity decreased, the pore pressure amplitude values decreased significantly, while the oscillations became violent. The values and spatial distributions of the displacements amplitudes and of the pore pressure varied with the Rayleigh modes. The higher Rayleigh modes resulted in longer wavelengths and lead to simpler spatial distributions. The responses around the valley are determined by the dynamic interaction between valley and layered halfspace, and the dispersion characteristics of the Rayleigh waves for a given layered halfspace.
References

Aki K., Larner K. L. Surface motion of a layered medium having an irregular interface due to incident plane SH waves. Journal of Geophysical Research, Vol. 75, 1970, p. 19211941.

Trifunac M. D. Surface motion of a semicylindrical alluvial valley for incident plane SH waves. Bulletin of the Seismological Society of America, Vol. 61, 1971, p. 17551770.

Wong H. L., Trifunac M. D. Surface motion of a semielliptical alluvial valley for incident plane SH waves. Bulletin of the Seismological Society of America, Vol. 64, 1974, p. 13891408.

Yuan X., Liao Z. Scattering of plane SH waves by a cylindrical alluvial valley of circulararc crosssection. Earthquake Engineering and Structural Dynamics Journal, Vol. 24, 1995, p. 13031313.

Todorovska M., Lee V. W. Surface motion of shallow circular alluvial valleys for incident plane SH waves: analytical solution. Soil Dynamics and Earthquake Engineering, Vol. 4, 1991, p. 192200.

Dravinski M. Scattering of plane harmonic SH wave by dipping layers or arbitrary shape. Bulletin of the Seismological Society of America, Vol. 73, 1983, p. 13031319.

Dravinski M., Mossessian T. K. Scattering of plane harmonic P, SV, and Rayleigh waves by dipping layers of arbitrary shape. Bulletin of the Seismological Society of America, Vol. 77, 1987, p. 212235.

Liang Ba J. W. Z. N. Surface motion of an alluvial valley in layered haftspace for incident plane SH waves. Earthquake Engineering and Engineering Vibration, Vol. 27, 2007, p. 19, (in Chinese).

SainchezSesma F. J., RamosMartinez 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 and Structural Dynamics Journal, Vol. 22, 1993, p. 279295.

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, Vol. 79, 1989, p. 13611382.

Boore D. M., Larner K. L., Aki K. Comparison of two independent methods for the solution of wave scattering problems: response of a sedimentary basin to incident SH waves. Journal of Geophysical Research, Vol. 76, 1971, p. 558569.

Bouchon M. Effect of topography on surface motion. Bulletin of the Seismological Society of America, Vol. 63, 1973, p. 615632.

Biot M. A. The theory of propagation of elastic waves in a fluidsaturated porous solid I: lowfrequency range. The Journal of the Acoustical Society of America, Vol. 28, 1956, p. 168178.

Biot M. A. The theory of propagation of elastic waves in a fluidsaturated porous solid II: Higherfrequency range. The Journal of the Acoustical Society of America, Vol. 28, 1956, p. 179191.

Biot M. A. Mechanics of deformation and acoustic propagation in porous media. Journal of Applied Physics, Vol. 33, 1962, p. 14821498.

Zhou X., Jiang L., Wang J. Scattering of plane wave by circulararc alluvial valley in a poroelastic halfspace. Journal of Sound and Vibration, Vol. 318, 2008, p. 10241049.

Liu Z., Liang J. Diffraction of plane P waves around an alluvial valley in poroelastic halfspace. Earthquake Science, Vol. 23, 2010, p. 3543.

Li W., Zhao C. Scattering of plane SV waves by circulararc alluvial valleys with saturated soil deposits. Soil Dynamics and Earthquake Engineering, Vol. 25, 2005, p. 981995.

Pedersen H., Sanchez Sesma F.J., Campillo M. Threedimensional scattering by twodimensional topographies. Bulletin of the Seismological Society of America, Vol. 84, 1994, p. 11691183.

De Barros F. C. P., Luco J. E. Amplication of obliquely incident waves by a cylindrical valley embedded in a layered halfspace. Soil Dynamics and Earthquake Engineering, Vol. 14, 1995, p. 163175.

Liang J., Wei X., Lee V. W. 3D scattering of plane SV waves by a circulararc alluvial valley. Chinese Journal of Geotechnical Engineering, Vol. 31, 2009, p. 13451353, (in Chinese).

Wolf J. P. Dynamic SoilStructure Interaction. Englewood Cliffs, PrenticeHall, NJ, 1985.

Liang J., You H. Green’s functions for uniformly inclined line distributed loads acting on an inclined line in a poroelastic layered site. Earthquake Engineering and Engineering Vibration, Vol. 4, 2005, p. 233241.

Liang J., You H., Lee W. W. Scattering of SV waves by a canyon in a fluidsaturated, poroelastic layered halfspace, modeled using the indirect boundary element method. Soil Dynamics and Earthquake Engineering, Vol. 26, 2006, p. 611625.

Ba Z., Liang J. 2.5D Scattering of incident plane SV waves by a canyon in layered halfspace. Earthquake Engineering and Engineering Vibration, Vol. 9, 2010, p. 587595.
About this article
This study is supported by National Natural Science Foundation of China under Grant No. 51578373 and 51578372, and Tianjin Research Program of Application Foundation and Advanced Technology under Grant No. 12JCQNJC04700, which is gratefully acknowledged.