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

HAL
open science

# Inversion of receiver functions without deconvolution-application to the Indian craton

Thomas Bodin, Huaiyu Yuan, Barbara Romanowicz

► To cite this version:

Thomas Bodin, Huaiyu Yuan, Barbara Romanowicz. Inversion of receiver functions without deconvolution-application to the Indian craton. Geophysical Journal International, 2014, 196, pp.1025-1033. 10.1093/gji/ggt431. insu-03581149

HAL Id: insu-03581149

https://insu.hal.science/insu-03581149

Submitted on 19 Feb 2022

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_2.png" /></figure>

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

# Inversion of receiver functions without deconvolution—application to the Indian craton

Thomas Bodin,<fnref n="1" /> Huaiyu Yuan<fnref n="1" /> and Barbara Romanowicz<fnref n="1" />,<fnref n="2" />,<fnref n="3" />

Accepted 2013 October 11. Received 2013 August 15; in original form 2013 January 30

## SUMMARY

Receiver functions are a powerful tool to isolate and interpret receiver-side structure effects in teleseismic seismic records. They are easily constructed by deconvolving one component of a seismogram by another. Deconvolution is the inverse of convolution, and hence can be mathematically viewed as an inverse problem. It is a numerically unstable procedure that needs to be stabilized (i.e. regularized). This points to a recurring problem in geophysical imaging: there is a trade-off between variance and resolution, where the user needs to arbitrarily define a level of compromise. Here we propose a novel misfit function for inversion of converted phases that avoids deconvolution. In this way, the choice of regularization parameters (e.g. water level, width of a low pass filter) is avoided, and statistics of data errors can be correctly accounted for. We use this misfit measure to construct a likelihood probability function and carry out a transdimensional Bayesian inversion for shear wave structure. After illustrating the method with a synthetic test, a real data application is shown where teleseismic signals recorded at HYB station (Hyderabad, India) and surface wave dispersion measurements are jointly inverted to provide a probabilistic 1-D seismic model beneath the station. The results help address the debate on the thickness of the lithosphere in this region. We show that the sharp negative velocity jump at 110 km that was previously interpreted as the lithosphere–asthenosphere boundary (LAB) is actually a mid-lithospheric discontinuity. The actual LAB is seen deeper as a milder gradient between 150 and 200 km.

**Key words:** Time-series analysis; Inverse theory; Probability distributions; Body waves; Surface waves and free oscillations; Cratons.

## 1 INTRODUCTION

The coda of teleseismic P waves contains a large number of phases, generated at interfaces beneath the receiver, that contain a significant amount of information on seismic structure. However, these phases are difficult to identify as they are buried in microseismic noise, and convolved with the source–time function (Fig. 1). Hence, the vertical V(t) and horizontal (radial) H(t) components of a seismogram for a plane P wave can be written using convolutions:

$$V(t) = s(t) * v(t) \qquad (1)$$

$$H(t) = s(t) * h(t), \qquad (2)$$

where s(t) is the source–time function, which may be quite complicated since it is related to dislocation time history and source area reverberations, and v(t) and h(t) are the vertical and radial impulse response functions of the near receiver structure.

The problem of isolating the structure effect is overcome by a method developed in the 1970s following the pioneering work by

Phinney (1964) now widely used in seismology. The idea is to deconvolve the vertical component from the horizontal components to produce a time series called a 'receiver function' (RF; Vinnik 1977; Burdick & Langston 1977; Langston 1979).

$$R_{\text{obs}}(t) = \frac{\mathbf{H}(t)}{V(t)} * \mathbf{G}(t) = \frac{\mathbf{h}(t)}{\mathbf{v}(t)} * \mathbf{G}(t), \qquad (3)$$

where the fraction refers to a deconvolution (or spectral division). In this way the influence of source and distant path effects are eliminated, and hence one can enhance conversions from P to S generated at boundaries beneath the recording site. Note that the receiver function is smoothed with a Gaussian filter G(t) in order to eliminate high frequency errors introduced during the deconvolution (Fig. 2). The smoothed receiver function waveform can be directly interpreted by visual inspection (e.g. Li et al. 2002; Audet et al. 2009; Nikulin et al. 2009; Abt et al. 2010; Tauzin et al. 2010; Levander et al. 2011; Leahy et al. 2012; Rychert et al. 2013). However, the waveform needs to be migrated to depth (Kosarev et al. 1999; Bostock 2002) which requires knowledge of the velocity model,

<footnote n="1">¹Berkeley Seimological Laboratory, 215 McCone Hall, UC Berkeley, Berkeley CA 94720-4760, USA. E-mail: thomas.bodin@berkeley.edu</footnote>

<footnote n="2">²College de France, 11 place Marcelin Berthelot, F-75231 Paris, France</footnote>

<footnote n="3">³Institut de Physique du Globe de Paris (IPGP), 1 Rue Jussieu, F-75005 Paris, France</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 1. Synthetic waveform computed for an incident plane P-wave arriving with incidence angle of 27° (ray parameter = 0.066 s km⁻¹), and traversing a layered earth model, where the last layer is an infinite half-space. The source-time function is a box car function (1 s duration) and some white Gaussian noise has been added with standard deviation of 3 per cent the maximum amplitude. For simplicity, we have not convolved the synthetic waveform with an instrument response function, and acknowledge that the resulting signal might seem different from an observed seismogram.</figcaption></figure>

and interpretation might be subjective. Furthermore, in the case of P-receiver functions, the solution is blurred with multiple reflections from the surface that produce apparent discontinuities (Kind et al. 2012).

Alternatively, receiver functions can be directly inverted in the time domain for an S-wave velocity model of the crust and uppermost mantle beneath the receiver (Kind et al. 1995; Sandvol et al. 1998; Vinnik et al. 2006; Darbyshire et al. 2009; Piana Agostinetti & Malinverno 2010; Stipčević et al. 2011; Shen et al. 2013; Sriniyas et al. 2013). In this case a misfit function is constructed by comparing the observed receiver function with a receiver function predicted for some earth model.

‘Migration’ and ‘inversion’ of receiver functions are two different ways to interpret converted phases in terms of seismic structure. In both cases interpretation is non-unique, and some information on absolute velocities needs to be added. In the case of migration, a background velocity model is used to migrate waveforms. For inversion, surface wave dispersion measurements are typically inverted jointly with receiver functions. Both approaches have advantages and drawbacks. Inversion is more quantitative, as a forward model will take into account all possible reflections and conversions, whereas caution and experience is needed when qualitatively interpreting migrated waveforms. However, in the migration approach the user is directly looking at data, which for example allows one to get an intuitive feeling of the level of noise. Finally, migration techniques are computationally cheap and allow construction of 2-D or 3-D sections. Inversions are usually done with expensive sampling methods (e.g. Monte Carlo), and hence are limited for the moment to the 1-D case.

Although both migration and inversion of RFs have been extensively used for the last 30 yr, there are two well known drawbacks: (1) The deconvolution is an unstable numerical procedure that needs to be damped. This results in a loss of resolution and introduces errors in the receiver function, with a trade-off between the two effects as shown in Fig. 2. (2) It is difficult to estimate the

<figure><img src="image_5.png" /><figcaption>Figure 2. Receiver functions computed from waveforms in Fig. 1 for two Gaussian filters, with parameters a = 2 (top), and a = 6 (bottom), compared with ‘true’ receiver functions. True receiver functions are computed with a Dirac source function and with no noise added to the waveform. This illustrates the trade-off between resolution and variance. Smoothing the waveform stabilizes the deconvolution at the cost of losing resolution. There is no mathematical model for the noise introduced by deconvolution.</figcaption></figure>

nature and level of uncertainties in the observed receiver function. That is, there is no clear way to theoretically propagate the noise present in different components of the seismogram into errors in the deconvolved waveform (Di Bona et al. 1998). These two issues have been well documented in the literature, and have been the subject of much work in the last decades (Oldenburg 1981; Ligorria & Ammon 1999; Park & Levin 2000; Helffrich 2006; Gurrola et al. 2007). Following ideas of Menke & Levin (2003), here we present a novel approach to RF inversion that avoids deconvolution, and hence directly overcomes these two problems.

## 2 A CROSS-CONVOLUTION MISFIT FUNCTION

Here we use the cross-convolution misfit function that was proposed by Menke & Levin (2003) for inversion of SKS splitting measurements. If $v_p(t, m)$, and $h_p(t, m)$ are predicted structure response functions for a given earth model $m$, we can convolve eq. (1) by $h_p(t, m)$ and (2) by $v_p(t, m)$:

$$h_p(t, m) * V(t) = s(t) * v(t) * I(t) * h_p(t, m) \quad (4)$$

$$v_p(t, m) * H(t) = s(t) * h(t) * I(t) * v_p(t, m). \quad (5)$$

The misfit function is then defined as the difference between the left-hand sides of (4) and (5).

$$\Phi(m) = \|v_p(t, m) * H(t) - h_p(t, m) * V(t)\|^2 \quad (6)$$

which is minimized when $v_p(t, m) = v(t)$ and $h_p(t, m) = h(t)$.

This misfit function is equivalent to the distance between the observed and predicted receiver functions. It does not make the inverse problem more linear or more unique. However, (1) it does not
require any deconvolution procedure, no damping parameter needs to be chosen, and hence no processing errors are introduced (2) the chi-squared $\chi^2$ (or log- likelihood) probability density function can be easily derived from errors statistics in seismograms $V(t)$, and $H(t)$.

Since discrete convolution in time is a simple summation, and since seismograms can be seen as corrupted by random errors, each sample of the signal obtained after discrete convolution is then a sum of random variables, whose statistics are straightforward to calculate with algebra of random variables. This is not the case with deconvolution schemes.

For example, let us assume as shown in Fig. 1 a simple case where $V(t)$, and $H(t)$ contain independent, and normally distributed random errors (i.e. Gaussian white noise) with standard deviation $\sigma$. It can be shown easily that the chi-squared statistic for the misfit function $\Phi(\mathbf{m})$ is:

$$ \chi^2 = \frac{\Phi(\mathbf{m})}{\sigma^2 \left[ \sum_{i=1}^{n} (\mathbf{h}_p^i)^2 + \sum_{i=1}^{n} (\mathbf{v}_p^i)^2 \right]}, \quad (7) $$

where $n$ is the number of samples in the signal. Note that the chi-squared function can also be computed for more complicated noise models, such as correlated noise in time (which occurs when filtering waveforms), or correlation between components (Yuan et al. 2008).

### 3 BAYESIAN INVERSION

A correct mathematical form for the chi-squared distribution allows us to properly write the likelihood function which measures the probability that the predicted and observed data are consistent given a mathematical model for the random noise distribution (Box & Tiao 1973; Smith 1991; Gelman et al. 1995; Sivia 1996).

Assessment of variability in the data and in the misfit function is not indispensable in optimization based inversion, that is, where one only seeks the maximum likelihood model. This is because uncertainty estimates do not affect the best-fitting solution, or peak of the likelihood function. For example, the solution to a linear regression does not change when changing the length of error bars around data points (Bodin et al. 2012b).

However, noise levels and statistics become crucial in a Bayesian sampling framework, where the goal is to quantify the full range of possible models compatible with data. This is because data uncertainty estimates strongly determine the shape (and width) of the likelihood function in model space (Dettmer et al. 2007, 2009). In this way, the cross-convolution misfit function enables us to carry out a proper Bayesian inversion which correctly propagates errors in the data toward model uncertainties. Data are taken as they are, with no need of arbitrarily defined stabilization parameters (i.e. water level, or width of Gaussian filters), which may bias the solution in a statistical sense.

We use a Bayesian formalism, and tackle the problem probabilistically. The goal is to estimate the posterior probability distribution $p(\mathbf{m}|\mathbf{d})$, which describes the probability of having a discontinuous model $\mathbf{m}$ given the observed measurements $\mathbf{d}=[V(t), H(t)]$ (Tantotola & Valette 1982; Duijndam 1988a,b). Bayes' theorem (Bayes 1763) is used to combine prior information on the model with the observed data to give the posterior probability density function:

$$ \text{posterior} \propto \text{likelihood} \times \text{prior} \quad (8) $$

$$ p(\mathbf{m} | \mathbf{d}) \propto p(\mathbf{d} | \mathbf{m})p(\mathbf{m}), \quad (9) $$

where $p(\mathbf{m})$ is the a priori probability density of $\mathbf{m}$, and describes the level of knowledge we have about seismic structure before considering $\mathbf{d}$. In this study we set priors to uniform distribution with relatively wide bounds, and hence impose little constraints to the final solution. The likelihood function $p(\mathbf{d} | \mathbf{m})$ is the probability of observing the data given a particular model. The form of this probability density function is given by what we think about uncertainties on $\mathbf{d}$. That is, the form of the error statistics for $\mathbf{d}$ must be assumed to derive $p(\mathbf{d} | \mathbf{m})$. Assuming independent, and normally distributed random errors for the two components of the seismogram, the likelihood distribution for our cross-convolution misfit function in (6) writes

$$ p(\mathbf{d} | \mathbf{m}) = \frac{\exp\left(-\frac{\chi^2}{2}\right)}{\sqrt{[2\pi\sigma^2 (\sum_{i=1}^n (\mathbf{h}_p^i)^2 + \sum_{i=1}^n (\mathbf{v}_p^i)^2)]^n}} \quad (10) $$

Here we sample the posterior distribution with the reversible jump algorithm (Green 1995, 2003) as implemented in Bodin et al. (2012c). For examples of recent application of this algorithm in Earth sciences, see Malinverno (2002), Hopcroft et al. (2009), Gallagher et al. (2011), Luo (2010), Dettmer et al. (2010), Piana Agostinetti & Malinverno (2010), Minsley (2011), Ray & Key (2012), Iaffaldano et al. (2012), Bodin et al. (2012a), Tkalčić et al. (2013), Sambridge et al. (2013), Young et al. (2013) and Iaffaldano et al. (2013).

This class of algorithm treats the number of model parameters (i.e. number of layers) as an unknown in the problem. Of course the data can be fit better as we increase the dimension of the model (i.e. the number of variable parameters), but the procedure naturally prevents the data to be fitted more than the given level of noise. That is, between two models that fit the data equally well, the simpler one (in terms of number of parameters) will be favoured in Bayesian sampling (MacKay 2003). A Markov chain Monte Carlo (MCMC) scheme is used to generate a collection of S-wave velocity profiles $\mathbf{m}_i$ with variable number of layers, whose statistical distribution is proportional to the posterior distribution. For a detailed description of the model parameterization and the algorithm, we refer the reader to Bodin et al. (2012c).

This direct search algorithm requires solving the forward problem a large number of times. We use the Thomson-Haskell matrix method (Thomson 1950; Haskell 1953) to compute the spectral response of a stack of isotropic layers to an incident planar P wave. Multiple conversions and reflections are considered with this method, which enables us to properly model the full P code. Although this algorithm does not account for anisotropy or dipping layers, it is fast and has been widely used in Monte Carlo algorithms (e.g. Shibutani et al. 1996; Sambridge 1999).

We first inverted synthetic waveforms calculated after propagating a plane wave through the model in Fig. 1. Three cases with different ray parameters were tested. The same amount of noise was added in the three cases (i.e. 3 per cent of maximum amplitude of the vertical component). Results are shown in Fig. 3. For each case, transdimensional sampling was carried out allowing between 2 and 35 interfaces. Bounds for the uniform prior distribution were set to 2.3–4.9 km s⁻¹ for S-wave velocity values.

Posterior inference was made using an ensemble of $7 \times 10^6$ 1-D models with a density representing the probability of the model given the data. The density plot of the ensemble allows us to visualize the posterior probability for $V$s at each depth (left-hand panels of Fig. 3). The resolution of seismic discontinuities can be examined by plotting the marginal distribution on the location of
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Synthetic waveforms, obtained after propagating a plane wave through the model in Fig. 1, are inverted for a layered $V_s$ model with a transdimensional Bayesian Inversion using cross-convolution. Three separate inversions are carried out for different incidence angles. (A and B) Ray parameter at 0.08 s km$^{-1}$. (C and D) Ray parameter at 0.066 s km$^{-1}$. (E and F) Ray parameter at 0.04 s km$^{-1}$. Data used for the middle panels is shown in Fig. 1. The solution is a large ensemble of 1-D profiles representing the probability of the model given the data. Left-hand panels: the density plot of this ensemble of models shows the probability for $V_s$ at each depth (the true model is shown in black). Right-hand panels: probability of having a discontinuity at each depth.</figcaption></figure>

inter interfaces in the ensemble of sampled models (right-hand panels). Note how the number and location of seismic discontinuities are well recovered.

In the top panels the plane wave is arriving with a large incidence angle (large ray parameter), conversions are strong, and the structure is well recovered. As the incidence angle decreases (middle and bottom panels), conversions become weaker, and the ability to resolve structure vanishes to zero. This effect is well seen by the algorithm as posterior distributions become flatter and wider when the ray parameter decreases. This illustrates the ability of the Bayesian algorithm to quantify uncertainties.

# 4 APPLICATION TO DATA FROM HYB STATION, INDIA

## 4.1 Data

In order to demonstrate how this novel approach to receiver-side structure imaging fares on real data, we apply it to waveforms recorded at the Hyderabad station (HYB). The station belongs to the Geoscope network since 1987 (Romanowicz et al. 1991), and hence provides a large number of useful recordings (Fig. 4). The site is located on the late Archean granite-gneiss terrain of the eastern Dharwar craton, formed at 2.5 Gyr and stabilized in early Proterozoic (Kiselev et al. 2008; Sarkar et al. 2003). A total of $N = 323$ events collected between 1997 and 2007 were used with backazimuths between 255° and 310° and with epicentral distances between 37° and 70° (i.e. ray parameters between 0.056 and 0.075 s km$^{-1}$). Only events with magnitude 6.0 were used.

Seismograms were cut for the same time window (starting 12 s before first P arrival, and with a total length duration of 72 s), normalized to equal energy, and rotated to radial and tangential components. We note, as explained below, that events do not need to be aligned to theoretical arrival, or maximum amplitude.

The observed vertical and horizontal waveforms needed for inversion were obtained by simply stacking all events. Influence of the receiver structure is common to all records and is enhanced by summation (Shearer 1991; Kind et al. 2012). Hence we have:

$$V(t) = \sum_{i=1}^{N} V_i(t) \qquad (11)$$

$$H(t) = \sum_{i=1}^{N} H_i(t), \qquad (12)$$

where the recorded signal for each event [$V_i(t)$, $H_i(t)$] can be written as the product of the structure and the source $s_i(t)$ plus some random noise $n$:

$$V_i(t) = s_i(t) * v(t) + n_i(t) \qquad (13)$$

$$H_i(t) = s_i(t) * h(t) + n_i(t). \qquad (14)$$

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Map showing HYB station (blue square) and the 323 events used to construct the observed waveform [$V(t)$, $H(t)$].</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Waveforms obtained by stacking all events in Fig. 4, and used for Bayesian inversion in Figs 6 and 7.</figcaption></figure>

For simplicity we have dropped the instrument response term. Since the convolution is a distributive law, (11) and (12) can be written:

$$V(t) = v(t) * \sum_{i=1}^{N} s_i(t) + \sum_{i=1}^{N} n_i(t) \qquad (15)$$

$$H(t) = h(t) * \sum_{i=1}^{N} s_i(t) + \sum_{i=1}^{N} n_i(t) \qquad (16)$$

with the second term vanishing to zero as the number N of events increases. Indeed, if $n_i(t)$ is a vector of Gaussian random numbers with variance $\sigma^2$, the variance of our waveform [$V(t)$, $H(t)$] is $\sigma^2/N$. In this way, by directly averaging measured waveforms with different source-time functions, we obtain a signal with enhanced signal-to-noise ratio (Fig. 5), that would be recorded for a source equal to the average of all sources. Since the convolution operator is linear, there is no need to adjust timing of different events and align waveforms. However, it is clear that alignment of events by maximum amplitude will help increase the signal-to-noise ratio (e.g. Shearer 1991; Kumar et al. 2010). Although alignment (or not) of events does not change the nature of the inverse problem (non-linearity, resolution), this will change the data uncertainty, and hence the final model uncertainty.

In receiver function analysis, the time arrival of converted phases is a function of slowness (called the moveout). Here our forward model is computed for a given ray parameter, and hence we only invert rays arriving around this given slowness (i.e. 0.056–0.0755 s km⁻¹). We note that at increased computational cost, one could jointly invert different ‘bins’ for different ray parameters. Note also that (15) and (16) hold only if the receiver-side structure sampled by each event is identical. We therefore choose only events in a given range of back azimuths (i.e. 255–310°). In this way we neglect possible Moho dip and anisotropy, which are only second-order effects in the context of deriving 1-D models of the crust and upper mantle. As we increase the quadrant defining the range of ray parameters and backazimuths, more events are included which

results in ambient and instrumental noise reduction. However, 3-D and moveout effects also become more significant which may result in data incoherence, particularly for deep conversions. Hence there is a trade-off between ‘data noise’ and ‘theoretical noise’ due to 3-D structure (i.e. errors in the forward model). A compromise needs to be found when defining the range of incident rays. For each experiment, the optimal interval depends on the number and distribution of available events, as well as the level of 3-D heterogeneities and anisotropy. Here the range of ray parameters and backazimuths was chosen by an ad hoc trial and error procedure.

## 4.2 Results

### 4.2.1 Inversion of converted body waves

The posterior solution obtained with the transdimensional Bayesian algorithm is shown in Fig. 6. Results are in strong agreement with Kiselev et al. (2008) who jointly inverted RF data with traveltime residuals in the same region. The crust is well resolved with both positive and negative transitions. A sharp Moho is clearly visible at 30 km depth, with an S velocity positive jump within the crust at a depth of 25 km. The H discontinuity (Revenaugh & Jordan 1991) is well imaged at 75 km. This is most plausibly attributed the transition to from spinel to garnet facies in aluminous peridotite (Hales 1969). The inverted velocity structure correlates also well with both P–S and S–P migrated receiver functions shown in Kumar et al. (2007) and Saul et al. (2000) for the same station. They observed a high velocity mantle lid with an underlying low velocity zone and a

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Probabilistic solution obtained with the trans-dimensional Bayesian algorithm under HYB station, using P to S converted phases only. Density plot of the ensemble of models in the solution showing the probability for Vs at each depth, given the data.</figcaption></figure>
boundary between them at a depth around 100 km, also seen here. Multiple Moho reflections (e.g the $Pp_{Ms}$ phase commonly seen in receiver functions at 14 s) which usually pollute the signal in the lower lithosphere are well handled by the inversion scheme, and not shown as a spurious interface. Note that the width of the posterior probability distribution increases with depth. This is due to the fact that rays have more separated paths with depth, and hence 3-D and moveout effects become more important resulting in data incoherence and increased uncertainty.

### 4.2.2 Joint inversion with surface waves

Although $P$ to $S$ converted phases are particularly suited to constrain the depth of discontinuities, they are only sensitive to relative changes in $S$-wave velocities in different layers, and poorly constrain absolute values. More information about the Earth structure can be brought by independent data types, such as surface wave dispersion measurements. With increasing computational power, methods to jointly invert RF and dispersion data are gaining in popularity (e.g. Ozalaybey et al. 1997; Du & Foulger 1999; Julia et al. 2000; Chang et al. 2004; Lawrence & Wiens 2004; Tkalčić et al. 2006; Yoo et al. 2007; Moorkamp et al. 2010; Salah et al. 2011; Bailey et al. 2012; González et al. 2012). When these different data types are inverted together, the complementary constraints are likely to better resolve structure. However, most of these studies are based on a regularized optimization scheme: they provide a best fitting model which depends on user-defined parameters, and quantification of uncertainties is difficult.

Here we use the transdimensional Bayesian inversion algorithm described in Bodin et al. (2012c), and jointly invert our summed seismograms with fundamental mode phase velocity measurements of Rayleigh waves (25–150 s) given by Ekström (2011). The long period dispersion measurements allow us to invert structure down to 350 km depth, and results are shown in Fig. 7. Dispersion mea-

surements allow us to get the correct velocity amplitudes, which ‘migrates’ the Moho down to 35 km. This Moho depth is very consistent with Zhou et al. (2000), Rai et al. (2006) and Julià et al. (2009), who jointly inverted the same data types at the same station. We note that the probability has been degraded in the depth range 40–80 km. In particular the discontinuity at 75 km is no longer clear. We think that this reflects some level of incompatibility between the surface waves and converted phases in this depth range which is probably due to anisotropic structure, and to our fixing the $Vp/Vs$ ratio to 1.7 in this inversion. This however, can be addressed in the future by including other types of data, such as Sp conversion data, relaxing the constraint on the $Vp/Vs$ ratio, and inverting for anisotropy by incorporating the transverse component of the $P$ waveforms (e.g. Vinnik et al. 2012).

At depths greater than 160 km, the resolution of converted phases is lower due to 3-D and moveout effects, and hence the structure is mostly constrained by Surface wave dispersion. However, a clear advantage of the class of algorithm used here is that it allows to jointly invert in a consistent manner different data sets that have different depth sensitivities, and this without having to make arbitrary choices. In particular, no explicit ‘tuning’ is needed to weigh different data-sets.

The results shown in Fig. 7 help address the debate on the thickness of the lithosphere on that part of the Indian craton. Indeed, Kumar et al. (2007) observed a consistent negative velocity jump at 110 km for Hyderabad station. They argued that this was the signature of the lithosphere–asthenosphere boundary (LAB), and concluded that such a thin continental lithosphere could explain the relatively fast velocity of the Indian plate as reconstructed from palaeomagnetic data. They subsequently inferred that the plume that partitioned Gondwanaland might also have melted the lower half of the Indian lithosphere, thus permitting faster motion. However, this model strongly disagree with the study of Mitra et al. (2006) who inverted Rayleigh wave velocity dispersions measurements and

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Joint Inversion of converted body waves and surface wave dispersion data from Ekström (2011). (a) Posterior probability distribution for $Vs$ at each depth. (b) Grey line: expected earth model. Black: maximum of the distribution at each depth. Red: 95 per cent credible interval. Right: probability of discontinuities.</figcaption></figure>
found a much thicker lithosphere, with a clear low velocity zone at 200 km. Furthermore, Oreshin *et al.* (2011) recently jointly inverted P-RFs, S-RFs and traveltime residuals, and their results also imply a thicker lithosphere under the station.

In this study, we observe the sharp negative velocity jump at 110 km that was interpreted as the LAB by Kumar *et al.* (2007). However, this is clearly not the LAB but rather a mid-lithospheric discontinuity, and the actual LAB is seen as a milder gradient between 150 and 200 km. It is important to note that the LAB and associated low velocity zone was detected by inversion of surface waves alone (Mitra *et al.* 2006), but could not be imaged with joint inversion of surface waves and RFs by Zhou *et al.* (2000), Rai *et al.* (2006) and Julià *et al.* (2009). This is because they were using standard optimization algorithms, where the unique solution obtained is not fully representative of the information contained in the data. We also note that such a model of a thick cratonic lithosphere comprising a mid-lithospheric low velocity zone was recently observed in the cratonic region of central North America by Abt *et al.* (2010), Yuan & Romanowicz (2010) and Yuan *et al.* (2011), and in the Australian craton (Ford *et al.* 2010). With the flexibility of the proposed inversion, we may look for this feature within other cratons around the world. We refer the reader to Fischer *et al.* (2010) and Rychert *et al.* (2010), for a discussion about upper-mantle negative discontinuities.

Second order discontinuities are visible in the asthenosphere, which could potentially be interpreted in terms of known global discontinuities, such as the L discontinuity (Lehmann 1961) at the bottom of the low velocity zone. Another positive discontinuity at 300 km could be associated with the X discontinuity at 300 km, which is usually observed with SS and PP precursors (Deuss 2009).

Although our analysis of the geological and geodynamical implications is rather limited, the aim here is to show that the ensemble solution produced with the cross-convolution misfit function can reconcile different observations reported in previous studies beneath the receiver.

## 5 DISCUSSION AND FUTURE WORK

Analysis of receiver functions is now a well-established seismological technique, and a number of algorithms have been proposed to infer seismic structure beneath broad-band stations. In the last decade, the intrinsic non-uniqueness of the problem has been addressed by adding surface wave dispersion data into the problem. However, a recurring issue is the definition of the misfit function one tries to minimize, and particularly the role of the data noise in this function. Here we have presented a novel norm to quantify the distance between observed and estimated measurements. This function is based on cross-convolution, and hence avoids deconvolution and the need to arbitrarily choose stabilization parameters.

This allows us to pose the problem with a proper Bayesian formalism. A posterior probability distribution is produced that fully describes our knowledge about the Earth. The variance of the posterior is directly determined by the noise on seismograms, and can be used to assess uncertainty. The posterior can be examined from several points of view to infer different properties of the model (e.g. depth of transitions, mean or most probable Vs value at one depth, etc.).

We only inverted for S-wave velocity structure while considering Vp/Vs ratio constant throughout the velocity model. An obvious improvement of the algorithm would be, albeit at increased computational cost, to also consider Vp/Vs ratios, in each layer. In addition,

layers have been assumed homogeneous and horizontal and it would be possible to treat anisotropy, slope of discontinuities, and lateral variations as unknowns in the problem. These improvements could be achieved by using densely spaced arrays, by including earthquake waveforms from a wide range of back-azimuths, and using more sophisticated forward solvers. Teleseismic P traveltime residuals could be added as in Kiselev *et al.* (2008) to obtain P-wave velocity constraints. Finally, S to P conversions (e.g. Vinnik *et al.* 2004; Geissler *et al.* 2010; Miller & Piana Agostinetti 2012; Rychert *et al.* 2012) could also be used with a similar misfit measure.

Here we have limited this study to P to S conversion, and focused on demonstrating the benefits of a cross-convolution misfit measure. This study is more a proof of concept, or study of feasibility which opens a full range of potential applications. We illustrate its power in the context of the debate on the thickness of the lithosphere under cratons.

## ACKNOWLEDGEMENTS

Thomas Bodin wishes to acknowledge support from the Miller Institute for Basic Research at the University of California, Berkeley.

## REFERENCES

* Abt, D., Fischer, K., French, S., Ford, H., Yuan, H. & Romanowicz, B., 2010. North american lithospheric discontinuity structure imaged by ps and sp receiver functions, *J. geophys. Res.*, **115**, doi:10.1029/2009JB006914.

* Audet, P., Bostock, M., Christensen, N. & Peacock, S., 2009. Seismic evidence for overpressured subducted oceanic crust and megathrust fault sealing, *Nat.*, **457**(7225), 76–78.

* Bailey, I., Miller, M., Liu, K. & Levander, A., 2012. Vs and density structure beneath the Colorado Plateau constrained by gravity anomalies and joint inversions of receiver function and phase velocity data, *J. geophys. Res.*, **117**(B2), B02313, doi:10.1029/2011JB008522.

* Bayes, T., 1763. *An Essay Towards Solving a Problem in the Doctrine of Chances*, C. Davis, Printer to the Royal Society of London.

* Bodin, T., Salmon, M., Kennett, B. & Sambridge, M., 2012a. Probabilistic surface reconstruction from multiple data sets: an example for the australian Moho, *J. geophys. Res.*, **117**(B10), doi:10.1029/2012JB009547.

* Bodin, T., Sambridge, M., Rawlinson, N. & Arroucau, P., 2012b. Transdimensional tomography with unknown data noise, *Geophys. J. Int.*, **189**, 1536–1556.

* Bodin, T., Sambridge, M., Tkalčić, H., Arroucau, P., Gallagher, K. & Rawlinson, N., 2012c. Transdimensional inversion of receiver functions and surface wave dispersion, *J. geophys. Res.*, **117**, B02301, doi:10.1029/2011JB008560.

* Bostock, M., 2002. Kirchhoff-approximate inversion of teleseismic wavefields, *Geophys. J. Int.*, **149**(3), 787–795.

* Box, G. & Tiao, G., 1973. *Bayesian Inference in Statistical Analysis*, Wiley.

* 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.

* Chang, S., Baag, C. & Langston, C., 2004. Joint analysis of teleseismic receiver functions and surface wave dispersion using the genetic algorithm, *Bull. seism. Soc. Am.*, **94**(2), 691–704.

* Darbyshire, F., Priestley, K., White, R., Stefánsson, R., Gudmundsson, G. & Jakobsdóttir, S., 2009. Crustal structure of central and Northern Iceland from analysis of teleseismic receiver functions, *Geophys. J. Int.*, **143**(1), 163–184.

* Dettmer, J., Dosso, S. & Holland, C., 2007. Uncertainty estimation in seismocoustic reflection travel time inversion, *J. acoust. Soc. Am.*, **122**, 161–176.

* Dettmer, J., Holland, C. & Dosso, S., 2009. Analyzing lateral seabed variability with bayesian inference of seabed reflection data, *J. acoust. Soc. Am.*, **126**, 56–69.
Dettmer, J., Dosso, S. & Holland, C., 2010. Trans-dimensional geoacoustic inversion, *J. acoust. Soc. Am.*, **128**, 3393–3405.

Deuss, A., 2009. Global observations of mantle discontinuities using ss and pp precursors, *Surv. Geophys.*, **30**(4), 301–326.

Di Bona, M. et al., 1998. Variance estimate in frequency-domain deconvolution for teleseismic receiver function computation, *Geophys. J. Int.*, **134**(2), 634–646.

Du, Z. & Foulger, G., 1999. The crustal structure beneath the northwest fjords, Iceland, from receiver functions and surface waves, *Geophys. J. Int.*, **139**(2), 419–432.

Duijndam, A., 1988a. Bayesian estimation in seismic inversion. Part I: principles, *Geophys. Prospect.*, **36**(8), 878–898.

Duijndam, A., 1988b. Bayesian estimation in seismic inversion. Part II: uncertainty analysis, *Geophys. Prospect.*, **36**, 899–918.

Ekström, G., 2011. A global model of Love and Rayleigh surface wave dispersion and anisotropy, 25–250 s, *Geophys. J. Int.*, **187**(3), 1668–1686.

Fischer, K., Ford, H., Abt, D. & Rychert, C., 2010. The lithosphere–asthenosphere boundary, *Ann. Rev. Earth planet. Sci.*, **38**, 551–575.

Ford, H., Fischer, K., Abt, D., Rychert, C. & Elkins-Tanton, L., 2010. The lithosphere–asthenosphere boundary and cratonic lithospheric layering beneath Australia from SP wave imaging, *Earth planet. Sci. Lett.*, **300**(3), 299–310.

Gallagher, K., Bodin, T., Sambridge, M., Weiss, D., Kylander, M. & Large, D., 2011. Inference of abrupt changes in noisy geochemical records using transdimensional change point models, *Earth planet. Sci. Lett.*, **311**, 182–194.

Geissler, W., Sodoudi, F. & Kind, R., 2010. Thickness of the Central and Eastern European lithosphere as seen by s receiver functions, *Geophys. J. Int.*, **181**(2), 604–634.

Gelman, A., Carlin, J., Stern, H. & Rubin, D., 1995. *Bayesian Data Analysis*. *Texts in Statistical Science*, Chapman & Hall.

González, O., Moreno, B., Romanelli, F. & Panza, G., 2012. Lithospheric structure below seismic stations in Cuba from the joint inversion of Rayleigh surface waves dispersion and receiver functions, *Geophys. J. Int.*, **189**, 1047–1059.

Green, P., 1995. Reversible jump MCMC computation and Bayesian model selection, *Biometrika*, **82**, 711–732.

Green, P., 2003. Trans-dimensional Markov chain Monte Carlo, *Highly Struct. Stoch. Syst.*, **27**, 179–198.

Gurrola, H., Baker, G. & Minster, J., 2007. Simultaneous time-domain deconvolution with application to the computation of receiver functions, *Geophys. J. Int.*, **120**(3), 537–543.

Hales, A., 1969. A seismic discontinuity in the lithosphere, *Earth planet. Sci. Lett.*, **7**(1), 44–46.

Haskell, N., 1953. The dispersion of surface waves on multilayered media, *Bull. seism. Soc. Am.*, **43**(1), 17–34.

Helffrich, G., 2006. Extended-time multitaper frequency domain cross-correlation receiver-function estimation, *Bull. seism. Soc. Am.*, **96**(1), 344–347.

Hopcroft, P., Gallagher, K. & Pain, C., 2009. A Bayesian partition modelling approach to resolve spatial variability in climate records from borehole temperature inversion, *Geophys. J. Int.*, **178**, 651–666.

Iaffaldano, G., Bodin, T. & Sambridge, M., 2012. Reconstructing plate-motion changes in the presence of finite-rotations noise, *Nat. Commun.*, **3**, 1048, doi:10.1038/ncomms2051.

Iaffaldano, G., Bodin, T. & Sambridge, M., 2013. Slow-downs and speed-ups of India–Eurasia convergence since: data-noise, uncertainties and dynamic implications, *Earth planet. Sci. Lett.*, **367**, 146–156.

Julia, J., Ammon, C., Herrmann, R. & Correig, A., 2000. Joint inversion of receiver function and surface wave dispersion observations, *Geophys. J. Int.*, **143**(1), 99–112.

Julià, J., Jagadeesh, S., Rai, S. & Owens, T., 2009. Deep crustal structure of the Indian shield from joint inversion of P wave receiver functions and Rayleigh wave group velocities: implications for Precambrian crustal evolution, *J. geophys. Res.*, **114**(B10), B10313, doi:10.1029/2008JB006261.

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.

Kind, R., Yuan, X. & Kumar, P., 2012. Seismic receiver functions and the lithosphere–asthenosphere boundary, *Tectonophysics*, **536**, 25–43.

Kiselev, S., Vinnik, L., Oreshin, S., Gupta, S., Rai, S., Singh, A., Kumar, M. & Mohan, G., 2008. Lithosphere of the Dharwar craton by joint inversion of P and S receiver functions, *Geophys. J. Int.*, **173**(3), 1106–1118.

Kosarev, G., Kind, R., Sobolev, S., Yuan, X., Hanka, W. & Oreshin, S., 1999. Seismic evidence for a detached Indian lithospheric mantle beneath Tibet, *Science*, **283**(5406), 1306–1309.

Kumar, P., Yuan, X., Kumar, M., Kind, R., Li, X. & Chadha, R., 2007. The rapid drift of the indian tectonic plate, *Nat.*, **449**(7164), 894–897.

Kumar, P., Kind, R. & Yuan, X., 2010. Receiver function summation without deconvolution, *Geophys. J. Int.*, **180**(3), 1223–1230.

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

Lawrence, J. & Wiens, D., 2004. Combined receiver-function and surface wave phase-velocity inversion using a niching genetic algorithm: application to Patagonia, *Bull. seism. Soc. Am.*, **94**(3), 977–987.

Leahy, G., Saltzer, R. & Schmedes, J., 2012. Imaging the shallow crust with teleseismic receiver functions, *Geophys. J. Int.*, **191**, 627–636.

Lehmann, I., 1961. S and the structure of the upper mantle, *Geophys. J. R. astr. Soc.*, **4**, 124–138.

Levander, A., Schmandt, B., Miller, M., Liu, K., Karlstrom, K., Crow, R., Lee, C. & Humphreys, E., 2011. Continuing colorado plateau uplift by delamination-style convective lithospheric downwelling, *Nat.*, **472**(7344), 461–465.

Li, A., Fischer, K., van der Lee, S. & Wysession, M., 2002. Crust and upper mantle discontinuity structure beneath eastern north america, *J. geophys. Res.*, **107**, doi:10.1029/2002JB001891.

Ligorria, J. & Ammon, C., 1999. Iterative deconvolution and receiver-function estimation, *Bull. seism. Soc. Am.*, **89**(5), 1395–1400.

Liu, K., Levander, A., Zhai, Y., Porritt, R. & Allen, R., 2012. Asthenospheric flow and lithospheric evolution near the mendocino triple junction, *Earth planet. Sci. Lett.*, **323**, 60–71.

Luo, X., 2010. Constraining the shape of a gravity anomalous body using reversible jump Markov chain Monte Carlo, *Geophys. J. Int.*, **180**(3), 1067–1079.

MacKay, D., 2003. *Information Theory, Inference, and Learning Algorithms*, Cambridge Univ. Press.

Malinverno, A., 2002. Parsimonious Bayesian Markov chain Monte Carlo inversion in a nonlinear geophysical problem, *Geophys. J. Int.*, **151**(3), 675–688.

Menke, W. & Levin, V., 2003. The cross-convolution method for interpreting SKS splitting observations, with application to one and two-layer anisotropic earth models, *Geophys. J. Int.*, **154**(2), 379–392.

Miller, M. & Piana Agostinetti, N., 2012. Insights into the evolution of the italian lithospheric structure from s receiver function analysis, *Earth planet. Sci. Lett.*, **345**, 49–59.

Minsley, B., 2011. A trans-dimensional Bayesian Markov Chain Monte Carlo algorithm for model assessment using frequency-domain electromagnetic data, *Geophys. J. Int.*, **187**(1), 252–272.

Mitra, S., Priestley, K., Gaur, V. & Rai, S., 2006. Shear-wave structure of the south indian lithosphere from Rayleigh wave phase-velocity measurements, *Bull. seism. Soc. Am.*, **96**(4A), 1551–1559.

Moorkamp, M., Jones, A. & Fishwick, S., 2010. Joint inversion of receiver functions, surface wave dispersion, and magnetotelluric data, *J. geophys. Res.*, **115**(B4), B04318, doi:10.1029/2009JB006369.

Nikulin, A., Levin, V. & Park, J., 2009. Receiver function study of the Cascadia megathrust: evidence for localized serpentinization, *Geochem. Geophys. Geosyst.*, **10**(7), Q07004, doi:10.1029/2009GC002376.

Oldenburg, D., 1981. A comprehensive solution to the linear deconvolution problem, *Geophys. J. R. astr. Soc.*, **65**(2), 331–357.

Oreshin, S., Vinnik, L., Kiselev, S., Rai, S., Prakasam, K. & Treussov, A., 2011. Deep seismic structure of the Indian shield, Western Himalaya, Ladakh and Tibet, *Earth planet. Sci. Lett.*, **307**(3), 415–429.

Ozalaybey, S., Savage, M., Sheehan, A., Louie, J. & Brune, J., 1997. Shear-wave velocity structure in the northern Basin and Range province from the combined analysis of receiver functions and surface waves, *Bull. seism. Soc. Am.*, **87**(1), 183–199.
Park, J. & Levin, V., 2000. Receiver functions from multiple-taper spectral correlation estimates, *Bull. seism. Soc. Am.*, **90**(6), 1507–1520.

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

Piana Agostinetti, N. & Malinverno, A., 2010. Receiver function inversion by trans-dimensional Monte Carlo sampling, *Geophys. J. Int.*, **181**(2), 858–872.

Rai, S., Priestley, K., Gaur, V., Mitra, S., Singh, M. & Searle, M., 2006. Configuration of the indian Moho beneath the NW Himalaya and Ladakh, *Geophys. Res. Lett.*, **33**(15), L15308.

Ray, A. & Key, K., 2012. Bayesian inversion of marine CSEM data with a trans-dimensional self parametrizing algorithm, *Geophys. J. Int.*, **191**, 1135–1151.

Revenaugh, J. & Jordan, T., 1991. Mantle layering from scs reverberations 3. The upper mantle, *J. geophys. Res.*, **96**(B12), 19811–19824.

Romanowicz, B. et al., 1991. The geoscope program: present status and perspectives, *Bull. seism. Soc. Am.*, **81**(1), 243–264.

Rychert, C., Shearer, P. & Fischer, K., 2010. Scattered wave imaging of the lithosphere–asthenosphere boundary, *Lithos*, **120**(1), 173–185.

Rychert, C. et al., 2012. Volcanism in the afar rift sustained by decompression melting with minimal plume influence, *Nat. Geosci.*, **5**(6), 406–409.

Rychert, C.A., Laske, G., Harmon, N. & Shearer, P.M., 2013. Seismic imaging of melt in a displaced hawaiian plume, *Nat. Geosci.*, **6**, 657–660.

Salah, M., Chang, S. & Fonseca, J., 2011. Crustal structure beneath the lower Tagus Valley, Southwestern Iberia using joint analysis of teleseismic receiver functions and surface-wave dispersion, *Geophys. J. Int.*, **184**(2), 919–933.

Sambridge, M., 1999. Geophysical inversion with a neighbourhood algorithm I. Searching a parameter space, *Geophys. J. Int.*, **138**(2), 479–494.

Sambridge, M., Bodin, T., Gallagher, K. & Tkalčić, H., 2013. Transdimensional inference in the geosciences, *Phil. Trans. R. Soc. A*, **371**, 20110547, doi:10.1098/rsta.2011.0547.

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.

Sarkar, D., Kumar, M., Saul, J., Kind, R., Raju, P., Chadha, R. & Shukla, A., 2003. A receiver function perspective of the Dharwar craton (India) crustal structure, *Geophys. J. Int.*, **154**(1), 205–211.

Saul, J., Kumar, M.R. & Sarkar, D., 2000. Lithospheric and upper mantle structure of the Indian shield, from teleseismic receiver functions, *Geo-phys. Res. Lett.*, **27**(16), 2357–2360.

Shearer, P., 1991. Imaging global body wave phases by stacking long-period seismograms, *J. geophys. Res.*, **96**(B12), 20353–20364.

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.

Shibutani, T., Sambridge, M. & Kennett, B., 1996. Genetic algorithm inversion for receiver functions with application to the crust and uppermost mantle structure beneath eastern Australia, *Geophys. Res. Lett*, **23**, 1829–1832.

Sivia, D., 1996. *Data Analysis: A Bayesian Tutorial*, Oxford Univ. Press.

Smith, A., 1991. Bayesian computational methods, *Phil. Trans.: Phys. Sci. Eng.*, **337**(1647), 369–386.

Srinivas, D., Srinagesh, D., Chadha, R. & Kumar, M. R., 2013. Sedimentary thickness variations in the Indo-Gangetic foredeep from inversion of receiver functions, *Bull. seism. Soc. Am.*, **103**(4), 2257–2265.

Stipčević, J., Tkalčić, H., Herak, M., Markušić, S. & Herak, D., 2011. Crustal and uppermost mantle structure beneath the external Dinarides, Croatia, determined from teleseismic receiver functions, *Geophys. J. Int.*, **185**(3), 1103–1119.

Tarantola, A. & Valette, B., 1982. Inverse problems = quest for information, *J. Geophys.*, **50**(3), 150–170.

Tauzin, B., Debayle, E. & Wittlinger, G., 2010. Seismic evidence for a global low-velocity layer within the Earth/'s upper mantle, *Nat. Geosci.*, **3**(10), 718–721.

Thomson, W., 1950. Transmission of elastic waves through a stratified solid medium, *J. Appl. Phys.*, **21**, 89–93.

Tkalčić, H., Pasyanos, M., Rodgers, A., Gok, R., Walter, W. & Al-Amri, A., 2006. A multistep approach for joint modeling of surface wave dispersion and teleseismic receiver functions: Implications for lithospheric structure of the Arabian Peninsula, *J. geophys. Res.*, **111**(B11), B11311, doi:10.1029/2005JB004130.

Tkalčić, H., Young, M., Bodin, T., Ngo, S. & Sambridge, M., 2013. The shuffling rotation of the Earth/'s inner core revealed by earthquake doublets, *Nat. Geosci.*, **6**, 497–502.

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

Vinnik, L., Reigber, C., Aleshin, I., Kosarev, G., Kaban, M., Oreshin, S. & Roecker, S., 2004. Receiver function tomography of the central Tien Shan, *Earth planet. Sci. Lett.*, **225**(1), 131–146.

Vinnik, L., Aleshin, I., Kaban, M., Kiselev, S., Kosarev, G., Oreshin, S. & Reigber, C., 2006. Crust and mantle of the Tien Shan from data of the receiver function tomography, *Izvest. Phys. Solid Earth*, **42**(8), 639–651.

Vinnik, L., Kiselev, S., Weber, M., Oreshin, S. & Makeyeva, L., 2012. Frozen and active seismic anisotropy beneath southern Africa, *Geophys. Res. Lett.*, **39**(8), doi:10.1029/2012GL051326.

Yoo, H., Herrmann, R., Cho, K. & Lee, K., 2007. Imaging the three-dimensional crust of the Korean Peninsula by joint inversion of surface-wave dispersion and teleseismic receiver functions, *Bull. seism. Soc. Am.*, **97**(3), 1002–1011.

Young, M. K., Rawlinson, N. & Bodin, T., 2013. Transdimensional inversion of ambient seismic noise for 3D shear velocity structure of the tasmanian crust, *Geophysics*, **78**(3), WB49–WB62.

Yuan, H. & Romanowicz, B., 2010. Lithospheric layering in the North American craton, *Nat.*, **466**(7310), 1063–1068.

Yuan, H., Dueker, K. & Schutt, D. L., 2008. Testing five of the simplest upper mantle anisotropic velocity parameterizations using teleseismic S and SKS data from the billings, Montana Passcal array, *J. geophys. Res.-Solid Earth* (1978–2012), **113**(B3), doi:10.1029/2007JB005092.

Yuan, H., Romanowicz, B., Fischer, K. & Abt, D., 2011. 3-D shear wave radially and azimuthally anisotropic velocity model of the North American upper mantle, *Geophys. J. Int.*, **184**(3), 1237–1260.

Zhou, L., Chen, W. & Özalaybey, S., 2000. Seismic properties of the central Indian Shield, *Bull. seism. Soc. Am.*, **90**(5), 1295–1304.