# Three-dimensional pre-stack depth migration of receiver functions with the fast marching method: a Kirchhoff approach

Cheng Cheng, Thomas Bodin, Richard M. Allen

► To cite this version:

Cheng Cheng, Thomas Bodin, Richard M. Allen. Three-dimensional pre-stack depth migration of receiver functions with the fast marching method: a Kirchhoff approach. Geophysical Journal International, 2016, 205 (2), pp.819-829. 10.1093/gji/ggw062. hal-02331135

HAL Id: hal-02331135

https://univ-lyon1.hal.science/hal-02331135

Submitted on 10 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Three-dimensional pre-stack depth migration of receiver functions with the fast marching method: a Kirchhoff approach

Cheng Cheng,<fnref n="1" /> Thomas Bodin<fnref n="1" />,<fnref n="2" /> and Richard M. Allen<fnref n="1" />

Accepted 2016 February 9. Received 2016 February 8; in original form 2015 May 15

## SUMMARY

We present a novel 3-D pre-stack Kirchhoff depth migration (PKDM) method for teleseismic receiver functions. The proposed algorithm considers the effects of diffraction, scattering and traveltime alteration caused by 3-D volumetric heterogeneities. It is therefore particularly useful for imaging complex 3-D structures such as dipping discontinuities, which is hard to accomplish with traditional methods. The scheme is based on the acoustic wave migration principle, where at each time step of the receiver function, the energy is migrated back to the ensemble of potential conversion points in the image, given a smooth 3-D reference model. Traveltimes for *P* and *S* waves are computed with an efficient eikonal solver, the fast marching method. We also consider elastic scattering patterns, where the amplitude of converted *S* waves depends on the angle between the incident *P* wave and the scattered *S* wave. Synthetic experiments demonstrate the validity of the method for a variety of dipping angle discontinuities. Comparison with the widely used common conversion point (CCP) stacking method reveals that our migration shows considerable improvement. For example, the effect of multiple reflections that usually produce apparent discontinuities is avoided. The proposed approach is practical, computationally efficient, and is therefore a potentially powerful alternative to standard CCP methods for imaging large-scale continental structure under dense networks.

**Key words:** Body waves; Computational seismology; Wave scattering and diffraction.

## 1 INTRODUCTION

Teleseismic receiver functions (e.g. Phinney 1964; Burdick & Langston 1977; Vinnik 1977; Langston 1979) have become a standard method to map seismic discontinuities in the crust and upper mantle beneath a broad-band seismometer. The standard process is to deconvolve the vertical component from the radial component in the coda of teleseismic *P* wave seismograms. This removes source and instrument effects, while preserving *P*-to-*S* converted phases generated at discontinuities beneath the receiver, which contain a significant amount of information on seismic structure.

There are various ways to analyse and interpret receiver functions. For example, one may invert receiver functions for 1-D profiles of velocity structure individually (Kind et al. 1995; Sandvol et al. 1998) or jointly with surface wave dispersion data (Julià et al. 2000; Shen et al. 2013; Bodin et al. 2014). Alternatively, Zhu & Kanamori (2000) and Chevrot & Van der Hilst (2000) introduced a grid search method over Moho depth ($H$) and $Vp/Vs$ ($\kappa$), the so-called $H-\kappa$ stacking technique, transforming time-domain receiver functions into the depth-$Vp/Vs$ domain. This is a widely used method in seismology to estimate crustal properties (e.g. Julià et al. 2003; Lombardi et al. 2008; Spieker et al. 2014, etc.).

Receiver functions observed at different stations of an array of seismometers may be combined into a 2-D or 3-D volume and interpreted structurally. This process often utilizes common conversion point (CCP) stacking (Dueker & Sheehan 1997), which stacks the amplitude of receiver functions from several station-event pairs sampling the same subvolume. Both the signal-to-noise ratio and the spatial resolution can be significantly improved by optimum CCP binning, moveout correction and subsequent stacking thus providing important constraints on crustal and upper-mantle structure. For example, the USArray seismic network has recently facilitated a number of CCP studies under the US continent (Abt et al. 2010; Schmandt et al. 2012; Tauzin et al. 2013; Lekić & Fischer 2014). CCP stacking has also been used in other regions where dense seismic networks are available, such as in Tibet, China (Caldwell et al. 2013; Shi et al. 2015); in the North China Craton (Ai et al. 2007; Zheng et al. 2014; Liu et al. 2015) and in Europe (Steckler et al. 2008).

However, early CCP stacking relies on the assumption of a 1-D Earth where discontinuities are horizontal. Despite the great success of CCP imaging, numerous studies document that complex 3-D structures as well as multiple reflections may result in significant artefacts in the final image. Both Abers (1998) and Morozov (2004) demonstrate the coherence and high amplitude of signals generated

<footnote n="1">¹Department of Earth and Planetary Science, University of California, Berkeley, CA 94720, USA. E-mail: cchengster@berkeley.edu</footnote>

<footnote n="2">²Laboratoire de Géologie de Lyon: Terre, Planètes et Environnement, CNRS, Université de Lyon 1, École Normale Supérieure de Lyon, F-69622 Villeurbanne, France</footnote>
by dipping, out-of-plane, structures. Later, a technique for correct-
ing for 3-D velocity perturbations along the 1-D ray path using
the linear tomography assumption with the traveltime corrections
calculated from a hybrid 3-D velocity model (Levander & Miller
2012) is introduced into CCP. But limited research has been done to
demonstrate the improvement over the earlier one in the presence of
heterogeneity. Also even if the smoothed background 3-D velocity
perturbation is ‘corrected’, CCP processing does not restore dip-
ping events or collapse diffractions because of intrinsic limitation
of CCP method, such as the horizontal layer assumption of the pre-
mapping stack (Rondenay 2009; Levander & Miller 2012). Overall,
CCP stacking presents three well-known issues: (1) it gives poor
results in the presence of lateral volumetric velocity variations, (2)
it fails at imaging dipping discontinuities and (3) final images are
strongly polluted by the signal of multiply reverberated waves.

Alternatively, more complex migration techniques can be used to image discontinuities (Claerbout 1985). Migration schemes are extensively used in exploration geophysics, and are an active field of research with a wide number of algorithms available, such as Kirchhoff migration (Schneider 1977), or wave-equation migration (Biondi & Sava 1999; Sava & Biondi 2004a,b). Compared with CCP, they require more computation and memory as each observed seismic phase is relocated on a large number of different gridpoints of the image, representing the location of a potential scatterer. The main advantage of migration techniques over CCP 1-D mapping is that they make fewer assumptions about the geometry of the subsurface structure. Instead of assuming that the Earth is made of planar, horizontal layers, migration algorithms treat scattering from 2-D and 3-D structures.

Pre-stack Kirchhoff migration of seismic records is a well-established technique in industry to image the shallow sedimentary structure (e.g. Schneider 1977; Wiggins 1984; Claerbout 1985). It was later introduced to seismology for imaging larger scale structures as local networks became dense enough (e.g. Rondenay et al. 2001; Levander et al. 2005). Teleseismic Kirchhoff migration is an imaging scheme applied to the scattered wavefield associated to incident body waves recorded at an array of stations. The data are weighted and stacked along diffraction hyperbolae for every potential scattering point in a regular grid defining the model space. Effectively, one sums over all the available data that have a traveltime consistent with a scatterer at the target point. All the teleseismic Kirchhoff migration techniques that have been developed to date consider only first-order interactions between the incident wavefield and subsurface perturbations (i.e. they make use of the single scattering, or Born approximation) (Rondenay 2009). Depending on how the scattering coefficient is treated, teleseismic migration schemes can be separated into two groups. One group is based on acoustic scattering, it treats the incoming *P* wave and the scattered *S* wave both as scalar, and involves stacking of singly scattered wavefields along diffraction hyperbolae to recover relative scattering intensity/potential at individual points through a 2-D or 3-D model space (Revenaugh 1995; Ryberg & Weber 2000; Sheehan et al. 2000; Levander et al. 2005). Another group of Kirchhoff migration is based on elastic scattering, which treats the incoming *P* wave as scalar and scattered *S* wave as a vector and combines stacking with inversion/backprojection operators (Bostock & Rondenay 1999; Bostock et al. 2001; Poppeliers & Pavlis 2003a,b).

Recently, other migration techniques have been borrowed from
industry and applied to teleseismic receiver functions, such as one-
way wave-equation migration (Chen et al. 2005), which is a 2-D
scheme that still relies on 1-D horizontal layer assumption for the
moveout correction; teleseismic shot profile migration (Shragge

et al. 2006), which is also a 2-D wavefield extrapolation scheme. Another example is reverse time migration (Shang et al. 2013), which highly depends on computation capability and image scale. These migration methods are based on numerical solver for the wave equation, where the full wavefield is accurately modeled and finite-frequency effects are accounted for.

Finally, full 3-D waveform inversion (Frederiksen & Revenaugh 2004) of scattered waves can be used to recover either scattering potential or estimates of localized material property perturbations relative to the background model. Both acoustic and elastic methods can be highly computationally expensive, especially when it comes to 3-D wavefield extrapolation using finite-difference schemes, or when a large amount of receiver functions are incorporated. For a complete review of available methods, we refer the reader to the review papers by Bostock (2007) and Rondenay (2009).

In this work, we present an efficient 3-D pre-stack Kirchhoff
depth migration method for teleseismic receiver functions. It highly
resembles the acoustic wave migration scheme, as it can account
for 3-D dipping structures and strong volumetric heterogeneities.
However, this migration scheme is easy to implement and as fast and
efficient as CCP depth mapping. It leads to an efficient algorithm for
performing migration in 3-D, which will allow for the possibility
of using it for studies where CCP stacking was previously the only
plausible approach. In Section 2, we describe the theoretical aspects
of our method. In Section 3, we illustrate its advantages with some
synthetic tests and compare it with the traditional CCP method.

## 2 METHODOLOGY

As summarized in Rondenay (2009), the general framework for
Kirchhoff teleseismic wave migration is:

$$
f(r) = \sum_{j=1}^{N} \sum_{r_0=1}^{M} \vartheta(r, r_0, p_j^0) \Delta \tilde{u} [r_0, p_j^0, t = T(r, r_0, p_j^0)] \quad (1)
$$

where *f(r)* is the scattering potential at gridpoint *r* and the stack
is performed on *N* incident waves, measured at *M* stations. *Δtilde u* is
the amplitude of the scattered wavefield measured at station *r*₀ at
time *t* = *T(r*, *r*₀, *p*<sub>j</sub><sup>0</sup>) predicted using a smooth model while *p*<sub>j</sub><sup>0</sup>
is the slowness term for the particular teleseismic event. The sum-
mation term over *r*₀ represents a stack over different stations of
the weighted scattered wavefield. *ϑ(r*, *r*₀, *p*<sub>j</sub><sup>0</sup>) is a weighting term
that accounts for amplitude effects, such as geometrical spreading,
that is, the amplitude of a converted phase observed at the surface
decays as the conversion point gets deeper. Fig. 1 depicts the princi-
ples of Kirchhoff migration for receiver functions. This is a general
approach and there are various ways to implement it. Different
methods apply different simplifications to the weights or filters.

The main advantage of pre-stack migration techniques over CCP
depth mapping is that they make much fewer assumptions about
the subsurface structure, and thus consider seismic scattering in
more general terms. However, they still rely on a smoothly varying
background velocity model, necessary to backproject the scattered
wavefield. In practice, we can get this absolute velocity model from
either global or regional tomographic images. 3-D volumetric large-
scale velocity structure is essential for proper image focusing.

## 2.1 Acoustic scattered-wave migration scheme

Following Levander et al. (2005), we implement an acous-
tic scattered-wave migration scheme. The receiver function
<figure><img src="image_3.png" /><figcaption>Figure 1. Schematic illustrations of how Kirchhoff migration works for a teleseismic receiver function. (a) Migrated parabola of one receiver function, constructed from a simple model made of a horizontal layer over a half-space. A vertical slice is shown. (b) Synthetic receiver function used for migration. (c) Ray paths for direct P wave and transmitted S wave in the same vertical plane as (a). The converted Ps phase is migrated to all possible conversion points that predict the observed traveltime t.</figcaption></figure>

time-domain migration equation writes:

$$ 
\begin{aligned}
f(r) &= S_{PS}(r) \\
&= \sum_{r_0=1}^{N} \sum_{r_e=1}^{M} \mathrm{RF}(r_e, r_0, t = (\tau_P + \tau_S - \tau_e)) \frac{A_S(r, r_0) \cos\theta(r_0)}{A_P(r, r_e) \beta(r_0)}
\end{aligned} 
\quad (2) $$

In which, r is the imaging gridpoint, r_e is the source position and RF is the receiver function data recorded at station r_0, which shall be migrated to depth. A_P(r, r_e) represents the geometric spreading of the P wave from the source to the imaging grid and A_S(r, r_0) the geometric spreading of the transmitted S wave from the imaging grid to the station. β(r_0) is the surface shear wave velocity and θ is the incidence angle. The time (t = (τ_P + τ_S - τ_e)) is the imaging condition for receiver function Kirchhoff migration, in which τ_P is the P traveltime from source to imaging grid, τ_S is the S traveltime from imaging grid to receiver and τ_e is the P traveltime from source to receivers (see Fig. 1).

In Levander et al.'s approach the traveltimes and amplitude of the P and S wavefields are calculated using an upwind finite-difference solver to account for heterogeneity in the imaging region. Due to heavy computational costs in 3-D, the finite-difference scheme used is typically a 2-D solver. So the whole formula should be considered as a pseudo-3-D equation. For example, in Levander et al. (2005), the correction for out of plane propagation with respect to a reference coordinate, y', using the ray parameter in the coordinate orthogonal to the image, P_y, is included in the phase term (i.e. the time t).

A clear advantage of acoustic scattered-wave migration is the simplicity of the processing and the ability to treat simultaneously a large number of input data to produce images of scattering potential. A limitation of the approach is the lack of formal relationship between the scattering potential and actual material property perturbations in the scattering expressions (except for the P-to-P scattering mode, which is a scalar function of P-wave velocity perturbation) (Bostock et al. 2001; Rondenay 2009). In particular, this limitation affects images generated by scattering modes involving

S waves (e.g. P-to-S), but can be addressed by solving the problem for elastic waves.

Our method can be viewed as a simplification on Levander et al.'s (2005) acoustic equation in order to consider fully 3-D structure, improved by incorporating elastic scattering patterns, which will be discussed in detail in the following sections. Overall, the main advantage of our method is its computational efficiency, which makes feasible applications at continental scales.

Instead of calculating the full wavefield at each gridpoint by using an expensive 2-D finite-difference solver, we make some simplifications on the formula in eq. (2), which allow us to include the full 3-D velocity structure. The term A_P(r, r_e) represents the amplitude of the wave propagating from the teleseismic source to the scattering gridpoint. It has relatively small variations across the model box, and hence can be simplified to a constant across the whole image. Another term A_S(r, r_0) represents the geometrical spreading of the scattered S wave, that is, the converted S wave traveling from the gridpoint to the receiver. Instead of calculating it using a finite-difference solver in a heterogeneous model, here we simplify it as a 1/d term in a 3-D case (Rondenay 2009), where d is the Euclidian distance between r and r_0. Assuming only geometrical spreading is like assuming a constant velocity model, and ignoring focusing and defocusing effects due to volumetric heterogeneities. This may be adequate with a smoothly varying velocity field. In this way, these simplifications leave us with only the traveltime calculation, which we implement with the fast marching method (FMM; Rawlinson & Cambridge 2004).

The simplified migration equation takes the form of:

$$ f(r) = S_{PS}(r) = \sum_{r_0=1}^{N} \sum_{r_e=1}^{M} \mathrm{RF}(r_e, r_0, t = (\tau_P + \tau_S - \tau_e)) \frac{1}{d} \times \mathcal{W} \quad (3) $$

where W represents some weighting coefficient which will be discussed in detail later.
## 2.2 Fast marching method

The FMM (Sethian 1996; Sethian & Popovici 1999) is a grid-based numerical scheme for tracking the evolution of monotonically advancing interfaces via finite-difference solution of the eikonal equation. The eikonal equation states that the magnitude of the travel-time gradient at any point along a wave front is equal to the inverse of the velocity at that point and is written as $|\nabla T = s(r)|$, where $s(r)$ is the local slowness. To date, FMM has been applied to a wide variety of problems including seismic wave propagation, photolithographic development, geodesics, deposition of sediments and medical imaging. In this work, we use FMM as a quick and efficient way to estimate traveltimes between source, receiver and gridpoint pairs.

The principal drawbacks of standard ray-tracing schemes are related to robustness, speed and ray selection. It can also be a time-consuming process, especially if a large number of sources and receivers are involved and the medium is 3-D, as in our migration case. FMM distinguishes itself by combining both unconditional stability and rapid computation, making it a truly practical scheme for velocity media of arbitrary complexity. Some synthetic tests (Rawlinson & Cambridge 2004) have shown that wave fronts can be accurately tracked with minimal computational effort, even in the presence of complex velocity fields and layer boundaries with high curvature. These features of FMM make our migration algorithm computationally effective, as it only requires traveltime computations for all source-receiver gridpoint pairs (see cartoon in Fig. 1). For a detailed discussion on FMM, we refer the readers to Rawlinson & Cambridge (2004).

## 2.3 Traveltime matrix and stacking

The FMM algorithm also allows us to compute traveltimes, and store them in a matrix for a later migration process. We compute three different types of traveltimes: (1) P wave traveltimes from sources to stations $\tau_e(r_e, r_0)$; (2) P wave traveltimes from sources to imaging gridpoints $\tau_p(r_e, r)$ and (3) S wave traveltimes from imaging gridpoints to stations $\tau_S(r, r_0)$. The last one is calculated using reciprocity principle and shooting rays from stations and propagating the wave front through the 3-D volume to be imaged. Three matrices are thus constructed in three steps:

(1) For each event, run FMM to propagate a wave front starting from a point source and store the traveltime between the source and each station $\tau_e(r_e, r_0)$. The wave front is tracked in a combination of a global model and local model for the 3-D volume to be imaged.

(2) For each event, run FMM to propagate a wave front starting from a point source, and store the time between the source and each gridpoint $\tau_p(r_e, r)$.

(3) For each station, run FMM to propagate a wave front starting from the station and store the traveltime between the station and each gridpoint $\tau_S(r, r_0)$.

Since the entire traveltime field is computed for each run of FMM, steps (1) and (2) can be simultaneously computed by only running FMM once for each source. In this way, the FMM scheme needs only be run $N+M$ times, where $N$ is the number of stations and $M$ the number of sources. Furthermore, since each run of FMM is independent, calculation can be easily parallelized (with linear speed-up), thus reducing computational costs. Here, the computers used were limited in memory and the entire traveltime could not be stored in memory. The imaging points and stations were separated

into different groups, and several runs of FMM were carried out for the same source.

Once traveltimes are computed and stored, the staking scheme is simple to implement, and only consists of two loops. It can be summarized by the following pseudo-code:

For each gridpoint $r$

For each observed receiver function associated with source $r_e$ and station $r_0$

Stack the amplitude of the receiver function at time

$$t = \tau_p + \tau_S - \tau_e$$

End second loop

End first loop

As shown in eq. (3), the stacked amplitudes are weighted with a number of terms that we detail below. Since each step within both loops can be done independently, the stacking procedure can also be conveniently parallelized (also with linear speed-up). For the synthetic example shown below (i.e. for ~650 000 gridpoints and 27 000 receiver functions), the stacking procedure took about 12 CPU hours.

## 2.4 Weighting terms

In this work, the data, that is backpropagated and stacked, are radial receiver functions, which are constructed by deconvolving the vertical component from the radial one. Scattered shear waves are observed both on vertical and radial components. Depending on their incidence angle, they will contribute more or less to the radial component. A maximum contribution is expected when the converted S wave arrives vertically under the receiver. That is, the amplitude of converted phases observed on the receiver functions depends on the location of the conversion point. We simply correct for this effect by applying two weighting coefficients in eq. (3). Fig. 2(a) shows $\theta_1$, which is the angle between the line connecting the imaging gridpoint and the station and vertical direction. Fig. 2(c) shows the migrated parabola after applying the $\cos(\theta_1)$ (see eq. 5 below). Compared with the migrated image shown in Fig. 1, there is a cosine distribution to the amplitude along the parabola.

Future work will include rotating each receiver function to an $SV/P$ receiver function. Given the ray parameter of the incoming P wave, the incidence angle of the P wave at the surface can be given from the P velocity at the surface with a free surface transform (Abt et al. 2010; Yu et al. 2013).

Another weighting coefficient is applied for the reason that the main energy in the radial component is from the great circle path. The energy backprojected at gridpoints that are not in the great circle plane are therefore downweighted by a $\cos(\theta_2)$ term (see eq. 5 below), where $\theta_2$ represent the angle between the line connecting imaging gridpoint to the station and the great circle in the horizontal map view, as indicated in Fig. 2(b).

## 2.5 Elastic wave migration and elastic scattering pattern

In the acoustic approach presented above, P and S waves are scalar fields, and only phase information (traveltime) is extracted from the data. This is because until now the weighting term in eq. (3) only comprises geometrical spreading, phase shifting and projection into the radial component, without accounting for amplitudes of elastic interactions. In this way, for an incident P wave impinging upon a scatterer, one assumes that all resulting scattered waves radiate energy uniformly in all directions. However, this is not valid for P-to-S conversions. For example, as illustrated in Fig. 4, the polarity
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Schematic illustration of the two weighting coefficients applied during the 3-D pre-stack Kirchhoff migration process. (a) This is a vertical section and shows $\theta_1$, which is the angle between the line connecting the imaging gridpoint and station in the vertical plane. (b) It is a map view and shows $\theta_2$, which is the angle between the line connecting imaging gridpoint and station, and the great circle path. (c) It is a vertical section through the migrated parabola after applying these two cosine terms.</figcaption></figure>

of a converted SV wave is reversed when the angle between the incoming P and converted S is greater than 180°. For teleseismic events arriving at near vertical incidence, this is likely to occur for strongly dipping discontinuities (greater than ~50°). In this case, a positive velocity jump will result in a negative phase in the receiver function. If we do not account for elastic interactions (radiation patterns), this negative phase will be misinterpreted as a negative velocity jump.

It is possible to account for the full elastic effects, by casting the problem in terms of classical linear inverse theory (Bostock & Rondenay 1999; Bostock et al. 2001). The cost of these fully elastic approaches though, is a substantial increase in requirements in terms of data sampling and computational cost compared to the acoustic case (Poppeliers 2001; Pavlis 2011).

In this work, we account for elastic interactions of forward scattered waves (transmissions) by adding a term in the weighting term in eq. (3) in addition to the elements introduced in the acoustic treatment. This factor corrects for the variation in amplitude of the scattered wave as a function of the angle $\theta$ between the incident and scattered (transmitted) wave. This factor is given by the radiation pattern described in Rondenay (2009):

$$ \varepsilon^{P \to S}(r, \theta) = \rho \left( \frac{\delta\beta}{\beta} \left( 2\frac{\beta}{\alpha}\sin2\theta \right) + \frac{\delta\rho}{\rho} \left( \sin\theta + \frac{\beta}{\alpha}\sin2\theta \right) \right) \quad (4) $$

where $\delta\beta/\beta$ is the shear wave velocity perturbation, $\delta\rho/\rho$ is the density perturbation and $\theta$ is the scattering angle (Fig. 5), which is calculated by projecting rays on the great circle plane, using the apparent slowness of the incoming P wave ray, and assuming a straight rays between projected gridpoint and station for the converted S wave. We note that a better estimation of this angle could be obtained by using the gradient of the traveltime field, at higher computational price. Future work will include accounting for fully 3-D elastic patterns (Dahlen et al. 2000). In elastic scattered wave migration/inversion schemes (Bostock & Rondenay 1999; Bostock et al. 2001), velocity and density perturbations are unknown parameters to be solved for. In this work, we simply are interested in the angular dependence, and we fix $\delta\beta/\beta$ to a constant value and assume the density perturbation is zero. The factor $\delta\beta/\beta$ we choose as no importance, as it will act as a scaling factor for the final image. This might also be reasonable when a smoothly changing velocity and density field is met across the whole model, which means the perturbation is small and would not change the pattern much across the whole model. The final migration equation writes:

$$ \begin{aligned} f(r) &= S_{PS}(r) \\ &= \sum_{r_0=1}^{N} \sum_{r_e=1}^{M} \mathrm{RF}(r_e, r_0, t = (\tau_P + \tau_S - \tau_e)) \\ &\quad \times \frac{1}{d} \cos \theta_1 \cos \theta_2 \varepsilon^{P \to S}(r, \theta) \end{aligned} \quad (5) $$
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) Schematic figure (not to scale) showing the earthquake source distribution map used for simulations. A circle of sources surrounds our grid of surface seismic stations (triangles), which is roughly 80° epicentral distance from the sources. A solid triangle corresponds to the seismograms showed in Fig. 4. (b) It is a vertical section though our synthetic velocity model. We test three cases with a horizontal or dipping interface. The angles of dip used were 0°, 30° and 60°.</figcaption></figure>

In the case of a gridpoint where the incident *P* and transmitted *S* are collinear ($\theta$ = 180°), no energy can be transmitted from the *P* wave to the *S* wave, and hence we do not migrate the amplitude of the RF at this gridpoint (the weight equals 0). As will be shown in the next section (Fig. 6), this scattering pattern significantly improves the image, especially for structures with large dipping angle.

# 3 SYNTHETIC TESTS

## 3.1 Model geometry and data

We compare our method with traditional CCP with synthetic seismograms generated from a known velocity model. The true model is defined as two layers of constant velocity separated by a dipping planar discontinuity (Fig. 3), a feature often targeted by RF studies in subduction zones. The top layer has seismic velocities of 3.9 km s⁻¹ for *S* waves and 7.2 km s⁻¹ for *P* waves. The bottom layer has velocities of 4.5 km s⁻¹ for *S* waves and 8.1 km s⁻¹ for *P* waves. In the migration process, the reference model used for traveltime calculation is a smoothed version of the true model around the discontinuity. We smoothed it over two gridpoints to mimic tomographic models.

We present results for three cases, with different angles of dip: 0°, 30° and 60°. At the surface, a 30 by 30 square stations array records the teleseismic events, with a station spacing of 30 km similarly to some dense networks in western North America (e.g. FACES and MENDOCINO experiments; Eakin *et al.* 2010). It can

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Synthetic receiver functions from the station showed with solid triangle in Fig. 3(a) for the three models studied here: (a) 0° dip, (b) 30° dip and (c) 60° dip. Receiver functions are shown between 5 and 80 s, and hence the first *P* pulse is not visible here. The effects of the dipping layer are clearly visible including the azimuthal dependence of traveltimes and polarities.</figcaption></figure>

be shown that ideally the station spacing required to avoid spatial aliasing approaches asymptotically $\lambda/2$ (Rondenay *et al.* 2005) where $\lambda$ is the wavelength. However, this effect tends to diminish with increasing depth and affects the resulting image between the surface and a depth equivalent to approximately twice the station spacing (Rondenay *et al.* 2005). In the synthetic tests shown here, the interstation spacing is 30 km and all the discontinuity structures starts deeper than 50 km, and hence we only stack receiver functions for gridpoints deeper than 50 km. The teleseismic sources used in this study are regularly distributed on a circle around the network, with an average epicentral distance of 80°, as showed in Fig. 3(a).

Synthetic waveforms are constructed with the RAYSUM package (Frederiksen & Bostock 2000), a ray theoretical scheme that allows one to model teleseismic waves in dipping anisotropic structures. Once synthetic seismograms are generated, receiver functions are constructed by deconvolving the vertical component from the radial component with a frequency-domain deconvolution method (Amon 1991). A low-pass Gaussian filter (corresponding to a time pulse of almost 4 s) is applied to stabilize the deconvolution, followed by a bandwidth filter (with cut-offs of 1–0.1 Hz). Fig. 4 shows the synthetic waveforms for all events at one station, for each of the three velocity models.

The main advantage of the proposed scheme is the ability to efficiently compute traveltimes in 3-D with FMM. However, we acknowledge that here we only test our approach in a 2.5-D scenario. This experiment should be seen only as a proof of concept, to demonstrate how our procedure outperforms CCP stacking. To better estimate the resolving power and full potential of our method, future work will consists of testing our scheme in a fully 3-D experiment (with 3-D synthetics computed with a Hybrid waveform modeling code; Monteiller *et al.* 2013). In such 3-D experiments, we will be able to compare results with other approaches based on a 2-D assumptions (e.g. Bostock & Rondenay 1999; Bostock *et al.* 2001; Levander *et al.* 2005). We also plan to assess the resolving power of the method in different situations, for example, by varying data sampling and data noise.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Receiver function images with and without application of the scattering coefficient as a weighting term. (a) Schematic figure shows events used in this simulation, for which all the ray path fall in the negative domain in (c). (b) Receiver functions from the same station as in Fig. 4, constructed with RAYSUM package for the events showed in (a). (c) Vertical section showing the model used in this synthetic test with a 60° dipping layer (dashed line). The positive and negative symbols show the polarity of *P* to *S* energy in receiver functions. (d) Migration result without the scattering polarity correction. (e) Migration result after applying the scattering correction. The solid black triangle in (d) and (e) shows where the array starts at the surface. Note the length of the array is ~900 km. Since there is almost no stacked energy beyond 500 km horizontally, here we only show the cross-section up to 500 km.</figcaption></figure>

## 3.2 The effect of scattering patterns

In order to illustrate the benefit of using elastic scattering patterns to weigh the amplitude of migrated waveforms at conversion points, we first compare results with and without accounting for scattering patterns, that is, acoustic versus elastic migration.

The amplitude of a *Ps* wave converted at a dipping positive discontinuity is reversed for events coming from the downdip direction (Figs 4c and 5b; Bianchi *et al.* 2008). This produces apparent negative discontinuities in acoustic migration stacks as illustrated in Fig. 5(c), for the 60° dipping case.

When all backazimuths are migrated together (Fig. 6), this negative energy coming from the downdip direction may interfere destructively with the positive energy coming from other backazimuths, which produces some artefacts in the final image (Fig. 6b).

The problem of reversed polarities can be addressed by introducing the factor $\epsilon^{P \to S}(r, \theta)$ accounting for elastic scattering patterns. Fig. 5(d) clearly shows that, after this correction, the polarity of the negative phases is flipped, and hence they correctly represent a positive velocity jump in the stacked migrated image.

Another synthetic test is shown in Fig. 6, for which we used events at all azimuths (Fig. 3a). Fig. 6(a) shows the stacked image using the acoustic migration (no polarity correction) for the 60° dipping case. Fig. 6(b) shows the elastic case with the scattering pattern polarity correction. Without the elastic correction (Fig. 6a), there is a negative energy artefact in the upper mantle, above the dipping interface, which disappears when the correction is included (Fig. 6b).

## 3.3 Pre-stack Kirchhoff Migration versus CCP

We compare the results of CCP depth mapping and final FMM pre-stack migration for different dipping interfaces (Fig. 7). Imaging dipping structures is one of the challenges to understand the dynamics of subduction zones (Rondenay 2009). During the CCP stacking process for the receiver functions, the stacking bins are designed to be equally spaced but the width of bins is allowed to vary along the stacking profile according to the data coverage. After the CCP stacking, the stacked receiver functions are simply mapped to
<figure><img src="image_8.png" /><figcaption>**Figure 6.** 3-D pre-stack Kirchhoff migrated image without (a), and with (b) the elastic scattering pattern polarity correction applied as a weighting coefficient. The 60° dipping interface is indicated by the solid black line. The solid black triangle in (a) and (b) shows where the array starts at the surface. Note the length of the array is ~900 km.</figcaption></figure>

the depth domain. The CCP depth images have similar appearance as those of the unmigrated time-domain CCP receiver functions (not shown here). The images obtained by 3-D pre-stack migration are shown along the same profile, which is perpendicular to the strike of the dipping plane. This allows us to directly compare the two approaches and to see their systematic differences.

As shown in Figs 7(a) and (d), both methods recover the discontinuity quite well in the 0° dipping case. However, the pre-stack Kirchhoff migration does not show any artefacts due to multiple reflections, as this energy does not interfere constructively during the stacking process. When interpreting receiver function CCP stacks, it is well known that the image is contaminated by spurious discontinuities due to multiples. Whether and how to interpret the negative energy below the Moho or other positive discontinuities is unclear. The lithospheric signal can be at the same depth as the multiples. Using our pre-stack Kirchhoff migration approach clearly images the discontinuity without the constructive energy from multiples.

In the case of dipping interfaces, the assumption of a horizontal and planar structure adopted in the CCP stacking leads to either wrong information on geometry or distortion of the shape of the discontinuities, such as the shallow dipping structure or defocusing of the diffracted energies (Figs 7b and c). In contrast, the pre-stack migration significantly reduces the unwanted stacking effects and properly accounts for the propagation effects of lateral heterogeneities by making no assumption on the shape of the discontinuities (Figs 7e and f).

In the 30° case, the CCP stacks yield a dip angle smaller than the real model, and a shallower interface, while the pre-stack Kirchhoff migration algorithm is still recovering the interface well. This is because the horizontal layer assumption made in the CCP is not valid, which results in underestimating the dip angle. We also note the presence of an artificial interface in the CCP image, due to multiple reflections.

In the 60° case, the CCP totally fails at recovering the dip angle and the position of the interface. The energy generated by the interface is highly defocused. In contrast, the pre-stack migration coherently recovers the interface down to ~400 km. The maximum depth to which the interface can be imaged is dependent on the array aperture. The signal-to-noise ratio in the resultant image is

proportional to the square root of the number of receiver functions used in stacking (Morozov & Dueker 2003). For real data with noise, since data are stacked, and that noise is not correlated between the different amplitudes stacked, noise will decrease with stacking, and hence effect of noise is also directly proportional to level of data sampling.

## 4 CONCLUSIONS

We have presented a practical and scalable 3-D Kirchhoff pre-stack depth migration scheme for receiver functions analysis of teleseismic arrivals. Receiver functions computed at an array of stations are backpropagated to depth through a 3-D background model with an eikonal solver (the fast marching method) that is significantly computationally cheaper than standard finite-difference schemes.

Synthetic experiments demonstrate that the proposed method allows us to accurately image 3-D complex structures, such as dipping discontinuities. Compared with traditional CCP depth mapping, our migration procedure presents three advantages that significantly improves the quality of the image:

1. By using a 3-D background model for migration, we account for 3-D volumetric heterogeneities, whereas CCP assume a 1-D reference Earth.

2. Dipping discontinuities are accurately mapped.

3. Multiple reflections interfere destructively and do not appear as spurious discontinuities in the final image.

At the same time, compared with more sophisticated finite-difference migration schemes (Bostock et al. 2001; Poppeliers & Pavlis 2003a,b), our method only requires computations of travel-times, and hence is conceptually simpler, easier to implement, and computationally cheaper. In this way, it has the potential to be applied to large arrays of stations (e.g. USArray) for seismic imaging at the continental scale.

We also acknowledge that Kirchhoff methods have drawbacks. Identifying appropriate imaging conditions (weighting terms) may be difficult. The effects of noise in the data and inaccurate velocity models (both inaccurate velocities and complexities of the real Earth that are not accounted for in the smooth model) and sparse
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Receiver function images using CCP depth migration (a)–(c), and 3-D pre-stack Kirchhoff depth migration (d)–(f) for our three synthetic models with dips of 0° (a) & (d), 30° (b) & (e) and 60° (c) & (f). Solid black lines mark the real location of the model discontinuity. The solid black triangle in all the figures shows where the array starts at the surface. Note the length of the array is ~900 km. Since there is almost no stacked energy beyond 500 km horizontally, here we only show the cross-section up to 500 km.</figcaption></figure>
(low-fold) coverage, lead to poor performance of Kirchhoff methods. A problem here is that most earthquake data are very low-fold, and these issues may be non-negligible for further applications. This first work only presents a proof of concept and future work will include testing our migration approach in fully 3-D structures, quantifying resolution and testing the limits in terms of data coverage, station spacing, data noise, etc.

Also one of the great advantages of Kirchhoff methods for earthquake applications is that the receivers do not need to be evenly spaced. However, large gaps between receivers will produce artefacts or low amplitudes or unrecognizable features in the model. This could possibly be addressed using least-squares migration scheme (Nemeth et al. 1999) in the future.

We expect this method to be a powerful tool for imaging subduction zones where the key tectonic structure has a large dip. Accurate imaging will put constraints on the slab location and material pathways into the mantle and also can provide prior information for geodynamic modeling of subducted plate. The method will likely be most helpful in regions with large slab dips, such as the Cascadia subduction zone, Japan, Central America and the Aegean.

## ACKNOWLEDGEMENTS

This work was supported by the University of California, Berkeley. Thomas Bodin wishes to acknowledge support from the Miller Institute for Basic Research at the University of California, Berkeley. We thank the editor, Prof Randy Keller and the two anonymous reviewers for their constructive comments which improve the quality of the paper.

## REFERENCES

Abers, G., 1998. Array measurements of phase used in receiver-functions calculations: importance of scattering, *Bull. seism. Soc. Am.*, **88**, 313–318.

Abt, D., Fischer, K.M., French, S.W., Ford, H.A., Yuan, H.Y. & Romanowicz, B., 2010. North American lithospheric discontinuity structure imaged by Ps and Sp receiver functions, *J. geophys. Res.*, **115**, B09301, doi:10.1029/2009JB006914.

Ai, Y., Chen, Y., Zeng, F., Hong, X. & Ye, W., 2007. The crust and upper mantle structure beneath southeastern China, *Earth planet. Sci. Lett.*, **260**(2007), 549–563

Ammon, C.J., 1991. The isolation of receiver effects from teleseismic P waveforms, *Bull. seism. Soc. Am.*, **81**, 2504–2510.

Bianchi, I., Agostinetti, N.P., Gori, P.D. & Chiarabba, C., 2008. Deep structure of the Colli Albani volcanic district (central Italy) from receiver functions analysis, *J. geophys. Res.*, **113**, B09313, doi:10.1029/2007JB005548.

Biondi, B. & Sava, P., 1999. Wave-equation migration velocity analysis, in *69th Annual International Meeting, Society of Exploration Geophysicists, Expanded Abstracts*, pp. 1723–1726.

Bodin, T., Yuan, H.Y. & Romanowicz, B., 2014. Inversion of receiver functions without deconvolution, *Geophys. J. Int.*, **196**(2), 1025–1033.

Bostock, M.G., 2007. Teleseismic body-wave scattering and receiver-side structure, in *Treatise on Geophysics*, Vol. 1, pp. 219–246.

Bostock, M.G. & Rondenay, S., 1999. Migration of scattered teleseismic body waves, *Geophys. J. Int.*, **137**, 732–746.

Bostock, M.G., Rondenay, S. & Shragge, J., 2001. Multiparameter two-dimensional inversion of scattered teleseismic body waves, 1. Theory for oblique incidence, *J. geophys. Res.*, **106**, 30 771–30 782.

Burdick, L. & Langston, C., 1977. Modeling crustal structure through the use of converted phases in teleseismic body-wave forms, *Bull. seism. Soc. Am.*, **67**(3), 677–691.

Caldwell, W.B., Klemperer, S.L., Lawrence, J.F. & Ashish, S.S.R., 2013. Characterizing the Main Himalayan Thrust in the Garhwal Himalaya,

India, with receiver function CCP stacking, *Earth planet. Sci. Lett.*, **367**, 15–27.

Chen, L., Wen, L.X. & Zheng, T.Y., 2005. A wave equation migration method for receiver function imaging, (I) Theory, *J. geophys. Res.*, **110**, B11309, doi:10.1029/2005JB003665.

Chevrot, S. & van der Hilst, R.D., 2000. The Poisson ratio of the Australian crust: geological and geophysical implications, *Earth planet. Sci. Lett.*, **183**(1), 121–132.

Claerbout, J.F., 1985. *Imaging the Earth's Interior*, Blackwell Scientific Publications.

Dahlen, F.A., Nolet, G. & Hung, S.-H., 2000. Fréchet kernels for finite-frequency traveltime—I. Theory, *Geophys. J. Int.*, **141**, 157–174.

Dueker, K.G. & Sheehan, A.F., 1997. Mantle discontinuity structure from midpoint stacks of converted P to S waves across the Yellowstone hotspot track, *J. geophys. Res.*, **102**, 8313–8327.

Eakin, C.M., Obrebski, M., Allen, R.M., Boyarko, D.C., Brudzinski, M.R. & Porritt, R., 2010. Seismic anisotropy beneath Cascadia and the Mendocino triple junction: interaction of the subducting slab with mantle flow, *Earth planet. Sci. Lett.*, **297**, 627–632.

Frederiksen, A.W. & Bostock, M.G., 2000. Modelling teleseismic waves in dipping anisotropic structures, *Geophys. J. Int.*, **141**, 401–412.

Frederiksen, A.W. & Revenaugh, J., 2004. Lithospheric imaging via teleseismic scattering tomography, *Geophys. J. Int.*, **159**, 978–990.

Julià, J., Ammon, C.J., Herrmann, R.B. & Correig, A.M., 2000. Joint inversion of receiver function and surface wave dispersion observations, *Geophys. J. Int.*, **143**, 1–19.

Julià, J., Ammon, C.J. & Herrmann, R.B., 2003. Evaluation of deep sediment velocity structure in the New Madrid Seismic Zone, *Bull. seism. Soc. Am.*, **94**, 334–340.

Kind, R., Kosarev, G. & Petersen, N., 1995. Receiver functions at the stations of the German Regional Seismic Network (GRSN), *Geophys. J. Int.*, **121**(1), 191–202.

Langston, C.A., 1979. Structure under Mount Rainier, Washington, inferred from teleseismic body waves, *J. geophys. Res.*, **84**, 4749–4762.

Lekić, V. & Fischer, K.M., 2014. Contrasting lithospheric signatures across the western United States revealed by Sp receiver functions, *Earth planet. Sci. Lett.*, **402**, 90–98.

Levander, A. & Miller, M., 2012. Evolutionary aspects of lithosphere discontinuity structure in the western U.S., *Geochem. Geophys. Geosyst.*, **13**(1), Q0AK07, doi:10.1029/2012GC004056.

Levander, A., Niu, F. & Symes, W.W., 2005. Imaging teleseismic P and S scattered waves using the Kirchhoff integral, in *Seismic Earth: Ar-ray Analysis of Broadband Seismograms*, Vol. 157 in AGU Geophysical Monograph, pp. 149–169, eds Levander, A. & Nolet, G., AGU, Washington, DC.

Liu, Z. *et al.*, 2015. Receiver function images of the mantle transition zone beneath NE China: new constraints on intraplate volcanism, deep subduction and their potential link, *Earth planet. Sci. Lett.*, **412**, 101–111.

Lombardi, D., Braunmiller, J., Kissling, E. & Giardini, G., 2008. Moho depth and Poisson's ratio in the Western-Central Alps from receiver functions, *Geophys. J. Int.*, **173**, 249–264.

Monteiller, V., Chevrot, S., Komatitsch, D. & Fuji, N., 2013. A hybrid method to compute short period synthetic seismograms of teleseismic body waves in a 3-D regional model, *Geophys. J. Int.*, **192**, 230–247.

Morozov, I.B., 2004. Crustal scattering and some artefacts in receiver function images, *Bull. seism. Soc. Am.*, **94**(4), 1492–1499.

Morozov, I.B. & Dueker, K.G., 2003, Signal-to-noise ratios of teleseismic receiver functions and effectiveness of stacking for their enhancement, *J. geophys. Res.*, **108**(B2), 2106, doi:10.1029/2001JB001692.

Nemeth, T., Wu, C. & Schuster, G.T., 1999. Least-squares migration of incomplete reflection data, *Geophysics*, **64**, 208–221.

Pavlis, G.L., 2011. Three-dimensional, wavefield imaging of broadband seismic array data, *Comput. Geosci.*, **37**(8), 1054–1066.

Phinney, R., 1964. Structure of the Earth's crust from spectral behavior of long-period body waves, *J. geophys. Res.*, **69**, 2997–3017.

Poppeliers, C., 2001. Prestack planewave migration of teleseismic P-to-S converted phases, *PhD dissertation*, Indiana University, Bloomington.
Poppeliers, C. & Pavlis, G.L., 2003a. Three-dimensional, prestack, plane wave migration of teleseismic P-to-S converted phases: 1. Theory, *J. geophys. Res.*, **108**, 2112, doi:10.1029/2001JB000216.

Poppeliers, C. & Pavlis, G.L., 2003b. Three-dimensional, prestack, plane wave migration of teleseismic P-to-S converted phases: 2. Stacking multiple events, *J. geophys. Res.*, **108**, 2267, doi:10.1029/2001JB001583.

Rawlinson, N. & Sambridge, M., 2004. Wave front evolution in strongly heterogeneous layered media using the fast marching method, *Geophys. J. Int.*, **156**(3), 631–647.

Revenaugh, J., 1995. A scattered-wave image of subduction beneath the transverse ranges, *Science*, **268**, 1888–1892.

Rondenay, S., 2009. Upper mantle imaging with array recordings of converted and scattered teleseismic waves, *Surv. Geophys.*, **30**, 377–405.

Rondenay, S., Bostock, M.G. & Shragge, J., 2001. Multiparameter two-dimensional inversion of scattered teleseismic body waves, 3. Application to the Cascadia 1993 data set, *J. geophys. Res.*, **106**, 30795–30808.

Rondenay, S., Bostock, M. & Fischer, K., 2005. Multichannel inversion of scattered teleseismic body waves: practical considerations and applicability, in *Seismic Earth: Array Analysis of Broadband Seismograms*, Vol. 157, of AGU Geophysical Monograph, pp. 187–204, eds Levander, A. & Nolet, G., AGU, Washington, DC.

Ryberg, T. & Weber, M., 2000. Receiver function arrays: a reflection seismic approach, *Geophys. J. Int.*, **141**, 1–11.

Sandvol, E., Seber, D., Calvert, A. & Barazangi, M., 1998. Grid search modeling of receiver functions: implications for crustal structure in the Middle East and North Africa, *J. geophys. Res.*, **103**(B11), 26899–26917.

Sava, P. & Biondi, B., 2004a. Wave-equation migration velocity analysis—I: Theory, *Geophys. Prospect.*, **52**, 593–606.

Sava, P. & Biondi, B., 2004b. Wave-equation migration velocity analysis—II: Subsalt imaging examples, *Geophys. Prospect.*, **52**, 607–623.

Schmandt, B., Dueker, K., Humphreys, E. & Hansen, S., 2012. Hot mantle upwelling across the 660 beneath Yellowstone, *Earth planet. Sci. Lett.*, **331**, 224–236.

Schneider, W.A., 1977. Integral formulation for migration in two and three dimensions, *Geophysics*, **43**, 49–76.

Sethian, J.A., 1996. A fast marching level set method for monotonically advancing fronts, *Proc. Natl. Acad. Sci. USA*, **93**, 1591–1595.

Sethian, J.A. & Popovici, A.M., 1999. 3-D traveltime computation using the fast marching method, *Geophysics*, **64**, 516–523.

Shang, X., de Hoop, M. & van der Hilst, R., 2012. Beyond receiver functions: passive source reverse time migration and inverse scattering of converted waves, *Geophys. Res. Lett.*, **39**, 1–7.

Sheehan, A.F., Shearer, P.M., Gilbert, H.J. & Dueker, K.G., 2000. Seismic migration processing of P–SV converted phases for mantle discontinuity structure beneath the Snake River Plain, western United States, *J. geophys. Res.*, **105**, 19055–19065.

Shen, W., Ritzwoller, M.H., Schulte-Pelkum, V. & Lin, F.-C., 2013. Joint inversion of surface wave dispersion and receiver functions: a Bayesian Monte-Carlo approach, *Geophys. J. Int.*, **192**(2), 807–836.

Shi, D., Wu, Z., Klemperer, S.L., Zhao, W., Xue, G. & Su, H., 2015. Receiver function imaging of crustal suture, steep subduction, and mantle wedge in the eastern India–Tibet continental collision zone, *Earth planet. Sci. Lett.*, **414**, 6–15.

Shragge, J., Artman, B. & Wilson, C., 2006. Teleseismic shot-profile migration, *Geophysics*, **71**(4), SI221–SI229.

Spieker, K., Wölbern, I., Thomas, C., Harnafi, M. & Moudnib, L.E., 2014. Crustal and upper-mantle structure beneath the western Atlas Mountains in SW Morocco derived from receiver functions, *Geophys. J. Int.*, **198**(3), 1474–1485.

Steckler, M.S., Agostinetti, N.P., Wilson, C.K., Roselli, P., Seeber, L., Amato, A. & Lerner-Lam, A., 2008. Crustal structure in the southern Apennines from teleseismic receiver functions, *Geology*, **36**(2), 155–158.

Tauzin, B., van der Hilst, R.D., Wittlinger, G. & Ricard, Y., 2013. Multiple transition zone seismic discontinuities and low velocity layers below western United States, *J. geophys. Res.*, **118**, 2307–2322.

Vinnik, L., 1977. Detection of waves converted from P to SV in the mantle, *Phys. Earth planet. Inter.*, **15**(1), 39–45.

Wiggins, J.W., 1984. Kirchhoff integral extrapolation and migration of non-planar data, *Geophysics*, **49**, 1239–1248.

Yu, C.Q., Chen, W.P. & van der Hilst, R.D., 2013. Removing source-side scattering for virtual deep seismic sounding (VDSS), *Geophys. J. Int.*, **195**(2013), pp. 1932–1941.

Zheng, T.Y., Zhao, L., He, Y.M. & Zhu, R.X., 2014. Seismic imaging of crustal reworking and lithospheric modification in eastern China, *Geophys. J. Int.*, **196**(2), 656–670.

Zhu, L.P. & Kanamori, H., 2000. Moho depth variation in southern California from teleseismic receiver functions, *J. geophys. Res.*, **105**, 2969–2980.