3-D scattering of obliquely incident Rayleigh waves by a saturated alluvial valley in a layered half-space

2.1. Biot’s theory

The equations of motion of a poroelastic medium [13-15] can be expressed as:

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 =(1-n){\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 density-like 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:

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:

where:

and $A_{P1}$, $B_{P1}$, $A_{P2}$, $B_{P2}$, $A_{SV}$, $B_{SV}$, $A_{SH}$ and $B_{SH}$ are the amplitudes of the up-going and down-going $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:

2.2. 3-D stiffness matrix of saturated, layered half-space and solutions of the free field responses

The constitutive equations for a poroelastic medium have the following form [13-15]:

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 3-D dynamic stiffness matrix of the saturated soil layer ${\mathbf{S}}_{P1-P2-SV-SH}^L$. Applying loads at the free surface of the saturated half-space, only down-going 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 up-going 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 half-space matrix ${\mathbf{S}}_{P1-P2-SV-SH}^R$ can be obtained. Considering the continuity of displacements and equilibrium of stresses and pore pressure at the layer interfaces, the global 3-D stiffness matrix of a multi-layered site ${\mathbf{S}}_{P1-P2-SV-SH}$ are obtained by assembling the layer matrices ${\mathbf{S}}_{P1-P2-SV-SH}^L$ and the half-space stiffness matrix ${\mathbf{S}}_{P1-P2-SV-SH}^R$. Once ${\mathbf{S}}_{P1-P2-SV-SH}$ is obtained, the discretized dynamic-equilibrium 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 free-field responses of the incident Rayleigh waves in a layered half-space, the relationship between the apparent velocity $c$ and the frequency $\omega$ must be established. Using the 3-D dynamic stiffness matrix of the layered site, the relationship between $c$ and $\omega$ can be established by setting $\left|{\mathbf{S}}_{P1-P2-SV-SH}\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 dynamic-stiffness matrix ${\mathbf{S}}_{P1-P2-SV-SH}$, which corresponds to the so-called 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 up-going waves and down-going waves in the layers can be determined using Eq. (8a). Then the free-field 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 free-field 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 free-field 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 half-space (Fig. 3). These versions of green’s functions for a layered half-space 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):

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:

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 half-space 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 half-space 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 half-space 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:

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 half-space and for the valley with an undrained free surface, respectively.

2.4. Boundary conditions

If the free surface of the layered half-space and the interface of the valley have a drained boundary condition, the boundary conditions at the valley interface $S$ can be expressed as:

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:

where:

Combining the solutions of Eqs. (15) and (12), the displacements and stresses outside of the valley will be:

and the displacements and the stresses inside of the valley will be:

If the free surface of the layered half-space 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]

4.1. Responses of a valley in a uniform half-space

A semi-circular saturated valley embedded in a uniform fluid-saturated half-space is first analyzed for both drained and undrained boundary conditions. The material parameters of the half-space 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 half-space, 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 free-field horizontal displacements of the obliquely incident Rayleigh waves.

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 3-D responses (i.e., for $\psi \ne$ 90°) are different from the 2-D case (i.e., for $\psi =\mathrm{}$90°). These results indicate that the obliquely incident wave fields cannot be simply decomposed into in-plane wave fields and anti-plane wave fields, and then be calculated as the total in-plane and anti-plane 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 2-D case ($\psi =$90°), only the in-plane 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 half-space. 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 half-space 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 half-space

A semi-circular valley embedded in a single elastic layer overlying an elastic half-space 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³, 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 half-space are as follows: ${\rho }^R=$1855 KN/m³, shear modulus ratio of the half-space 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 free-field horizontal displacements of the obliquely incident Rayleigh waves in the layered half-space.

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 half-space 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. 8-9 and in Figs. 10-11, 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 half-space (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)

Fig. 13Contours of pore pressure amplitudes for the first three modes (GR/GL= 5, H/a= 4.0)