## Abstract

This study focuses on the theoretical aspect of topographic scattering induced by a shallow asymmetric V-shaped canyon under plane shear horizontal-wave incidence. An analytical approach, based on the region-matching technique, is applied to derive a rigorous series solution, which is more general than that in a previous study. For the wave functions constrained in two angular directions, a novel form of Graf's addition formula is derived to arbitrarily shift the local coordinate system. Barrier geometry, angle of incidence and wave frequency are taken as the most significant parameters in exploring the topographic effects of localized concave free surfaces on ground motions. Both surface and subsurface motions are presented. Comparisons with previously published results and boundary-element solutions show good agreement. Frequency-domain results indicate that, for the high-frequency case at a low grazing angle (corresponding to the potential case in teleseismic propagation), the high levels of amplified motions occur mostly on the illuminated side of the canyon. When the windward slope is steeper, the peak amplitude values, at least 2.4 times larger than those of free-field responses, tend to increase. Time-domain simulations display how a sequence of scattered waves travel and attenuate at regional distances.

## 1. Introduction

Topography-induced scattering has been known to affect the propagation direction of wavefields at regional distances [1–3]. The area of rugged terrain may intensify diversity of seismic motions and give rise to the generation of complex seismic signals. Consequently, the term ‘topographic effect’, which is closely related to the localized geometric character of the free surface (e.g. canyons, slopes, cliffs and ridges), is in wide use among much of the literature [4–7]. In various branches of applied science, such as shallow geophysics, seismology and civil engineering, the measurement and assessment for potential topographic effects on earthquake ground shaking have been of continued interest [8–11].

V-shaped canyons are a common type of natural terrain on the ground surface. Since the early 1970s, some simplified models have been built to approximate the cross-sectional profiles of real canyons. For shear horizontal (SH) waves, which are the horizontally out-of-plane polarized shear waves, the symmetric V-shaped model seems to be the most popular. Various numerical methods have been developed to investigate topography-induced scattering [12–18]. These studies promoted a thorough understanding of the effect of topographic irregularities on seismic ground motions. Nowadays, substantial improvement in computing performance (fast processors, parallel processing, etc.) has boosted the growth of mesh discretization techniques. Representatives of grid-based approaches are the finite-difference, finite-element, spectral-element and boundary-element methods (BEMs) [19–21]. Although these approaches are flexible, their reliability has to be strictly verified by the existing canonical solutions.

In general, for a few simple geometric shapes [22], the method of separation of variables (MSV), which is sometimes termed the method of wave function expansions in engineering seismology, is feasible to derive exact analytical solutions. For the two-dimensional half-plane problems under SH plane wave propagation, there are only two cases available so far, that is, the semicircular and semi-elliptic canyons [23,24]. However, the MSV has a quite limited range of applications because of its dependence on the property of separability of functions. The MSV cannot be employed straightforwardly even for some simple geometric cases whose boundaries partly conform to a certain separable coordinate system, such as circular sectorial canyons, U-shaped gorges, semicircular hills and semi-elliptic ridges. This means that applying the MSV to obtain the explicit closed-form expressions for expansion coefficients is nearly impossible. Therefore, the construction of a rigorous series solution may be one of the best alternatives when the exact analytical solution cannot be attempted using the MSV. Fortunately, in this respect, the use of the region-matching technique (RMT) is effective for the above-mentioned geometric shapes [25–28].

The implementation of the RMT generally involves the following steps: (i) decompose the analysed region into two (or more) subregions; (ii) expand the displacement field within each subregion in terms of adequate wave functions with unknown coefficients; (iii) unify the coordinate systems in two adjacent subregions, if necessary; and (iv) determine the unknown expansion coefficients by imposing the matching conditions. Step (i) is sometimes termed the domain decomposition and relies on the appropriate selection of auxiliary boundaries with which the MSV can be applied directly to each subregion. Step (iii) depends on the utilization of suitable coordinate-transformed relations. Nevertheless, two key points should be mentioned before using the RMT. First, step (ii) should be correctly tied in with step (i). Second, the addition theorems [29–31] should be recast in an appropriate form for specific arrangements of the polar, elliptic or spherical coordinate systems. Only when grasping simultaneously these two key points can a breakthrough in extending the ability of the RMT to cope with more problems of antiplane scattering be possible.

Chang [32] and Tsaur & Chang [33] tried to open up a new opportunity for the application of the RMT. Chang [32] concentrated on the surface/subsurface discontinuities of truncated circular types, which include the truncated semicircular canyons, the partially filled semicircular alluvial valleys, the symmetric V-shaped canyons, the underground truncated circular cavities and the underground tunnels with a truncated circular hole [34–38]. In the context of elastodynamics, Chang [32] promoted the use of the RMT, because the types of scatterers are broadened to those with contours purely composed of straight lines and those with contours consisting of circular-arc lines together with straight lines. It is noted that in most of the previous studies the boundaries of scatterers consist of a single curved line [25]. Tsaur & Chang [33] applied the RMT to derive a correct series solution for the case of circular-arc hills. Although Yuan & Liao [39] attempted to carry out a similar task, unfortunately they obtained an erroneous series solution due to an unsuitable conjunction of aforementioned steps (i) and (ii). Chang [32] and Tsaur & Chang [33] highlight the importance of the two key points for the utilization of the RMT. Numerical results obtained from the series solutions shown in Tsaur & Chang [35–37] and Tsaur *et al.* [38] validate the newly developed numerical/theoretical methodologies [40–43]. Obviously, it is worthy of further promoting the derivations of series solutions.

With regard to V-shaped canyon geometry, owing to inconformity with any one of the separable coordinate systems, the MSV cannot be applied straightforwardly to the governing wave equation to yield the complete expressions in terms of eigenfunctions (or wave functions), albeit a very simple case. Tsaur & Chang [35], Chang [32] and Tsaur *et al.* [38] used the RMT to derive the series solutions for both the shallow and deep symmetric canyons subjected to plane wave incidence. For deep symmetric canyons, Chang [32] proposed two RMTs, with two distinct auxiliary boundaries. The stress singularity at the bottom of the deep canyon is not taken into account in the first type and is incorporated inherently in the second one. Liu [44] followed the same line presented to deal with both the shallow and deep asymmetric cases. However, Liu [44] handled merely the steady states at very low frequencies. Later, Zhang *et al.* [40] applied the RMT, with an auxiliary boundary that touches the sharp bottom of the canyon, to conduct the deep asymmetric case. Their auxiliary boundary is similar to that chosen by Liu [44]. Nevertheless, the strategy of domain decomposition adopted by Liu [44] and Zhang *et al.* [40] for the deep asymmetric case cannot work anymore for the shallow asymmetric case. Zhang *et al.* [41] also derived a unified series solution to the problem of symmetric canyons. Gao & Zhang [43] investigated both the shallow and deep symmetric cases due to a line-source excitation. In order to improve the computational efficiency of coordinate transformation, Gao & Zhang [43] constructed properly the scattered wavefields by following the same ideas of domain decomposition proposed by Tsaur & Chang [35] and Tsaur *et al.* [38]. The intrinsic difficulty in applying the MSV to several specific cases of V-shaped canyons has been conquered by these authors based on the promising RMT.

In this work, we focus on the theoretical aspect of topographic scattering induced by a shallow asymmetric V-shaped canyon. The RMT is adopted to derive a series solution to the problem under plane SH wave excitation. Suitable wave functions are well used to characterize the out-of-plane motions. An adequate form of Graf's addition formula [31] is newly derived for Bessel functions of fractional order, as compared with those limited to Bessel functions of integer order [39]. Such a recast version of Graf's addition formula not only plays the leading role in converting the wave functions from local coordinate systems to global ones, but also keeps the radial and angular wave functions decoupled after coordinate transformation, which is a key feature to make the orthogonality of angular (trigonometric) functions applicable.

The present model can be seen as an improvement of the shallow asymmetric model treated by Liu [44] and as a further extension of the shallow symmetric model conducted by Tsaur & Chang [35] and Chang [32]. The newly derived form of Graf's addition formula is devoted to arbitrary offsets so that it is more general than those, with offsets along the vertical direction, proposed by Tsaur & Chang [35] and Chang [32]. Since the coordinate-transformed formula and the enclosed-subregion expression adopted herein are simpler than those in Liu [44], the equations of this study are more compact than those in Liu [44]. Consequently, the series solution derived herein enables a more effective computational performance, notably for high frequencies. Apart from the intrinsic significance of analytically solving the present boundary-value problem, the asymmetric model is more practical than the symmetric one because many natural canyons/valleys are indeed asymmetric, e.g. the Feitsui Canyon in northern Taiwan [45], the Yellow River Canyon in the northeastern Tibetan Plateau [46], the Gamperdona Valley in western Austria [47], the Braldu Valley in northern Pakistan [48], the fluvial valleys in northwest Syria [49], the Grenoble Valley in France [50] and the San Fernando Valley in southern California [51]. Actually, the asymmetric model can degenerate into the symmetric one; by contrast, it is infeasible for the latter to degenerate into the former.

Taking into account the terms shallow and deep canyons for the asymmetric V-shaped case, the classification criterion is based on the dimensionless ratio of the distance (i.e. *r*_{e} tagged in figure 1) between the canyon bottom and the centre of the canyon top to the half-width of the canyon top (i.e. *a* shown in figure 1). When the ratio is smaller than unity (i.e. *r*_{e}/*a*<1), the canyon is classified as a shallow one; otherwise, it is classified as a deep one (i.e. *r*_{e}/*a*≥1). Complying with such a classification criterion, the canyon bottom does not penetrate the auxiliary boundary when applying the RMT to the shallow asymmetric case. This feature leads to well-defined expressions of the wavefields in individual subregions. All the field solutions satisfy the governing equations and all the boundary conditions except for those at the artificial interface. Especially for the interior subregion enclosing the main part of the canyon, the field solution contains a singularity in the stress fields near the canyon bottom. Thus, at high frequencies the computed results obtained from the derived series solution remain desirably stable. This outcome echoes the fact that selecting the appropriate auxiliary boundary is of vital importance.

The main contributions of this work can be summarized as follows. First, the proposed series solution is novel *per se* because it is likely to be absent in the literature, and it enriches the catalogue of series solutions for canyon problems related to SH wave scattering. Second, the derived series solution fills the gap in previous studies, that is, those for the symmetric V-shaped case [32,35,38,41] and those for the deep asymmetric V-shaped case [40,44]. Third, the analytical nature of the derived series solution makes the high-frequency computation delightfully easy and efficient. Fourth, the computed results obtained herein are helpful to advance knowledge of the topographic effect pertaining to shallow asymmetric V-shaped canyons.

## 2. Theoretical formulations

Consider a homogeneous, isotropic, linearly elastic, semi-infinite medium (with shear modulus *μ* and shear-wave velocity *c*_{s}) bounded by the horizontal ground surface, inlaid with an infinitely long, shallow, asymmetric, V-shaped canyon (figure 1). An infinite train of plane SH waves (with an angular frequency *ω*) is incident upon this canyon at an angle *α* to the *y*-axis. The origin of global coordinate systems (*x*,*y*) and (*r*,*θ*) is set at the centre of the canyon top, while that of local coordinate systems (*x*_{1},*y*_{1}) and (*r*_{1},*θ*_{1}) is at the canyon bottom. Relative positions between the global and local coordinate systems are associated with (*r*_{e},*θ*_{e}). The half-width and depth of the canyon are *a* and *d*, respectively. The central angle of the canyon bottom is *β*_{1}+*β*_{2}, ranging from *π* to 3*π*/2. For the present shallow case, the constraint *r*_{e}/*a*<1 has to be met, and, therefore, the interior angle between two adjacent slopes is an obtuse angle.

As seen in figure 1, through introducing a semicircular auxiliary boundary *S*_{a}, the half plane is divided into two regions, an open region 1 and an enclosed region 2. In these two regions, the steady-state out-of-plane motions are required to satisfy the governing Helmholtz equations, namely,
^{2} is the two-dimensional Laplacian and *k*=*ω*/*c*_{s} is the shear wavenumber. The subscripts *j*, where *j*=1 and 2, stand for the total displacement fields in regions 1 and 2, respectively. Throughout this work, the time-harmonic factor

The zero-stress boundary conditions on the horizontal ground surface and the canyon surface are

For the half-plane medium without any surface/subsurface irregularities, the free-field displacement *u*^{F} can be expressed as a sum of the incident waves and their reflected waves from the horizontal ground surface, that is,
*ε*_{n} is the Neumann factor (equal to 1 if *n*=0 and to 2 if *n*≥1) and *J*_{n}(⋅) denotes the *n*th-order Bessel function of the first kind. Note that this expression intrinsically satisfies the traction-free condition on the horizontal ground surface (2.2).

The ‘total’ scattered field *u*^{S} in open region 1 may be separated into two parts, *u*^{S0} and *u*^{S2}, that is, *u*^{S}=*u*^{S0}+*u*^{S2}. The first part *u*^{S0} represents the scattered field that excludes the effect of region 2, and it corresponds to the scattered field derived by Trifunac [23] for the case of semicircular canyons, that is,
*n*th-order Hankel function of the second kind and the primes stand for differentiation with respect to the arguments of corresponding functions.

The second part *u*^{S2} is the contribution due to the presence of region 2, and its proper wave function is as follows:
*A*_{n} and *B*_{n} are unknown.

The displacement of the resultant wavefield *u*_{1} in the open region 1, which is the combination of free and scattered wavefields, can be expressed as
*u*_{2}, satisfying the Helmholtz equation (2.1) and the stress-free boundary conditions (2.3), is given by
*ν*=*π*/(*β*_{1}+*β*_{2}), *C*_{n} will be determined. Note that (2.9) inherently satisfies the requirement for the stress singularity existing in the vicinity of the canyon bottom (which is a salient corner).

In order to rewrite (2.9) in terms of (*r*,*θ*), the necessary coordinate transformation from (*r*_{1},*θ*_{1}) to (*r*,*θ*) is accomplished via Graf's addition formula for Bessel functions [31], which is recast in an appropriate form as follows:
*m* and *n* are integers.

Substituting (2.10) into (2.9) results in

Next, applying the orthogonal properties of sine/cosine functions to the stress continuity condition on the auxiliary boundary *S*_{a}, that is,
*π*/2,*π*/2] leads to the following relations for unknown expansion coefficients *A*_{n}, *B*_{n} and *C*_{n}:

Similarly, the enforcement of the displacement continuity condition across the auxiliary boundary *S*_{a} is essential, that is,
*π*/2,*π*/2] gives
*A*_{n} and *B*_{n}, using the Wronskian relation that involves both the Bessel and Hankel functions [53], and rearranging yields the following coupled systems of linear algebraic equations for unknown coefficients *C*_{n}:

After the infinite series in (2.22) and (2.23) are truncated properly, the expansion coefficients *C*_{n} can be evaluated by the standard matrix techniques. For numerical computations, truncating the infinite summation to a finite number of terms is necessary. The summation indices *n* and weighting indices *q* in (2.22) and (2.23) are truncated after 2*N*−1 and *N*−1 terms, respectively; consequently, (2.22) and (2.23) constitute a system of 2*N* equations with 2*N* unknowns. The number of truncation terms under consideration depends only on the requirement for accuracy. Once the coefficients *C*_{n} are found, the expansion coefficients *A*_{n} and *B*_{n} can be evaluated straightforwardly via (2.15) and (2.16).

For the specific location, the displacement amplitude |*u*| is evaluated from the expressions of corresponding wavefields (2.8) and (2.9), that is,
*u*^{F} is evaluated directly from (2.4) to obtain the total displacement amplitude of region 1.

For convenience, in the following analysis the dimensionless frequency *η* is defined as:

The canyon model considered herein can be extended to the basin model with alluvium. When the semicircular auxiliary boundary *S*_{a} is the interface between two distinct materials, only slight modifications to some parameters involving material properties (e.g. shear modulus, mass density and shear-wave velocity) are needed in the above equations. In addition, the proposed model is applicable to the linearly viscoelastic model. For time-harmonic problems, the aforementioned formulism can be straightforwardly applied to consider the effect of damping and attenuation of filling materials. However, a few real-valued parameters (e.g. wavenumbers) have to be changed into the complex-valued ones, implying that the evaluation of Bessel functions with complex arguments is indispensable. Some simple spring-dashpot models may be used to simulate the dynamic behaviour of viscoelastic media [54].

## 3. Numerical results and discussions

### (a) Convergence test

At the initial stage of calculations, to get the proper truncation value *N* in (2.22) and (2.23), a number of convergence tests are performed. It is worth emphasizing that the inner series with indices *m* (see (A 1) and (A 2)) are truncated and bounded from −*M*+1 to *M*−1 terms. These internal sums should be calculated accurately via numerical testing for their convergence, thereby leaving only one parameter (i.e. summation indices *n*) to remove the problem of relative convergence. Following extensive numerical experiments, the value of *M*=200 is adequate to produce all the graphs hereafter.

Table 1 shows the variations in displacement amplitudes and maximum relative errors when the outer truncation index *N* increases. The selected case is of *r*_{e}/*a*=0.75 and *θ*_{e}=−30°. Only the horizontally incident cases (*α*=90°) are checked because they take more truncation terms than others to reach convergence. The relative error is defined to be
*V*_{p} and *V*_{f} denote the two adjacent computed values of displacement amplitudes. For the five specific positions (*P*_{1}–*P*_{5} located at *x*/*a*=−0.8,−0.5,−0.1,0.2 and 0.5, respectively), the values of displacement amplitudes are displayed in table 1. The maximum values of relative errors are calculated from the displacement amplitudes at 800 equally spaced locations, ranging between *x*/*a*=−4 and 4.

As shown in table 1, the numerical results converge to at least four significant digits at *N*=14 for *η*=0.25 and at *N*=32 for *η*=3. Generally, more terms of *N* are required when the wave frequency increases. The convergence criterion for *N* is that the maximum relative error falls below the threshold value of 0.01%. Based on further numerical tests for other geometries, *N*=18 is adequate for *η*≤0.5, *N*=20–30 for 0.5<*η*≤2 and *N*=32–42 for 2<*η*≤4. At higher frequencies (12≤*η*≤16), *N*=54–66 is sufficient to produce reliable results.

### (b) Verification for the shallow symmetric case

When the canyon bottom is located on the *y*-axis (i.e. *θ*_{e}=0°), the asymmetric case turns into the symmetric one. Since Sánchez-Sesma & Rosenblueth [14] applied the point-source discretization and a least-squares scheme to treat the problem of symmetric V-shaped canyons, two cases shown in their fig. 9 are taken as verification examples. Furthermore, Tsaur & Chang [35] derived a series solution to the shallow symmetric V-shaped case so that the 12 cases shown in their fig. 2*a*–*d* can be used as a check. For brevity, only a few results are illustrated here. Figure 2 displays the surface displacement amplitudes |*u*| versus the dimensionless horizontal distance *x*/*a*: figure 2*a* for *η*=0.7547 and *α*=0° at *d*/*a*=0.53, and figure 2*b* for *η*=1 and *α*=60° at *d*/*a*=0.75. Specific locations on the canyon surface are displayed by a bold black line ranging between *x*/*a*=−1 and 1. As shown in figure 2, the agreement is fairly good.

### (c) Verification for the shallow asymmetric case

Taking the asymmetric canyon into account, figure 3 shows two cases selected from fig. 2.9(a) and 2.13(d) of Liu [44] for surface displacements at *η*=2: figure 3*a* for *α*=0°, *d*/*a*=0.5 and *θ*_{e}=−56.3099°, and figure 3*b* for *α*=90°, *d*/*a*=0.75 and *θ*_{e}=−33.6901°. Figure 3*a*,*b* exhibits consistency between the present results and those of Liu [44].

To reveal the sensitivity of verification checks with respect to frequency, figure 4*a*–*d* displays the computed results for *η*=6, 8, 10 and 12, respectively, when *r*_{e}/*a*=0.75, *θ*_{e}=−30° and *α*=30°. Included also are the computed results obtained from the classical BEM because only the low-frequency (*η*≤4) results are presented by Liu [44]. The adopted BEM is based on the direct formulation (e.g. Dominguez [55]) and the so-called ‘constant’ element is used. This means that unknown quantities are the nodal values of the field variables and their normal derivatives that have clear physical meanings (i.e. displacements and tractions, respectively). Employing the collocation method and enforcing the corresponding boundary conditions, a global square system of equations can be constructed. In figure 4, an overall good agreement can be found between the series solution and the boundary-element solution.

As expected, owing to the analytical nature of the adopted wave functions, the computational storage of the derived series solution is smaller than that of the classical BEM, especially for the high-frequency band. Following a general discretization criterion (i.e. around 10–20 nodal points per wavelength), a large number of elements are required in the BEM analysis to get reasonable results; for example, when *η*=12 (figure 4*d*), 1516 piecewise-constant elements are used in the BEM and *N*=54 in the derived series solution.

### (d) Surface and subsurface motions in the frequency domain

In order to demonstrate the effect of canyon asymmetry on the displacement amplitudes, figure 5 shows the computed results for five angles of incidence (*α*=0°, 30°, 90°, −45° and −90°) at *η*=2. Three canyon geometries whose bottoms are located at *d*/*a*=0.75 are selected for *θ*_{e}=0°, −15° and −35° (figure 5*a*). As displayed in figure 5*b*, the asymmetric feature of motions under symmetric excitation is expected owing to the canyon asymmetry. At vertical incidence (*α*=0°), the maximum displacement amplitude for the asymmetric case occurs on the horizontal ground surface near the steep slope of the canyon (figure 5*b*). At non-vertical incidence (*α*≠0°), the peak amplitude shifts in the same direction (to the left) when the canyon bottom shifts (figure 5*c*–*f*). At horizontal incidence (*α*=90° and −90°), the motion patterns in the illuminated region behave much like those of the standing waves (figure 5*d*,*f*). The windward slope is the left-hand surface of the canyon for 0°<*α*≤90°, and the right-hand surface for −90°≤*α*<0°. When the windward slope is steeper, the local maximum amplitude on the canyon surface tends to increase (figure 5*c*,*d*); on the contrary, it is inclined to decrease when the windward slope becomes gentler (figure 5*e*,*f*). Comparing figure 5*c*,*e* with figure 5*d*,*f*, one can see that the peak amplitudes for oblique incidence are smaller than those for grazing incidence.

In order to reveal how the surface displacement amplitude changes with *r*_{e}/*a* when *θ*_{e}=−45°, figure 6*a* displays the three canyon cross sections (*r*_{e}/*a*=0.4, 0.6 and 0.8), whereas figure 6*b*–*f* gives the results for *η*=2 at *α*=0°, 45°, 90°, −60° and −90°, respectively. In figure 6*b*,*c*, the local maximums of displacement amplitudes on the horizontal ground surface close to the steep slope of the canyon (see those within the range from *x*/*a*=−4 to −1) decrease with decreasing *r*_{e}/*a*. These local maximums tend to shift to the left when *r*_{e}/*a* becomes larger. If the wave propagation direction is parallel (or nearly parallel) to the horizontal ground surface (figure 6*d*–*f*), then the peak amplitudes of motions within the canyon increase with increasing *r*_{e}/*a*. The location of peak amplitude shifts gradually towards the position of the canyon bottom.

Figure 7 exhibits the variations in surface motions around the canyon for −90°≤*θ*_{e}≤0° at *η*=6. Corresponding curves for the maximum values of displacement amplitudes (i.e. *θ*_{e} are also included in figure 7. For the case of *θ*_{e}=90°, the canyon topography disappears so that the displacement amplitudes are equal to 2. For vertical incidence (*α*=0°) in figure 7*a*, the free-field response may be amplified 1.5–1.75 times at −68°≤*θ*_{e}≤−15°. As seen in the horizontally incident cases (figure 7*d*,*f*), the large-amplitude motions take place on the illuminated side of the canyon. In figure 7*d*, the peak values may reach 4.75 at least for −70°≤*θ*_{e}≤0°. Figure 7*f* shows that, approximately in the range of *θ*_{e} from −32° to −90°, the maximum amplitude somewhat monotonously decreases with decreasing *θ*_{e}. The overall maximum amplification of figure 7 is about 5.2 at *α*=−90° and *θ*_{e}=−32.11°.

Figure 8 displays the displacement amplitudes versus the dimensionless horizontal distance and the incident angle for *r*_{e}/*a*=0.75 at *η*=6, so as to manifest the influence of incident angles on the surface motions. Figure 8*a*–*d* corresponds to *θ*_{e}=−10°, −30°, −50° and −70°, respectively. As expected, for the nearly symmetric case in figure 8*a*, the motion pattern seems to be antisymmetric about *α*=0°. In figure 8*a*–*d* for 80°≤*α*≤90°, the peak amplitude of motions tends to increase as the incident angle bends towards the horizontal ground surface. The amplified motions are 1.5–1.75 times larger than those on the free field at −40°≤*α*≤−5° for *θ*_{e}=−10°, at −30°≤*α*≤10° for *θ*_{e}=−30° and at −20°≤*α*≤20° for *θ*_{e}=−50°. The overall maximum amplification of figure 8 is about 5.2 at *θ*_{e}=−30° and *α*=−90°.

Figure 9*a*–*f* shows the variations in spectral displacement amplitudes at the surface of a canyon with parameters *r*_{e}/*a*=0.75 and *θ*_{e}=−30°, at *α*=0°, 30°, 60°, 90°, −45° and −90°, as a function of dimensionless frequency. As shown in figure 9*a*–*f*, when the wave frequency increases, the surface motions in the illuminated region become more oscillatory than those in the shielded region. Especially for grazing incidence (figure 9*d*,*f*), one can find that the relatively high levels of ground shaking (in the red region) are on the impact side of the canyon. The peak amplitudes approach 3.5 at 5≤*η*≤16 for vertical incidence (*α*=0°) in figure 9*a* and 4.75 for horizontal incidence (*α*=90°) in figure 9*d*. For the other horizontally incident case (figure 9*f*), a similar trend can be found at 8≤*η*≤16. From the obliquely incident cases shown in figure 9*b*,*e*, the average amplification factors are between 1.2 and 1.4 at about 1≤*η*≤16. The overall maximum amplification of figure 9 is about 5.2 at *η*=5.8 and *α*=−90°.

For the high-frequency case (*η*=12), figure 10*a*–*d* demonstrates the subsurface motions around the canyon with *r*_{e}/*a*=0.75 and *θ*_{e}=−30° at *α*=0°, 45°, 90° and −90°, respectively. Figure 10*b* shows that the local maximum motions occur on the right-hand side of the horizontal ground surface (1.5≤*x*/*a*≤1.75). For grazing incidence (figure 10*c*,*d*), the subsurface motions in the illuminated region are at least two times higher than those in the nearby region. This might be attributed to the strongly constructive interference between the incident, reflected and scattered waves. The highly shadowed regions imply that the canyon is resistant enough to the wave propagation, and, therefore, it can act as an energy barrier.

Overall, the surface and subsurface motions in the vicinity of the canyon are dependent not only on the geometric shape of the canyon, but also on the frequency and incident angle of arriving waves.

### (e) Surface and subsurface motions in the time domain

Using the fast Fourier transform algorithm, one may pursue the time-domain responses in the frequency domain. The incident signal is a symmetric Ricker wavelet, which is defined to be
*f*_{c} is the characteristic frequency. In figures 11 and 12, *f*_{c} is set to be 2 Hz so that the related computations are made at discrete frequencies ranging from 0 to 8 Hz with an interval of 0.0625 Hz. The time window is 16 s. The half-width of the canyon top *a* and the shear-wave velocity *c*_{s} are 1 km and 1 km s^{−1}, respectively. The location of the canyon bottom is selected at *r*_{e}/*a*=0.75 and *θ*_{e}=−30°. The reference point for *t*=0 when the plane waves pass by is specified at the position (*x*,*y*)=(0 km, 4 km) for vertical incidence and at (*x*,*y*)=(−4 km, 0 km) for grazing incidence.

To reveal the scattering effect in the vicinity of the canyon, figure 11*a*,*b* displays the synthetic seismograms for vertical (*α*=0°) and oblique (*α*=30°) incidence, respectively; both of them contain 101 time series equally distributed along the ground surface between *x*=−4 and 6 km. On the top and bottom of figure 11*a*,*b*, several signals (see arrows *L*1, *L*2, *R*1, *R*2, etc.) received after the direct-wave signals are labelled in turn. The sources of these marked traces may be explicitly visualized in figure 12. In figure 11*a*, the direct waves impinge on the horizontal ground surface at *t*=4 s, whereas in figure 11*b* they reach *x*=−4 km at *t*=2 s and *x*=6 km at *t*=7 s. When the locations are far away enough from the canyon, the amplitudes of recorded waveforms become exactly two times higher than those of incident pulses. In figure 11*a*,*b*, the amplifications of surface motions can be up to about 1.4 times higher than those of the free-field response. Figure 11*b* shows that the shielding effect of the canyon makes the displacement amplitudes on the right-hand slope of the canyon (0 km<*x*<1 km) smaller than those of direct waves.

In order to show the discrimination between different signal sources, eight snapshots of transient subsurface motions at *α*=0° are given in figure 12. Furthermore, the arrows marked in figure 11*a* are tagged appropriately in figure 12. Figure 12*a*–*d* shows that, owing to the interaction of incoming pulses, a scattered wavetrain *B* spreads out from the canyon bottom. When two reflected waves from adjacent slopes of the canyon pass through, respectively, the upper left and right corners of the canyon, two scattered waves *L*1 and *R*1 arise and meet each other on the right-hand slope of the canyon (figure 12*d*). The re-radiated waves *L*2 and *R*2 come from the scattered waves *B*. Figure 12*e*–*h* shows that a chain of scattered waves (see *L*3, *L*4 and *R*3) is generated from the bottom and two upper corners of the canyon because they behave as new wave sources. Also, it is obvious that the amplitudes of early-generated scattered waves (*L*1 and *R*1) are greater than those of late-generated scattered waves (*L*3, *L*4 and *R*3). Moreover, a key feature of the diffraction is that a part of the reflected wavefront from the horizontal ground surface is cut off by the canyon and then continuously regenerates.

## 4. Concluding remarks

The problem of SH wave scattering induced by a shallow asymmetric V-shaped canyon has been successfully addressed and solved. The derivation of a rigorous series solution has been achieved via an analytical approach centring on the use of the robust RMT. An appropriate selection of the auxiliary boundary has promoted the application of the MSV and has also made the derived wave functions inherently include the stress singularity near the canyon bottom. A new form of Graf's addition formula has been derived for the wave functions constrained in two angular directions, so that it is more general than those given in the literature. With the aid of such a recast version of Graf's addition formula, the coordinate transformation for Bessel functions between two polar coordinate systems arbitrarily located can be performed conveniently. Both the steady-state and transient alterations in the surface and subsurface motions have been evaluated and analysed. There is good consistency between the computed results obtained from the derived series solution and those from the BEM and existing literature. Spectral properties of ground-motion responses show that for grazing incidence the motions on the illuminated side of the canyon are amplified significantly when the dominant lengths of topographic features are about 1.5 times longer than the incident wavelengths. Topographic geometry strongly influences the surface and subsurface responses with dramatic changes at higher frequencies.

## Funding statement

K.-H.C. and J.-H.W. are sincerely grateful for the financial support granted by the Ministry of Science and Technology (former National Science Council), Taiwan, Republic of China (project nos. NSC 101-2811-M-001-167 and NSC 101-2116-M-001-015, respectively).

## Acknowledgements

The authors are much obliged to the editor and two reviewers for their encouragement and the insightful and valuable suggestions and comments that have improved this paper.

## Appendix A

In (2.22) and (2.23), the associated functions used for brevity are listed as follows:

- Received March 17, 2014.
- Accepted December 5, 2014.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.