Surface Motion for P-Wave Scattering by an Arbitrary-Shaped Canyon in Saturated Half-Space

: This paper obtained a semianalytical solution for the P-wave scattering problem by an arbitrary-shaped canyon in a saturated half-space by using Biot’s theory, the wave function expansion method, and the moments method. Firstly, based on the Biot ﬂuid-saturated porous media theory and the wave function expansion method, the wave potentials which automatically satisfy the zero-stress boundary condition on the surface of the half-space are obtained. Then, the boundary value problem is transformed into an algebraic problem by the method of moments according to the boundary conditions, and then solved numerically by truncation. By adjusting the parameters, the saturated medium in the original model approximately degenerates into a single-phase elastic medium, and the correctness of the proposed method is veriﬁed by comparing it with published results. Finally, the effects of the canyon shape, the porosity of the soil, and the incidence angle and frequency of the incident wave on the amplitude of ground surface motion are investigated. The results show that the incidence angle has a signiﬁcant effect on the ground surface motion, while the porosity of the soil has little inﬂuence on the amplitude of surface motion. The inﬂuence of canyon shape on surface motion is mainly reﬂected in the shielding effect on the incident wave. Although P-wave scattering by a canyon is a traditional problem, most of the analytical solutions are limited to solving the scattering of seismic waves in geometric regular canyons. However, the actual shape of the canyon is not regular, which limits the application of closed analytical solutions in practical engineering. In this paper, the scattering of P-waves in arbitrary-shaped canyons is successfully solved by using a semianalytical method combined with a numerical method-moments method, which provides a possibility for engineering application.


Introduction
The scattering of seismic waves by canyon topography in half-space is a classic model for studying site seismic effects [1,2]. In the past, over half a century ago, a large number of scholars engaged in this research area and achieved rich research results [3][4][5]. These results provide a solid theoretical basis for the study of ground motion characteristics and structural seismic design.
The analytical method simplifies the physical model and uses the wave function expansion method to solve the analytical solution of the wave scattering problem. For SH-waves, Trifunac et al. first used the wave function expansion method [6] to obtain the analytical solution to the SH-wave scattering problem of a semicircular canyon [3] in half-space, and analyzed the influence of wavelength, incident angle, and other factors on canyon response. Since then, different scholars have carried out further research on seismic wave scattering according to its terrain (V-shape canyon [7,8], semicircular canyon group [9], and a radially multilayered semicylindrical canyon with a rockfill dam [10]). As the research deepens, some scholars have begun to pay attention to the inhomogeneity of the soil. Zhang et al. [11] derived an exact solution for a surficially inhomogeneous semicylindrical canyon under oblique incident SH-waves. Panji [12] presented the solution of a linear corrugated orthotropic stratum above a semi-infinite substratum based on a half-space time-domain boundary element method (TD-BEM) subjected to oblique incident plane SH-waves. Similar to Panji [12], based on TD-BEM, Mojtabazadeh-Hasanlouei et al. [13] successfully obtained the surface motion of an orthotropic Gaussian-shaped sedimentary basin embedded in a linear elastic half-space under transient SH-wave propagation. Parvanova et al. [14] deals with seismic response evaluation of a finite multilayered geological structure with nonparallel boundaries rested on a semi-infinite elastic half-plane subjected to transient seismic waves.
Compared to SH-waves, the analytical solution for P/SV wave scattering by local terrain is complicated due to the coupling of P and SV wave propagation, making it challenging to satisfy stress-free boundary conditions on the half-space surface [15,16]. For this problem, many scholars have provided relevant methods. Lee and Cao [4,5] were the first to approximate the horizontal half-space surface by introducing a circular surface with a sufficiently large radius (namely the big-arc method). They obtained an analytical solution for the scattering problem of P and SV waves by circular canyons by combining this method with the wave function expansion method [6]. The essence of this method is to transform the wave problem in the infinite region into a finite region, thus avoiding the problem of satisfying the zero-stress condition on the horizontal ground surface. Zhang et al. [16] applied the complex variable function theory [17] and the mirror image method [6] to map the horizontal surface into a semicircle, solve the horizontal surface stress-free boundary conditions on the image plane, and then obtain the analytical solution of the plane wave scattering problem in a semicircular canyon. Yao et al. [18] proposed a wave function combination method based on the Lamb solution [19] and the wave function expansion, and obtained the series solutions of P, SV, and Rayleigh wave propagation and scattering around a semicircular canyon. Lee et al. used the wave potential function [15] that automatically satisfies the surface stress-free boundary condition, and combined it with the method of moments [20] to solve the P-wave scattering problem by an almost circular arbitrary-shaped canyon [21], a shallow arbitrary-shaped canyon [22], and an arbitrary-shaped alluvial valley [23] in elastic half-space. With the development of research, many scholars have extended the scattering problem of P/SV waves from uniform elastic media to saturated porous elastic media. Based on Biot theory [24][25][26] and the bigarc method [4,5], Li et al. [27] proposed an approximate analytical solution of SV wave scattering by a cylindrical canyon in a saturated elastic half-space. Liang et al. [28] obtained an analytical solution to the scattering problem of SV waves by a shallow circular-arc canyon in a saturated elastic half-space. Satyam et al. [29] applied Biot's theory to simulate the earthquake response of a saturated soil site and develop hazard maps against seismic hazards for Vijayawada city, India. Zhang et al. [30] applied the Hankel integral transform method to solve the analytical solution of the P-wave scattering problem in a semicircular canyon in a complex plane. Some researchers have begun to pay attention to seismic wave propagation in partially saturated soil; Barak et al. [31] investigated inhomogeneous plane harmonic wave propagation in dissipative partially saturated soils. Then, Kumari et al. [32] presented the procedure for the exact evaluation of the horizontal and vertical motion at the surface of partially saturated soils induced by faster waves. With the development of ocean engineering, some scholars began to pay attention to the scattering of seismic waves by canyons in the sea. Martínez-Calzada et al. [33] showed numerical results of seismic amplification on the seabed where there is the presence of deposits formed by soft soils. Feng et al. [34] studied the effects of underwater canyons on plane seismic waves (SV wave) propagation.
It is worth pointing out that the analytical solution is limited in that it can only solve some terrain with regular geometric shapes, such as circular arc-shaped canyons [27,28,30]. At present, there is little research on the analytical solution of seismic wave scattering by complex terrain in saturated half-space. In this paper, based on Biot fluid-saturated porous media theory [24][25][26], the wave function expansion method [6], and the moments method [20], a semianalytical solution is obtained for the P-wave scattering problem of an arbitrary-shaped canyon in saturated half-space. Then, the effects of canyon shape, the porosity of the soil, and the incident angle and frequency of the incident wave on the site seismic response are analyzed, which can provide a theoretical basis for seismic design of engineering construction of canyon sites in saturated half-space. Figure 1 presents a computational model for solving the seismic effects of an arbitraryshaped canyon site under P-wave incidence. As shown in Figure 1, Ω is the half-space region that is a uniform, isotropic, and fluid-saturated porous elastic medium. λ and µ are Lame's constants, and ρ is the total density of the fluid-saturated medium. There is a plane P-wave with incident angle γ 1 (the angle between the propagation direction and the horizontal direction) and the circular frequency ω, which causes the canyon response. The width of the canyon is 2a. The cartesian coordinate system x-o-y and the polar coordinate system r-o-θ are introduced, where o is the midpoint of the canyon. There exists any point on the boundary wheren is the normal line passing through the point, and α is the angle between the normal line and the polar diameter r passing through the point. According to the definition in [21], when r turns clockwise ton, α is negative; when r turns counterclockwise ton, α is positive.

Biot's Theory
Based on Biot's theory [24][25][26], without considering the viscous coupling of solid and liquid two-phase media in saturated soil and the dissipative characteristics of the media, the wave equations of skeleton and fluid in saturated elastic media can be written as follows [35] (λ + 2µ) where φ and Φ are P-wave potential functions in soil skeleton and fluid, and ψ and Ψ are SV wave potential functions in soil skeleton and fluid. λ, µ, P, and Q are the elastic constants of the saturated porous media, which are determined by Poisson's ratio v, porosity n s , critical porosity n cr , critical bulk modulus K cr , solid particle bulk modulus K g , fluid bulk modulus K f , and soil skeleton bulk modulus K dry . K dry = Kcr + (1 − n s /n cr )(K g − Kcr), µ = 3(1 − 2v)K dry /(2 + 2v), Q = n s (1 − n s − K dry /K g )K g /((1 − n s − K dry /K g ) + n s K g /K f ), P = n s 2 K g /((1 − n s − K dry /K g ) + n s K g /K f ), λ = 3v/(1 + v)K dry + Q 2 /P. ρ 11 , ρ 12 , and ρ 22 are the dynamic mass coefficients. ρ 11 = (1 − n s )ρ s +γ(1 − n s )ρ f , ρ 12 = −γ(1 − n s )ρ f , ρ 22 = n s ρ f +γ(1 − n s )ρ f , where ρ s is the density of dry soil; ρ f is the fluid density;γ is the interaction coefficient between the reaction soil skeleton and the fluid, depending on the shape of the soil particles, when the soil particles are spherical,γ = 0.5 [35].
By solving Equations (1) and (2), we can derive the wave equations about φ 1 , φ 2 , and ψ t [36,37] In the above formula, φ 1 and φ 2 are longitudinal waves, corresponding to the potentials of fast and slow P-waves, respectively; ψ t is the transverse wave, representing the potential of the SV wave. k 1 , k 2 , and k t are the fast P-wave number, slow P-wave number, and SV wave number, respectively, k i = ω/c i (I = 1, 2), k t = ω/c t . c 1 and c 2 are fast P-wave velocity and slow P-wave velocity, respectively, and c t is SV wave velocity [35].
where A = (λ + 2µ)P − Q 2 , B = ρ 11 P + ρ 22 (λ + 2µ) − 2ρ 12 Q, C = ρ 11 ρ 22 − ρ 2 12 . The P-wave and SV wave potentials in the soil skeleton are expressed as follows and the P-wave and SV wave potentials in the fluid can be expressed as where The stress-strain relationship of fluid-saturated, poroelastic media is as follows [36] where σ ij and e ij are the effective stress and strain tensors of the soil skeleton, respectively; δ ij is the Kronecker function; σ p is fluid stress.

Boundary Conditions
In this problem, the stress-free boundary condition is exhibited along the horizontal surface of half-space. In the cartesian coordinates x-o-y and polar coordinates r-o-θ, the stress boundary conditions are expressed as follows in Equations (14) and (15), respectively. σ yy = 0, τ xy = 0 and σ p = 0 at y = 0 (14) At the canyon boundary, the stress-free boundary condition also needs to be satisfied, and the fluid stress for the drained case will vanish. Thus, the stress and pore pressure boundary conditions result in the form As shown in Figure 1, T r is the radial stress component; T θ is the angular stress component [21] T

Free-Field
Plane longitudinal fast P-wave φ (i) 1 enters the half-space with circular frequency ω, incident angle γ 1 , and amplitude 1. Under the rectangular coordinate system x-o-y, the incident plane fast P-wave is characterized as When the incident wave propagates to the horizontal surface, it can produce a reflected fast P-wave φ (r) t . The expression of its potential is as follows in which γ 2 and γ t are the angle between the reflected slow P-wave and the reflected SV wave and the horizontal line, respectively, which satisfies Snell's Law k 1 cosγ 1 = k 2 cosγ 2 = k t cos γ t . A P 1 , A P 2 , and A P t are the amplitudes of the corresponding reflected fast and slow P-waves and SV waves, respectively.

Scattering Waves
According to [21], the scattering waves of the soil skeleton in half-space can take the following forms In the above formulas, φ t are the potentials of the scattered fast P-wave, slow P-wave, and the scattered SV wave, respectively, and A n , B n , and C n are the coefficients to be determined.
The total wave field in saturated half-space can be obtained as follows

Expression of Stress and Displacement
For the horizontal surface boundary condition, according to [15,21], the stress-free boundary condition can be automatically satisfied by using the expression of the wave field potential function containing only the sine function. Therefore, the problem can be solved only by solving the undetermined coefficient of the scattered wave field satisfying the stress-free boundary condition on the canyon boundary.
The displacement in the soil skeleton is expressed as follows The effective stress in the soil skeleton σ r , σ θ , and τ rθ are calculated by submitting Equations (35)-(37) into Equation (12), and is shown as In addition, by submitting Equations (35) and (36) into Equation (13), the fluid stress σ p can be gained as follows

The Solution
According to the above explanation, submitting the stress expression Equations (40)-(43) contained the undetermined coefficients into the corresponding boundary condition Equations (16)-(18), we can obtain where X 1 n = A n , X 2 n = B n , X 3 n = C n ; the expressions of L ij n and R j are shown in Appendix A. To solve Equation (44), the moments method [20] is applied. Multiply both sides of Equation (44) by the weight function sin(mθ) [21], and then integrate on [0, π] to obtain Thus, a set of infinite linear equations Equation (45) with undetermined coefficients is obtained. Solving this system of equations gives A n , B n , and C n . After the unknown coefficient is determined, the wave potentials can be obtained, and the surface motion can be solved.

Verification
When solving Equation (45), it is necessary to truncate the infinite series to obtain the convergent surface motion result. The series expansion term n and the number of weight functions m need to be tested in detail when truncating, and the test method can be seen [21]. Define the truncation term of the number of series expansion terms n of the wave potential function as N, and the truncation term of the number of weight functions m as M.
To ensure that the coefficient matrix obtained is square, M = N is set in this paper.
To consider the influence of canyon width and the wavelength Λ, a dimensionless frequency is defined [3] To convert the displacement of the soil skeleton from u r and u θ in the polar coordinate system to u x and u y in the rectangular coordinate system, the conversion formula in [21] is applied as follows In addition, the displacement in the analysis will be normalized. For an incident P-wave with an amplitude of 1, the normalized displacement is defined as follows The correctness of the calculation is verified by comparing the degradation model with published results. Take the parameters n cr = 0.36, K cr = 200 Mpa, K g = 36,000 Mpa, K f = 2000 Mpa, ρ g = 2650 kg/m 3 , ρ f = 1000 kg/m 3 , v = 0.25 [35]. To degrade this arbitraryshaped canyon model in a saturated half-space to a semicircular canyon model in a singlephase medium [21], the porosity of saturated soil is taken as a small value, n s = 0.01, in which the saturated soil can be approximately considered as a single-phase medium soil. Taking the dimensionless frequencies of 2 and 4, and considering both vertical and oblique incidences, the comparative results are illustrated in Figure 2. The degenerated solution of this paper has a good coincidence with the model [21], which verifies the correctness of the calculation in this paper.

Parameter Analysis
The method used in this paper can be used to solve the scattering problem of seismic waves by an arbitrary-shaped canyon. In this section, we will take semicircular and trapezoidal canyons as examples to research the effects of canyon shape, porosity n s , incidence angle γ 1 , and dimensionless frequency η on the ground surface motion. The parameters analyzed in the two cases are as follows: incidence angles γ 1 = 5 • , 30 • , 60 • , and 90 • , where 90 • is the vertical incidence and 5 • is the approximate horizontal incidence; the porosity of saturated soil n s = 0.1, 0.3, 0.34, corresponding to 'hard soil', 'medium hard soil', and 'soft soil', respectively [35]; the dimensionless frequencies η = 1, 2, 4, and 6, representing low frequency, mid frequency, mid-high frequency, and high frequency, respectively. For each canyon shape and dimensionless frequency, the number of truncated terms N in the equations is tested when calculating the surface displacement response. The value of N is continuously increased until the result converges or numerical stability is reached, and the test method can be seen [21]. In general, the parameters are discussed from three angles: incident wave property (incidence angle and dimensionless frequency), canyon site property (canyon shape), and half-space site property (porosity of saturated soil). As shown in Figure 3, the cases of semicircular canyon Figure 3a and trapezoidal canyon Figure 3b are calculated, respectively. For the two cases, plot the ground surface displacements of the xand y-components in the range of x/a ∈ [−4, 4]. For convenience, the whole drawing region is divided into three categories: the front side (in which the waves arrive first x/a ∈ [−4, −1]), the backside (x/a ∈ [1,4]), and the surface of the canyon (x/a ∈ [−1, 1]). Please see Section 4 for calculation parameters.

Semicircular Canyon
The semicircular canyon is shown in Figure 3a. The width of the canyon is 2a, and a is the radius of the semicircular canyon. The number of truncated terms for each dimensionless frequency is N = 50, 76, 86, and 102.

Medium Porosity
For dimensionless frequency η = 2, Figure 4 shows the effects of incidence angle and porosity on ground surface displacement response. The influence of porosity on the amplitude of x-component and y-component displacement is different, but the peak value corresponding to the porosity is the same. Compared with y-component displacement, the amplitude of x-component displacement is more affected by porosity, especially on both sides of the canyon, and the amplitude of surface displacement decreases with the increase of porosity on the whole. It should be noted that the surface x-component displacement peak corresponds to the porosity. Except for the incidence γ 1 = 5 • , the peak displacement occurs when the porosity n s = 0.1, and the peak displacement occurs when n s = 0.34 at all other incidence angles. The influence of porosity on y-component displacement is not so obvious, especially at the front side and backside of the canyon; the y-component displacement amplitude curves are almost similar under different porosity, but the influence of porosity on the canyon boundary is more obvious. However, the porosity corresponding to the peak value of y-component displacement is consistent with that of x-component displacement.  Figure 5 shows the influence of different incidence angles and dimensionless frequencies on surface displacement when the porosity ns = 0.3. The incidence angle has a significant influence on the amplitude of surface displacement. When the P-wave is oblique incident at the angle of γ1 = 30°, the amplitude of x-component displacement is the largest, followed by 60°, 5°, and 90°. For y-component displacement, the amplitude of  Figure 5 shows the influence of different incidence angles and dimensionless frequencies on surface displacement when the porosity n s = 0.3. The incidence angle has a significant influence on the amplitude of surface displacement. When the P-wave is oblique incident at the angle of γ 1 = 30 • , the amplitude of x-component displacement is the largest, followed by 60 • , 5 • , and 90 • . For y-component displacement, the amplitude of displacement is the largest when the P-wave is vertically incident at the incident angle γ 1 = 90 • , followed by 60 • , 30 • , and 5 • . At the canyon boundary, the difference between the effects of the incident angle on the horizontal displacement amplitude is not obvious, and the effect of the incident angle on the y-component displacement amplitude is the same as that on both sides, but the displacement oscillation of the two components is more obvious than that on both sides, and the displacement peak value appears at near of x/a = ±1. In addition, at oblique incidence, the displacement curve of the backside of the canyon is relatively stable and the oscillation is weak, which reflects the shielding effect of the canyon.  Figure 3b shows the trapezoidal canyon model. The upper base span of the canyon is 2a, and the height of the trapezoid is h, in this case, h/a = 1. The slope angle of the trapezoid is the angle between the waist and the horizontal direction, and δ = 60 • . The number of truncated terms for each dimensionless frequency is N = 20, 28, 32, and 36. Figure 6 shows the effects of different incidence angles and porosity on surface displacement in the case of dimensionless frequency η = 2. The influence of porosity on the amplitude of x-component and y-component surface displacement is different. Compared with y-component displacement, the amplitude of x-component displacement is more affected by porosity, and decreases with the increase of porosity. It should be noted that when the oblique incidence occurs, the peak surface displacement occurs near x/a = −1, and the porosity is n s = 0.1. The influence of porosity on y-component displacement is not so obvious, especially at the front side and backside of the canyon; the y-component displacement amplitude curves are almost similar under different porosity, but the influence of porosity on the canyon boundary is more obvious. The position of the y-component displacement peak value is consistent with that of the x-component displacement, and the porosity at this time is n s = 0.1. placement in the case of dimensionless frequency η = 2. The influence of porosity on the amplitude of x-component and y-component surface displacement is different. Compared with y-component displacement, the amplitude of x-component displacement is more affected by porosity, and decreases with the increase of porosity. It should be noted that when the oblique incidence occurs, the peak surface displacement occurs near x/a = −1, and the porosity is ns = 0.1. The influence of porosity on y-component displacement is not so obvious, especially at the front side and backside of the canyon; the y-component displacement amplitude curves are almost similar under different porosity, but the influence of porosity on the canyon boundary is more obvious. The position of the y-component displacement peak value is consistent with that of the x-component displacement, and the porosity at this time is ns = 0.1.  Figure 6. Influence of porosity on the displacement response of a trapezoidal canyon (η = 2). Figure 6. Influence of porosity on the displacement response of a trapezoidal canyon (η = 2). Figure 7 shows porosity n s = 0.3 in the saturated half-space, and the influence of different incidence angles and dimensionless frequencies on the surface displacement. The oscillation of the displacement amplitude curve increases obviously with the increase of dimensionless frequency, and the incidence angle has a significant effect on the surface displacement. Compared with the low frequency η = 1, the amplification effect of the medium and high frequency η = 4 and 6 on the displacement amplitude is significantly increased, and the oscillation peak value appears near x/a = ±1. The incidence angle has an obvious influence on the x-component displacement amplitude, and the displacement amplitude is the largest when γ 1 = 30 • , followed by 60 • , 30 • , and 5 • . Consistent with horizontal displacement, the influence of incident angle on vertical displacement amplitude is also obvious. When γ 1 = 90 • , the vertical displacement amplitude is the largest, followed by γ 1 = 60 • , 30 • , and 5 • . For the cases of medium and high frequency η = 4 and 6, when the oblique incidence γ 1 = 5 • , 30 • , and 60 • , the oscillation of the displacement amplitude curve of the backside is significantly smaller than that of the front side of the canyon, which obviously reflects the shielding effect of the canyon. In addition, at the same dimensionless frequency, the oscillation of the displacement curve is more obvious at the boundary of the trapezoidal canyon than that of the semicircular canyon.

Conclusions
First, based on Biot's fluid-saturated porous media theory, this paper applies the scattering wave field potential function, which can automatically satisfy the surface stress freedom in half-space, and then converts the boundary value problem into a series of algebraic equations with infinite terms by applying the moments method according to the stress-free boundary conditions on the boundary of an arbitrary-shaped canyon. The semianalytic solution of P-wave scattering by arbitrary-shape canyon in saturated halfspace is obtained by the numerical solution of truncation. On this basis, the influence of the shape of the canyon, porosity, incidence angle, and frequency on surface displacement is investigated. The main conclusions are:

1.
The incidence angle has an obvious influence on the displacement amplitude. For the solved cases, almost all of them meet the following rules: for x-component displacement, when the P-wave is oblique incident at the incidence angle γ 1 = 30 • , the displacement amplitude is the largest, followed by 60 • , 5 • , and 90 • ; for y-component displacement, the amplitude of displacement is the largest when the P-wave is vertically incident at the incident angle γ 1 = 90 • , which is 60 • , 30 • , and 5 • , successively.

2.
Compared with the incidence angle, porosity has little influence on the amplitude of displacement, especially for y-component displacement. The displacement curves of the three kinds of porosity almost coincide on the front side and backside of the canyon, but there is a slight difference in the inner part of the canyon. The effect on the displacement amplitude in the horizontal direction is more significant, and the displacement amplitude is smaller and smaller as the porosity increases on the whole. 3.
The influence of the shape of the canyon on displacement amplitude is mainly reflected in the shielding effect on the incident wave. At the same dimensionless frequency, the ground surface displacement of the trapezoidal canyon is more oscillating than that of the semicircular canyon, indicating that the trapezoidal canyon has a more shielding effect on the incident wave.

Data Availability Statement:
The raw data supporting the conclusions of this article will be made available by the authors.

Conflicts of Interest:
The authors declare that they have no known competing financial interest or personal relationships that could have appeared to influence the work reported in this paper.