2.5-D Time-Domain Finite-Difference Modelling of Teleseismic Body Waves

Analysing the teleseismic waveforms, we often calculate the synthetic seismograms. For complex structures such as subduction zones, however, it may be difficult to calculate accurate synthetic waveforms because of strong lateral heterogeneity. The laterally varying features such as steep sea-bottom topography and thick sedimentary layers can have a large effect even on long-period teleseismic body waveforms. For example, we can expect the large-amplitude later phases as the result of the structural effect, which cannot be predicted by the flat-layered model structure usually assumed in teleseismic-waveform analysis.


Introduction
Teleseismic-waveform analysis is one of the most effective approaches for study of the earthquake source process.It is also useful for investigation of the subsurface structures since the teleseismic seismograms have much information on the structures beneath the stations as well as on those near source and around the paths between the source and the stations.
Analysing the teleseismic waveforms, we often calculate the synthetic seismograms.For complex structures such as subduction zones, however, it may be difficult to calculate accurate synthetic waveforms because of strong lateral heterogeneity.The laterally varying features such as steep sea-bottom topography and thick sedimentary layers can have a large effect even on long-period teleseismic body waveforms.For example, we can expect the large-amplitude later phases as the result of the structural effect, which cannot be predicted by the flat-layered model structure usually assumed in teleseismic-waveform analysis.
A full treatment of such effect requires a three-dimensional (3-D) model of the structure and a 3-D calculation for the wavefield, which requires 3-D numerical techniques such as the 3-D finite-difference or finite-element method.Recent advances in high performance computers have already brought full 3-D elastic modelling for seismic wave propagation within reach.Even a single CPU computer could now be used for full 3-D numerical simulations by exploitation of a single or multi-GPU (Graphics Processing Units) computing (e.g., Okamoto et al., 2010).However full 3-D modelling of large-scale seismic wave propagation is still computationally expensive due to its requirements for large memory and a large number of fast processors, and would be too costly even on parallel hardwares for solutions of large-sized problems in routine-like real data analyses because of many case computations.Nevertheless, in order to provide a quantitative analysis of real seismic records from complex regions such as subduction zones, we need to be able to calculate the 3-D wavefields.
An economical approach to modelling of seismic wave propagation which includes many important aspects of the propagation process is to examine the three-dimensional response of a model where the material parameters vary two-dimensionally.Such a configuration in which a 3-D field is calculated for a 2-D medium is sometimes called two-and-a-half-dimensional (2.5-D) problem (e.g., Eskola & Hongisto, 1981).As a compromise between realism and computational efficiency, 2.5-D methods for calculating 3-D wavefields in 2-D varying structures have been developed.Bleistein (1986) developed the ray-theoretical implications of

2.5-D elastodynamic equation for a plane-wave incidence
We first use the physical properties of the wavefield to derive a 2.5-D elastodynamic equation in the time domain for the situation of an incident plane wave.Throughout this article we employ a Cartesian coordinate system [x, y, z], where the x and y are the horizontal coordinates and z is the vertical one.
For an isotropic linear elastic medium, the source-free 3-D elastodynamic equation in the time domain is given by , y, z, t) are the displacements at a point (x, y, z) at time t, and the stress components are τ rs = τ rs (x, y, z, t), (r, s = x, y, z).W e have used a contracted notation for derivatives ∂ tt ≡ ∂ 2 /∂t 2 , and ∂ r ≡ ∂/∂r, (r = x, y, z).The stress and displacement components are related by the 3-D Hooke's law through the Lamé constants λ = λ(x, y, z) and µ = µ(x, y, z) as follows: Numerical modelling schemes such as the finite-difference and pseudospectral methods can compute directly discretised versions of equations ( 1) and ( 2), where the bounded computational domains are usually represented by grids.
We now derive a 2.5-D equation of motion for a 3-D wavefield in a 2-D medium which is constant with respect to one coordinate and varies with the other two coordinates.We will assume the medium is constant in the y-direction throughout the rest of this article, so that the material properties take the form Furthermore we assume the medium includes a homogeneous half-space underlying the 2-D heterogeneous region whose top may be bounded by a free surface.
Consider an upgoing plane wave with horizontal slowness [p x , p y ], which passes a point [x 0 , y 0 , z 0 ] in the homogeneous half-space at a time t = t 0 .The 3-D wavefield at arbitrary time and position in the medium including the 2-D heterogeneous region has the characteristic of repeating itself with a certain time delay for different observers along the medium-constant axis (i.e., y-axis).For instance, the wavefield in the vertical plane y = y 0 at the time t = t 0 must be identical to that in the vertical plane y = 0 at the time t = t 0 − p y y 0 .This means u(x, y, x, t)=u(x,0,z, t − p y y), τ rs (x, y, x, t)=τ rs (x,0,z, t − p y y), (r, s = x, y, z).(4) If the structure is invariant in both of the horizontal (x-and y-) directions so that the material properties depends only on the vertical (z-) direction (i.e., 1-D heterogeneous medium), equation ( 4) might reduce to u(x, y, x, t)=u(0, 0, z, t − p x x − p y y), τ rs (x, y, x, t)=τ rs (0, 0, z, t− p x x− p y y), (r, s = x, y, z).
(5) Note that this is 'Snell's law' for plane-wave propagation in a 1-D heterogeneous medium.Equation ( 4) is thus a 2.5-D version of the 'Snell's law' which is also mentioned below.Equation ( 5) could be used for modelling three-component seismic plane waves in vertically heterogeneous media in the time domain (JafarGandomi & Takenaka, 2007;2010;Tanaka & Takenaka, 2005).
Let us be back to the 2.5-D problem.From relations (4), the derivatives of the displacement and the stress with respect to y can be expressed as ∂ y u(x, y, x, t)=−p y ∂ t u(x, y, z, t), ∂ y τ rs (x, y, x, t)=−p y ∂ t τ rs (x, y, z, t), where ∂ t ≡ ∂/∂t.The equivalent relations for the stress in ( 4) and ( 6) can be derived directly from those for the displacement in ( 4) and ( 6) through equations ( 2) and (3).
Substituting ( 6) into ( 1) and ( 2) we obtain the equation of motion and the stress representations This set of equations represents the 2.5-D elastodynamic response of a medium in the absence of source.Note that all the variables in this set of equations are real-valued.When we solve equations ( 7) and ( 8), we can set y = y 0 , so that these equations are reduced to 2-D ones.Once equations ( 7) and ( 8) have been solved for y = y 0 , we can deduce the wavefield at any y from that at y = y 0 by shifting the time origin by p y (y − y 0 ) (see equation ( 4)).
We next give an alternative derivation of these time-domain 2.5-D equations ( 7) and ( 8), from the 2.5-D equations in the frequency-wavenumber domain, and recover the characteristic of the 2.5-D wavefield, equation (4), in the process of deriving these equations.Fourier-transforming the 3-D equations ( 1) and ( 2) with respect to t and y, we obtain the following source-free 2.5-D elastodynamic equation in the frequency-wavenumber domain: −ω where we have used the y-invariance of the medium, i.e. equation (3), and have used the for the transform to the frequency-wavenumber domain.For a fixed value of the wavenumber k y , equations ( 9) and (10) depend on only two space coordinates, i.e. x and z.For each value of k y , these equations can therefore be solved as independent 2-D equations.The invariance of the medium in the y direction means that there is no coupling between different k y components.Whereas for full 3-D problems there would be coupling between different k y wavenumber components.For 2.5-D problems with an incident plane wave, we need to consider only one k y for each ω, which is that of the incident plane wave.
The inverse transform of the double Fourier transform (11) is Changing the order of the integrations, and inserting the following relation between the wavenumber k y and the slowness p y : equation ( 12) leads to where ĝ in the time-slowness domain has been defined as We then find and 1 2π For an incident plane wave in a 2.5-D situation, the horizontal wavenumber k y of all wavefields is constant for each ω, and equal to that of the incident plane wave.Further from (13) we require p y to be invariant and equal to the y-component of the slowness of the incident plane wave, which represents 'Snell's law' for 2.5-D problems as mentioned above.Thus, when the slowness of the incident wave is p y0 ,û(x, p y , z, t − p y y) can be represented as where δ(x) is the Dirac delta function.Applying the inverse transform ( 12) to the displacement in the frequency-wavenumber domain ũ(x, k y , z, ω), and using equations ( 14), ( 16) and ( 18), we obtain u(x, y, z, t)= û(x, p y0 , z, t − p y0 y). Then, We can recover ( 4) from ( 19) by appropriate substitutions: setting y to 0 we obtain an equation at time t and then making the particular choice t − p y0 y gives In a similar way, we can obtain the corresponding equation for the stress.Applying the partial Fourier inversions ( 16) and ( 17) to ( 9) and ( 10), and using equations ( 19) and ( 20), we recover the earlier forms ( 7) and ( 8).
In equations ( 7) and ( 8) the time derivatives appear on both sides of these equations, which may be inconvenient for direct discretisation with the finite-difference method.Here instead of direct use of equations ( 7) and ( 8), we employ the 2.5-D equation in terms of velocity-stress that is well suited to the use of the staggered-grid finite-difference technique.After some manipulation of ( 7) and ( 8) (Takenaka & Kennett, 1996b), we get the following velocity-stress formulation of the 2.5-D elastodynamic equation: y, z, t) is the particle velocity, and , with P-wave velocity α and S-wave velocity β.When we solve (22), we can set y = 0 as well as the case of ( 7) and (8) .

Finite-difference scheme
We use a staggered-grid finite-difference scheme (e.g., Hayashida et al., 1999;Levander, 1988;Virieux, 1986), which is stable for any values of Poisson's ratio, making it ideal for modelling marine problems.The finite-difference grid is staggered in time and two-dimensional space (x-z plane) as shown in Fig. 1.The y-components of particle velocity v locates at the same grid points as the normal stresses both in time and space.We should note the y coordinate is not discretised but continuous.Letting x = i∆x or x =( i ± 1/2)∆x, z = j∆z or z =(j ± 1/2)∆z, and t = l∆t or t =( l ± 1/2)∆t ; ∆x and ∆z are the grid spacings in the xand z-direction, respectively, and ∆t is the time step, and using Levander's notation (Levander, 1988), the difference equations for (22) are, for example, where D + t is forward difference operator in time, and D ± x and D ± z are the forward-or backward-difference operators in space, with sign chosen to center the difference operator about the quantity being updated.For example, in case of a second-order accurate in time and fourth-order accurate in space scheme we used, where The fourth-order spatial finite-difference scheme usually needs more than five grid points per wavelength (Levander, 1988).The finite-difference region (computational domain) is divided into two parts: a model zone for the upper part and an initial wave zone for the lower part.The model zone fully includes the target structural model and may be heterogeneous.The initial wave zone locates under the model zone and should be homogeneous.It is needed for the incident wave, and an upgoing plane wave is input as the initial condition within it.The initial wave zone should have no velocity contrast at the interface contacting the model zone to prevent artificial reflections and conversions.The size of the computational domain is set sufficiently large so that artificial noises, such as noises from the both ends of the input plane wave, can be ignored.As Okamoto (1994) we here select the time step as 62 % of the maximum allowed time step by the stability condition for the 2-D P − SV finite-difference scheme of a second-order accurate in time and fourth-order accurate in space (Levander, 1988).

Teleseismic waveform synthesis for source inversion
Seismic displacement due to a point source of earthquake may be expressed as where M pq (t) is moment tensor with time varying components, operation * is convolution with respect to time t, and G n;p (x, t; ξ, τ) is displacement Green's tensor representing the nth component of elastic displacement at a receiver position x and time t caused by a unit point force in the p-direction at a source position ξ and time τ (e.g., Aki & Richards, 2002).We have used the convention of summation over repeated suffices.In source inversions derivative of the Green Ąfs tensor ∂ q G n;p (x, t; ξ,0) is empirically called just "Green's function".We here follow this custom.
Applying the spatial reciprocity: equation ( 26) reduces to where E pq;n (ξ, t; x,0) is the strain tensor at the source position corresponding to G p;n (ξ, t; x,0).This equation shows as follows.The reciprocity of the elastodynamic theory is exploited 312 Seismic Waves, Research and Analysis www.intechopen.comto calculate the Green's functions.The displacement at receiver due to the moment tensor at source can be calculated by solving the reciprocal problem; the reciprocal relation is applied to the response displacement at the source position due to a unit body force acting at the receiver position.In the teleseismic problem, the virtual force is applied at infinity (i.e., station location).Since we are only concerned with the teleseismic body waves, the wavefield due to this virtual force can be approximated by a plane wave.The response of the near-source structure to an incident plane wave is then calculated and converted to the far-field displacement.This approach using the reciprocity for teleseisimic body wave synthesis was employed by Bouchon (1976) where a simple method was presented, which combines the reciprocity theorem and the flat layer theory (Thomson-Haskell matrix formulation: Haskell, 1953;Thomson, 1950) to yield teleseismic body wave radiation from seismic sources embedded in the horizontally layered crustal models.
For the conversion to the far-field displacement, based on the assumption of ray propagation in the mantle, the response of the near-source structure mentioned above is multiplied by the geometrical spreading effect within the mantle upon the amplitude of the body wave and convolution with the source time function, with atenuation operator for the path in the mantle, with the response at the surface of the receiver crust for the incident impulsive teleseismic body wave, and with the instrumental response, where the actually employed response of the receiver crust often accounts for only the free surface effect at the receiver or is calculated from a very simple crustal model (this issue may be related to the subject of the next subsection).This method has been extensively used to calculate teleseismic waveforms for studies on earthquake source processes (e.g., Miyamura & Sasatani, 1986).Even now it is one of the most popular methods for calculation of teleseismic Green's functions in routine-like source inversions (e.g., Kikuchi & Kanamori, 1991;Nakamura et al., 2009).
The reciprocal formulation has the following advantages: for a single station the Green's functions for many point sources at different positions can be obtained simultaneously in a single numerical calculation, which facilitates the waveform analysis to find the best position for the earthquake source location without repeating the time-consuming numerical calculations.Although Bouchon (1976) actually treated calculation for horizontally layered structures, he mentioned in the paper that the reciprocal approach is also applicable to the case of a source in a layered structure with irregular interfaces.Wiens (1987;1989) and Yoshida (1992) applied this approach to planar-dipping structures including the sea floor by using a ray-theoretical technique developed by Langston (1977).Furthermore Okamoto & Miyatake (1989) and Okamoto (1993) extended it to arbitrary 2-D heterogeneous structures by using the finite-difference method in time domain.However, their calculation is based on the 2-D elastodynamic equation, and is limited to two dimenions, so that the available stations are restricted to those located in the direction perpendicular to the medium-constant axis (y-axis in the previous sections).This restriction in the azimuthal coverage makes it difficult to examine the source process in detail.
In order to get teleseismic synthetics for arbitrary azimuth Okamoto (1994) proposed a method for calculating the 2.5-D telseismic body waves, which solves the time-domain version of 3-D elastodynamic equation in the mixed coordinate (x and z)-wavenumber (k y ) domain ( 9) and (10) for each of a number of discretised wavenumber k y with the finite-difference time-domain scheme and perform an inverse Fourier transform over wavenumber k y (i.e.wavenumber summation) to obtain the synthetic seismograms in the spatial domain.His method requires the computation time more than hundreds times of the corresponding 2-D computations.et al. (1997) presented a method without wavenumber summation so that 2.5-D teleseismic synthetics requires only computation time similar to the corresponding 2-D ones.This method uses the 2.5-D elastodynamic equation for a plane-wave incidence ( 22) and its finite-difference time-domain scheme proposed by Takenaka & Okamoto (1997) which was described in the previous section.It has been employed for calculating teleseismic Green's functions in several source inversion studies (e.g., Okamoto & Takenaka, 2009a;b) and for modelling the teleseismic waveforms including a near-source scattering inside a subducted plate (Kaneshima et al., 2007).We here show two results for source inversion among them as examples.
Figure 2 2(a)) are allowed to vary with respect to the trench-perpendicular axis, while they are assumed to be invariant with respect to the trench-parallel axis.This model was constructed from the results of seismic surveys in the nearby area (Kopp et al., 2002) and global reference models (Dziewonski & Anderson, 1981;Kennett & Engdahl, 1991;Laske et al., 2001).In Fig. 2(a) the point sources for Green's function computations are located along the dip of the main-shock fault plane.The along-dip interval of source points is 8.0 km for the section from S1 to S8 and 8.1 km for the section from S8 to S15.The rigidity for sources S1-S7 is 16.3 GPa and for sources S8-S15, 38.6 GPa.The best double couple of the Global CMT solution (http://www.globalcmt.org)shown in Fig. 2(b) is employed for each Green's function computation.Following the standard 1-D teleseismic wave computations, mantle attenuation is incorporated by choosing t * = 1.0 for P-waves and 4.0 for SH-waves, while anelastic attenuation is not included in the finite-difference computations that evaluate near-source response.The 1-D Green's functions were computed by the method of Kroeger & Geller (1983).The comparison between waveforms of the 2.5-D and flat-layered Green's functions (Fig. 2(c)) clearly illustrates the large effect of the heterogeneous structure on the body waves.The waveforms of 2.5-D Green's functions have prolonged, large amplitude later phases.They appear irrespective of station azimuth, and are not reproduced by 1-D model.At the oceanic trench regions large effect of laterally heterogeneous structure is expected to appear on the teleseismic body waveforms: thick water layer, dipping ocean bottom, and thick sediments near the source distort ray paths to teleseismic stations and often cause large later phases on the teleseismic body waveforms.This effect must be evaluated carefully before a detailed source process analysis is applied to real earthquake records.
Okamoto & Takenaka (2009b) studied strong effect of near-source structure on teleseismic body waveforms from two well-recorded aftershocks of the 2006 Java tsunami earthquake.
Figure 3 shows the results of one of the two events: M W 6.1, 2006/07/17 15:45:59.8,9.420 • S 108.319 • E. They applied a "waveform relocation technique" which combines a waveform inversion of source parameters with a grid search procedure to correct possible systematic bias in hypocentral parameters.In the waveform inversion 2.5-D teleseismic Green's functions are used.The grid spacing for grid search is 2 km horizontally and 1 km vertically.In Fig. source are close to those of the Global CMT.The residual contour distribution in Fig. 3(a) and the RMS error plots in Fig. 3(c) indicate that the source location could be well constrained both vertically and horizontally.Use of the 2.5-D modelling makes it possible to obtain improved source parameters at the trench regions where only teleseismic data are available.

Modelling for receiver function analysis
Receiver function analysis is one of the effective and popular methods for study of the crust and upper mantle structures using teleseismic waveform data (e.g., Ammon, 1991;Cassidy, 1992;Dugda et al., 2005;Farra & Vinnik, 2000;Kanao & Shibutani, 2011;Langston, 1979;Owens et al., 1988;Saita et al., 2002;Suetsugu et al., 2004;Zhu & Kanamori, 2000).In the analysis it is often necessary to calculate synthetic waveforms for the structure models.For this purpose horizontally layered structure models have been assumed, because the response  of such a simple structure model to teleseismic P wave (plane wave) can be calculated easily and accurately by a semi-analytical method such as the Thomson-Haskell matrix method (Haskell, 1953;Thomson, 1950).However, for structures with strong lateral heterogeneity such as subdunction zones, it is often difficult to consider horizontally layered media for modelling the teleseismic body waves that propagate through such complex structures.Full 3-D modelling of seismic wave propagation is still computationally intensive.Takenaka & Okamoto (1997) used the 2.5-D finite-difference method described in the previous sections to simulate teleseismic seismograms at ocean-bottom stations for assessing the effect of sea-bottom topography.Ando et al. (2003) applied this approach to a profile across a realistic model of subduction zone structure to simulate the effect of a subducting slab on the receiver function waveforms observed at subduction zone, and Takenaka et al. (2004) then clearly demonstrated the azimuthal dependence of the slab-converted phases in the receiver functions.We here illustrate their results.Fig. 4. Pand S-wave velocity profiles of a subduction zone model.The slab consists of two layers: the upper layer of 7 km thickness corresponding to the oceanic crust has velocities of 6 % lower relative to the mantle of the ak135 model (Kennett et al., 1995), while the lower layer corresponding to the oceanic slab mantle has velocities of 5 % higher relative to the ak135 mantle.Figure 4 shows a cross-section of a subduction zone model which is an analogue to that of Tohoku, Japan.The shown area is the target where all the seismic phases for the receiver functions are modelled.The actual computational domain is set sufficiently larger so that artificial reflection noises from the bottom and both sides of the domain do not contaminate the synthetic seismograms.The synthetic seismograms of teleseismic P wave at the ground surface between A and B in Fig. 4 were calculated for events of epicentral distance 80 • with backazimuths of 0 • to 180 • in the interval of 22.5 • .The definition of the backazimuth (φ) is indicated in Fig. 5.The epicentral distance 80 • corresponds to the P-wave incident angle of around 17 • at the surface.Three-components of the synthetic seismograms for the teleseismic events of backazimuths of 0 • and 90 • are shown in Fig. 6.The signal of the incident wave (source wavelet) was assumed to be a simple Gaussian pulse.The radial and transverse receiver functions for each backazimuth event can be calculated from the synthetic seismograms (Fig. 7).

Conclusion
Two-and-half-dimensional approach in the time domain has been considered for calculating three-dimensional teleseismic body wave propagation in two-dimensionally varying medium.It is an economical approach for calculating 3-D wavefields in real problems, and requires storage only slightly larger than, and computation time only slightly longer than those of the corresponding 2-D calculation.A finite-difference scheme solving the 2.5-D elastodynamic equation for 3-D seismic response of a 2-D structure model due to an obliquely incident plane wave has been described.The modelling of such seismic wavefields for a 2.5-D situation with an incident plane wave is of considerable practical importance.For instance, this approach can be applied to modelling of teleseismic body waves observed on complex crustal structures or radiated from shallow earthquakes occurring in subduction zones, where the laterally heterogeneous media can have large effects on the waveforms.We showed some numerical examples which include modelling of teleseismic body waves for the earthquake source analysis and the receiver function analysis.The method described here is very efficient, so that we expect it could be used for waveform inversions for routine-like source retrievals 322 Seismic Waves, Research and Analysis www.intechopen.comand subsurface structure reconstructions from teleseismic seismograms in the near future, which need many forward modelling computations.

310
Seismic Waves, Research and Analysis www.intechopen.com

Fig. 2 .
Fig. 2. (a) 2-D model of the Java trench.P-wave velocity is shown in colour scale (S-wave velocity is zero in the ocean).Green circles denote point source positions along the dip for Green's function computations.(b) Global CMT solution of the main shock and the teleseismic station coverage.(c) Examples of synthetic P-waveforms (Green's functions) for station MA2.2.5D denotes those computed for the 2-D model of the near-source structure, and 1D denotes those for the flat-layered model.Attached indexes (S2-S14) denote the source positions.Numbers attached to the 1-D waveforms denote cross-correlation coefficients between 1-D and corresponding 2.5-D Green's functions for a dulation from 0 to 90 s.The 1-D model consists of a standard crust(Kennett & Engdahl, 1991) additionally overlain by a 3-km-thick ocean and 2-km-thick sediment.(d) Same as (c) but for station MBAR.(Reproduced fromOkamoto & Takenaka, 2009a).

Fig
Fig. 3. (a) Open star indicates the best point source position of the 17-July-2006 event (M W 6.1).Locations of Global CMT (triangle) and PDE (diamond) are also projected.The contour lines denote the residual distribution of the grid search relocation by the waveform inversion.The contour interval is 0.02.(b) Observed (top) and 2.5-D synthetic (bottom) waveforms.Attached number denotes the maximum amplitude of the observed waveform in µm.Also plotted are the source time function (STF) and the focal mechanism.A time window of 70 s after the onset (indicated by vertical lines) is used for the inversion.The estimated moment tensor components in unit of 10 17 Nm are: M rr = −1.81,M θθ = 4.63, M φφ = −2.81,M rθ = 11.1,M rφ = −2.04,M θφ = 0.67, which yield a scalar moment of 1.19 × 10 18 Nm (M W 6.0).(c) RMS error in travel time analysis plotted versus the distance with respect to trench-parallel axis (positive toward N116 • E with the origin placed on the cross section through the PDE epicenter).Most of the travel times listed in USGS NEIC Monthly Earthquake Data Report are used.(Reproduced from "Effect of near-source trench structure on teleseismic body waveforms: an application of a 2.5D FDM to the Java trench" by T. Okamoto & H. Takenaka, in Advances in Geosciences, Vol. 13 (Solid Earth), Ed.Kenji Satake, Copyright (C) 2009 by World Scientific Publishing.)

Fig. 5 .
Fig. 5. Geometrical definitions.θ is the incident angle of a plane P wave, and φ is the backazimuth.

Figure 8 (
Figure 8(a) displays circular plots of the radial receiver functions over backazimuth of events for three stations.The receiver functions for backazimuths of more than 180 • can be obtained from the synthetic seismograms for backazimuths of less than 180 • through symmetrical

Figure 8 (
b) shows linear plots of the radial and transverse receiver functions for the left station among the three stations shown in Fig.8(a).In both components of the receiver functions the slab Ps phases generated at dipping interfaces are clearly seen as convex arrival patterns with the latest arrival around 17 s at backazimuth of 180 • .The slab-converted phases exhibit the amplitude and arrival-time variations as a function of the backazimuth: the arrival is the latest for the backazimuth (φ) which is equal to the dip direction of the slab (φ 0 = 180 • ); and the amplitude variation shows a pattern with the backazimuth like cos[(φ − φ 0 )/2] for the radial receiver function, and cos(φ − φ 0 ) for the transverse one, respectively.