# Receiver functions from seismic interferometry: a practical guide

Benoit Tauzin, Thanh-Son Pham, Hrjoje Tkalcic

► To cite this version:

Benoit Tauzin, Thanh-Son Pham, Hrjoje Tkalcic. Receiver functions from seismic interferometry: a practical guide. Geophysical Journal International, 2019, 217 (1), pp.1-24. 10.1093/gji/ggz002. hal-02293275

HAL Id: hal-02293275

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

Submitted on 31 Aug 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>

# Receiver functions from seismic interferometry: a practical guide

Benoit Tauzin,<fnref n="1" />,<fnref n="2" /> Thanh-Son Pham<fnref n="2" /> and Hrvoje Tkalčić<fnref n="2" />

Accepted 2019 January 6. Received 2018 May 5; in original form 2018 December 21

## SUMMARY

This paper reviews the concepts underlying the well-documented receiver functions (RFs) method, and places it in the conceptual framework of seismic interferometry. We first propose a simple and efficient approach for isolating the receiver-side seismic response (i.e. the record of reflections and conversions in the stratification beneath receivers): this method makes use of the *P*-wave coda recorded on the radial and vertical components of three-component stations, applies spectral whitening, which is followed by auto- and cross-correlation. The interferometric principle underpinning RFs analysis is shown theoretically and illustrated in practice using earthquake records and synthetic waveforms computed from simple structures. We point out to a major limitation, which is the contamination of the receiver-side response by propagation effects in the source-side structure. We then apply our approach to teleseismic earthquake data recorded in California. We show that the reconstructed vertical and horizontal seismic responses can be back-projected to illuminate the crustal and mantle structure. We build comparable ~300-km-long seismic reflectivity profiles from pure *P*-wave reverberations and from the converted wavefield across the forearc and arc of the southern Cascadia subduction zone. Then, we show a case of processing data from narrow bandpass, short-period and single-component sensors, usually unsuitable for RFs analysis. Finally, through the same interferometric principle, we attempt to demonstrate a link between event- and noise-based seismic interferometry. We demonstrate that it is possible to extract approximate responses from the records of low-magnitude—down to 4.5—teleseismic earthquakes. We make a comparison of these estimates with the result from the autocorrelation of the continuous ambient noise seismic wavefield. While the amplitudes of the extracted receiver-side responses are mutually different, their phases are in a relative agreement. This development opens a way to the use of small magnitude teleseismic earthquakes to characterize the receiver-side structure.

**Key words:** Crustal imaging; Seismic interferometry; Seismic noise.

## 1 INTRODUCTION

Teleseismic receiver functions (RFs) represent the seismic response near a recording station to an incoming teleseismic plane wave. The major part of the incoming energy is directly recorded by the surface sensors whereas a significant part of it is reverberated in the form of multiple-order reflections within the underlying stratification. Due to the partitioning of energy at seismic boundaries, these later arriving seismic phases in the coda project onto both the vertical and horizontal components of the three-component stations. A standard processing method for extracting RFs is to calculate the deconvolution of one component by another (Clayton & Wiggins 1976; Langston 1979; Ammon 1991), the underlying hypothesis being that the incoming wavefield (the source) and the medium response can be separated in projection on the different components (Vinnik 1977; Langston 1979). The deconvolution acts as a matched filter,

correlating a source signal, with a complex signal, the seismogram, to detect the repetition of the source in the seismogram (e.g. Ligorria & Ammon 1999).

Claerbout (1968) showed that the reflection response of a horizontally layered acoustic medium can be retrieved from the autocorrelation of a plane-wave transmission response (Fig. 1a). This theory is a component of a more general theory of seismic interferometry that has been developed for arbitrarily heterogeneous and elastic media (e.g. Wapenaar 2003). Building on this theory and pointing to the similarities between the operations of cross-correlation, convolution or deconvolution and summation of the recorded wavefields, Galetti & Curtis (2012) proposed that seismic interferometry could be seen as a generalization of RF analysis. Sun & Kennett (2016) exploited stacked correlograms of three-component records of teleseismic events at the Australian stations. They found a good match between the recovered *P* and *S* reflection responses, and the

<footnote n="1">¹Université de Lyon, Université Claude Bernard Lyon 1, ENS, CNRS, Laboratoire de Géologie de Lyon: Terre, Planètes, Environnement, 69622 Villeurbanne, France. E-mail: benoit.tauzin@univ-lyon1.fr</footnote>

<footnote n="2">²Research School of Earth Sciences, The Australian National University, Canberra, Australian Capital Territory 2601, Australia</footnote>
(a) Acoustic case

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

(b) P-S, P-P-S, P-S-S modes

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

(c) P-P-P mode

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

**Figure 1.** Generalized receiver functions (RFs) from seismic interferometry and the different modes of scattering. (a) Principle of *P*-wave coda seismic interferometry as originally proposed by Claerbout (1968) in the acoustic case. The near-vertical incidence *P* wave and all of its reverberations in the subsurface at left give the transmission response recorded at the seismometer. The reflection response at right is obtained from the autocorrelation of the transmission response. (b) The key crustal phases recorded on teleseismic radial RFs. These phases are associated with three modes, *P*-S, which corresponds to the direct conversion from *P* into SV, and P-P-S, P-S-S (S-P-S), which are reverberations in the structure below the station with mixed legs in *P* and S. Based on Crotwell *et al.* (1999), we note Pms, PPms and PSms the interactions of the modes with the Moho discontinuity (see text). (c) The crustal phase recorded on the vertical component and of attention in this paper is PPmp, a first-order *P*-wave reverberation within the stratification associated with the mode P-P-P. This mode can be extracted from *P*-coda auto-correlation. At right, we show that the seismic interferometry principle from Claerbout (1968) applies to extract the different modes. The reflection at the free surface acts as a virtual source, and because of this, the interferometric principle is sometimes referred as virtual deep seismic sounding (Yu *et al.* 2012, 2013).

reflectivity obtained from *P*-to-*S* and *S*-to-*P* RFs. Pham & Tkalčić (2017) devised a new processing technique and applied it to the coda of transmitted *P* waves to infer the solid rock-ice interface in Antarctica and the Moho in South Africa. Their method is based on spectral whitening, followed by auto- and cross-correlation of the radial and vertical components of the teleseismic records. The recovery of the vertical and radial responses of the receiver-side subsurface is exceptionally good. The above studies can be seen as a confirmation of the idea that seismic interferometry is a generalization of RF analysis (Galetti & Curtis 2012).

Schuster *et al.* (2004) extended the derivations of the Earth’s elastic response through seismic interferometry to include an imaging step after the correlations. They termed this method ‘interferometric seismic imaging’. The method of interferometric imaging was applied to several case studies in global seismology with the use of passive earthquakes. For examples, Abe *et al.* (2007) extracted surface-related multiple reflections from cross-correlations of teleseismic waves among pairs of receiver stations on a profile along the Itoigawa-Shizuoka Tectonic Line in Japan. They found the Moho and lower crustal reflectors consistent with RF imaging. Tseng *et al.*
(2009), using a method they referred to as ‘virtual deep seismic sounding’, reconstructed the P-to-P reflection response originating from the conversion of teleseismic shear-waves into a P wave at the free surface (S-P-P) below the Hi-CLIMB array in Tibet. They successfully imaged the Moho along the profile. Ruigrok & Wapenaar (2012) used teleseismic P waves travelling through the Earth’s core to obtain a quasi zero-offset reflection response below this same Hi-CLIMB experiment. The processing gave good results in recovering the Moho, and possibly illuminating lithospheric fault zones. Yu et al. (2012) applied virtual deep seismic sounding to a dense seismic array in the North China Craton and recovered the signature of the Moho and of an intracrustal discontinuity. More recently, Sun & Kennett (2017) exploited stacked event-based correlograms at another dense array in the western and central North China Craton and potentially inferred the signature of a mid-lithospheric discontinuity.

An alternative way to extract the body-wave reflectivity beneath seismic recorders is to correlate continuous seismic ambient noise records in the high-frequency band (0.5–4 Hz; e.g. Poli et al. 2012; Gorbatov et al. 2013). The underlying assumption is that a component of the seismic ambient noise is associated with a wavefield propagating at near-vertical incidence, assumption used in the original seismic interferometry principle from Claerbout (1968). Poli et al. (2012), Tibuleac & von Seggern (2012), Gorbatov et al. (2013), Kennett (2015), Kennett et al. (2015), and more recently, Oren & Nowack (2017) and Becker & Knapmeyer-Endrun (2018), made use of this principle to extract the reflectivity of the crust and mantle underneath seismic arrays of variable nature (broad-band, short-period) and geometry (linear, two-dimensions, regularly and heterogeneously distributed stations). The origin of this vertically propagating wavefield in the ambient noise remains however debated (e.g. Bonnefoy-Claudet et al. 2006; Ruigrok et al. 2011; Poli et al. 2012; Lin & Tsai 2013).

Noise is the generic term used to denote ambient vibrations of the ground caused by sources such as tide, water waves striking the coast, turbulent wind, effects of wind on trees or buildings, industrial machinery, cars and trains, etc. (Bonnefoy-Claudet et al. 2006). It has been shown that P waves can be extracted from noise records at arrays far from offshore storms in small time-windows and in the 0.1–1.0 Hz frequency band (e.g. Roux et al. 2005; Ruigrok et al. 2011). However, the comprehensive review of Bonnefoy-Claudet et al. (2006) suggests that a major contribution to the noise are surface sources prone to generate horizontally propagating surface waves. Any wavefield propagation can be seen as the superimposition of frequency and slowness components (e.g. Kennett & Pham 2018), and from theoretical and empirical analysis of seismic noise interferometry, we know that a favorable slowness and backazimuth distribution is required to extract a meaningful response from the noise (Claerbout 1968; Wapenaar 2003; Paul et al. 2005; Galetti & Curtis 2012). The imperfect knowledge of such components can make the extraction of the subsurface response unreliable. This doubt leads for instance to questions about including or not the contribution of large earthquakes in noise correlation studies for reconstructing the global Earth’s response (Poli et al. 2012; Lin & Tsai 2013).

This article reviews the concepts underlying the well-known RF method, and points to the analogy with seismic interferometry. We demonstrate this analogy in theory in Section 2, and in practice in Section 3. We use the processing outlined by Pham & Tkalčić (2017), and the record of a good-quality earthquake at a broad-band three-component station. We also illustrate the processing on synthetic waveforms computed from simple 1-D structures, pointing to

the major limitation. In Section 4, we show that the reconstructed vertical (and horizontal) seismic responses can be back-projected to illuminate the crustal and mantle structure beneath a dense array of seismic stations (the ‘interferometric seismic imaging’ principle from Schuster et al. 2004). Then in Section 5, we apply these principles to the records of teleseismic earthquakes at a dense broad-band temporary array in northern California, the Flexible Array Mendocino Experiment (FAME). We image the gently subducting Gorda Plate beneath the Cascadia arc north of the Mendocino Triple Junction. We find quite comparable images using the pure P-wave mode of reverberations, and the more conventional P-to-S converted modes of vibrations. In Section 6, we discuss other potential applications. We show a case of processing data from narrow bandpass, short-period and single-component sensors, usually unsuitable for RFs analysis. We envision a promising application on long-duration short-period single-component seismic deployments in California. Finally, using the same principle (Claerbout 1968), we attempt to demonstrate a link between event- and noise-based seismic interferometry. We show that it is possible to extract an approximate receiver-side response from the records of low-magnitude (4.5–5.4) teleseismic earthquakes. We compare these results with the result obtained from processing continuous seismic ambient noise, such as in the studies of Gorbatov et al. (2013), Kennett (2015), Kennett et al. (2015), Becker & Knapmeyer-Endrun (2018) and others. We find that the phase in the approximate response is in a relative agreement, but the amplitude is not, probably due to differences in the seismic wavefield directivity. Because the RFs represent approximations for both the horizontal and vertical receiver-side responses, we refer to them as ‘generalized RFs’ (Galetti & Curtis 2012). This generalization applies to any type of seismic sources, including the seismic ambient noise.

# 2 THEORETICAL BACKGROUND

## 2.1 Receiver-side seismic response

Because at teleseismic distances the spherical wavefront has a large curvature, the seismic excitation $S(t, p)$ of the receiver-side structure is well-described with a plane-wave with slowness $p$. The seismograms corresponding to a teleseismic P plane-wave and its coda recorded as functions of time $t$ on the vertical and radial components $Z(t, p)$ and $R(t, p)$ result from the convolution (*) of this ‘source function’ $S(t, p)$, with the instrumental response $I(t)$, and $E_z(t, p)$ and $E_r(t, p)$, the vertical and radial seismic response of the stratification beneath the receiver, respectively.

$$ 
\begin{aligned}
Z(t, p) &= S(t, p) * I(t) * E_z(t, p) \\
R(t, p) &= S(t, p) * I(t) * E_r(t, p).
\end{aligned}
\qquad (1) $$

In this formulation, $E_z$ and $E_r$ terms are receiver terms, that is ‘RFs’ (Langston 1979). They do not represent the actual Green’s function from the source to the receiver (Wapenaar 2003; Galetti & Curtis 2012), as the effects of propagation away from the receivers are included in the source function.

Consider a layer of thickness $H$ over a half space. The seismic responses $E_{z,r}(t, p)$ of the medium are the sequential records of phases, caused by the original P plane-wave, that are reflected and converted within the layer,

$$ E_{z,r}(t, p) = \sum_{m,n}^{\infty} a_{z,r}^{m,n} \delta(t - t_p - m\tau_p - n\tau_s). \qquad (2) $$
In eq. (2), $t_p$ is the observed arrival time of the direct P wave at the station, $\tau_p$ and $\tau_s$ are the times taken by the transmitted P and converted S wave fronts, both retaining the original ray parameter $p$, to sweep the layer vertically, and $m, n$ represent the number of times that the P and S wave fronts translate over the layer.

The vertical sweeping times are related to the vertical slownesses $\eta_p$ and $\eta_s$ in the layer, and the thickness $H$,

$$ 
\begin{aligned}
\tau_p &= H\eta_p = H(1/v_p^2 - p^2)^{1/2} \\
\tau_s &= H\eta_s = H(1/v_s^2 - p^2)^{1/2}.
\end{aligned} 
\quad (3) $$

The superscript indices $m$ and $n$ satisfy that $m+n$ is non-negative even, $m \ge 0$, and $n \ge -1$. For instance, when $m = n = 0$, the term $a^{0,0}\delta(t - t_p)$ represents the peak of the direct P phase. When $m = -1$ and $n = 1$, $\delta[t - t_p - (\tau_s - \tau_p)]$ is the arrival associated with the direct P-S conversion, and when $m = 2$ and $n = 0$, the phase is a P-wave reverberated one time in the layer (Fig. 1c).

$a_{z,r}^{m,n}$ are the amplitudes of the converted and reflected phases in projection on the vertical and radial components. Without considering anelastic attenuation, these amplitudes result from the multiple interactions of the phases with the interfaces in the stratification,

$$ a^{m,n} = \prod_{i}^{m+n+1} a_{i}^{m,n}. \quad (4) $$

Each interaction $i$ affects the seismic amplitudes with a proportion $a_i^{m,n}$ that depends on the contrast in elastic properties (density and velocity) across the interface and the incidence angle of the wave, namely the transmission and reflection coefficients (Section 2.6). Few interfaces such as the free surface and the Moho have sufficient contrast to generate large coefficients, and even for these interfaces only a limited number of interactions are likely to be observed ($m + n \le 4$).

In the frequency domain, eqs (1) and (2) express as

$$ 
\begin{aligned}
Z(\omega, p) &= S(\omega, p)E_z(\omega, p) \\
R(\omega, p) &= S(\omega, p)E_r(\omega, p),
\end{aligned} 
\quad (5) $$

and

$$ E_{r,z}(\omega, p) = e^{-i\omega(t-t_p)} \sum_{m,n}^{\infty} a_{z,r}^{m,n} e^{i\omega(m\tau_p+n\tau_p)}, \quad (6) $$

where $\omega$ is the angular frequency, and the instrument response $I(\omega)$ is dropped for the sake of simplicity.

## 2.2 Receiver functions

The radial receiver function, $RF_r(\omega, p)$, in frequency domain is defined as

$$ RF_r(\omega, p) = \frac{R(\omega, p)}{Z(\omega, p)}. \quad (7) $$

This division in the frequency domain is known as a deconvolution operation. We can multiply the numerator and denominator in eq. (7) by the complex conjugate of the denominator and

$$ RF_r(\omega, p) = \frac{R(\omega, p)Z^*(\omega, p)}{Z(\omega, p)Z^*(\omega, p)}. \quad (8) $$

Here, as pointed out by Galetti & Curtis (2012) (also e.g. Ammon 1991), a connection exists between the deconvolution and cross-correlation operations. The cross-correlation of two time-series is defined in the frequency domain as the multiplication of one with the conjugate of the other. Therefore, the numerator of the second

term in eq. (8) is the cross-correlation of the radial and vertical components.

If we substitute eq. (5) into eq. (8),

$$ 
\begin{aligned}
RF_r(\omega, p) &= \frac{S^2}{||Z||^2} E_r(\omega, p) E_z^*(\omega, p) \\
&= \frac{1}{||E_z||^2} E_r(\omega, p) E_z^*(\omega, p).
\end{aligned} 
\quad (9) $$

This formulation does not make assumption on the nature of the source $S(\omega, p)$ or the vertical response $E_z(\omega, p)$, in opposition with the classical RF method (Section 2.4). The cross-correlation term $E_r E_z^*$ includes the phase information in the RF. The term $||E_z||^2$, which is the autocorrelation of the vertical response, has a null-phase.

## 2.3 Virtual sources at the free-surface

Substituting eq. (6) into the cross-correlation term in (9)

$$ 
\begin{aligned}
E_r(\omega, p) E_z^*(\omega, p) &= \sum_{m,n} a_r^{m,n} e^{i\omega(m\tau_p+n\tau_s)} \\
&\quad \times \sum_{p,q} a_z^{p,q} e^{-i\omega(p\tau_p+q\tau_s)}
\end{aligned} 
\quad (10) $$

where the common absolute delay term, $e^{-i\omega(t-t_p)}$, cancels out. The phases extracted through correlation of the seismograms have a delay relative to the direct P wave.

This result underlies the principle of teleseismic interferometry (e.g. Schuster et al. 2004; Tseng et al. 2009; Ruigrok & Wapenaar 2012; Galetti & Curtis 2012). Because all the phases are obtained with respect to the direct P-wave arrival, the P-wave reflection under the free surface provides a virtual source for subsequent reflections off the Moho or any other interface in the crust and mantle (Fig. 1, right). Because of this virtual source, the method consisting in extracting through correlation the reflections response is also sometimes referred as ‘virtual deep seismic sounding’ (Yu et al. 2012, 2013).

After some algebra, we can write eq. (10) in a simplified form,

$$ E_r(\omega, p) E_z^*(\omega, p) = \sum_{\Phi} a_{\Phi} e^{i\omega\tau_{\Phi}}, \quad (11) $$

where $a_\Phi$ and $\tau_\Phi$ are the amplitudes and delay times of peaks associated with summed products of pairs of transmitted and reflected phases $\Phi$ on the seismograms. The coefficients $a_\Phi$ are the products of the amplitudes of these phases on the radial and vertical components.

## 2.4 Classical RFs

Given the small incidence angle of the teleseismic wavefield (typically below 30°), the direct P-wave ($m = 0, n = 0$ in eq. 6) is the dominant phase recorded on the vertical component $Z(\omega, p)$. In the original version of the RF method, Langston (1979) made the assumption that the vertical seismic response can be neglected

$$ 
\begin{aligned}
E_z(\omega, p) &\approx a_z^{0,0} e^{-i\omega(t-t_p)} \quad \text{or} \\
E_z(t, p) &\approx a_z^{0,0} \delta(t - t_p).
\end{aligned} 
\quad (12) $$

This assumption considers that the first and higher order reverberations of the P wave ($n = 0, m = 2$ or more in eq. 6) are not recorded on the vertical component, and that the vertical component can be used as an approximation for the source function, $S(\omega, p) \approx Z(\omega, p)$.
p). In this case, the radial RF from eq. (9) is the shifted version of the radial response of the structure,

$$RF_r(\omega, p) \approx \frac{1}{a_z^{0,0}} E_r(\omega, p) e^{i\omega(t-t_p)}, \quad (13)$$

where the term $e^{i\omega(t-t_p)}$ cancels the absolute delay $t_p$ in the radial component response.

## 2.5 Reflected and converted phases

Given the rapid decrease of amplitude through multiplication of reflection coefficients in the layering (eq. 4), the observed reverberations in eq. (6) have a small number of P or S legs in the stratification (typically $m+n \le 4$). Only the first (and possibly second) order reverberations are well observed on the RFs.

In the following, we will adopt the phase naming convention from Crotwell et al. (1999) where P represents a compressional wave, upgoing or downgoing, in the crust or mantle, and the lower case p is strictly an upgoing P wave within the crust or mantle. We make a similar use of capital and lower case letters for shear waves, S or s. The interaction with the Moho discontinuity is denoted with the letter ‘m’. This interaction can be generalized to any interface of the medium. For instance, we use ‘t’ and ‘b’ for interactions with the top and bottom of a low-velocity zone in a simplified model of the Cascadia forearc (Section 3).

Following eq. (11), the pairs of phases contributing to the phases $\Phi$ observed on the cross-correlation have respective indices $m' = m - p$ for P-legs, and $n' = n - q$ for S-legs. With this convention, we denote the direct conversion for ($m' = -1, n' = 1$) as Pms, the reverberation for ($m' = 1, n' = 1$) as PPms, and the reverberation for ($m' = 0, n' = 2$) as PSms. The P-P-P mode interacting with the Moho discontinuity for ($m' = 2, n' = 0$) is noted PPmp (Fig. 1c, left).

The traveltimes of these main phases are given by

$$\tau_{[\Phi=Pms]} = \tau_s - \tau_p \quad (m'=-1, n'=1), \quad (14)$$

$$\tau_{[\Phi=P Pms]} = \tau_s + \tau_p \quad (m' = 1, n' = 1),$$

$$\tau_{[\Phi=P Sms]} = 2\tau_s \quad (m' = 0, n' = 2),$$

$$\tau_{[\Phi=P Pmp]} = 2\tau_p \quad (m' = 2, n' = 0),$$

where the expressions of $\tau_p$ and $\tau_s$ are given in eq. (3). Generalized to a stratified velocity model, these equations provide a way to convert the time-domain converted and reflected signals to the depth domain, and an imaging principle for back-projecting the correlations into the spatial domain and reconstructing an image of the subsurface (Section 4).

## 2.6 Seismic amplitudes

The amplitude $a_{z,r}^{m,n}$ of converted and reflected phases brings complementary information on the seismic structure below the receivers. These amplitudes are function of the transmission and reflection coefficients ($a_i^{m,n}$ in eq. 4), themselves function of the contrasts in velocity ($\Delta v_s, \Delta v_p$) and density ($\Delta \rho$) across the interfaces, as well as the incidence angle (slowness) of the wave impinging on the stratification. For first-order seismic interfaces, analytical expressions of these transmission and reflection coefficients are derived for instance by Aki & Richards (1980, pp. 144–145), Lay & Wallace (1995, p. 98) and Chapman (2004, p. 218), with different notations, normalization and sometimes mistakes in their expressions.

In Fig. 2, we show the result of computing the solid-solid trans-

mission and reflection coefficients, as well as the free-surface reflection coefficients for the pure P-P and S-S modes, and the P-S converted modes. We use the formula from Aki & Richards (1980) for a first order discontinuity separating a one-layer crust (medium 1) from a half-space mantle (medium 2). The contrast in velocity and density at the Moho boundary is downward positive. We indicate with $2 \to 1, 1 \to 2$ and $1 \to air$ the directions of propagation of seismic waves from a medium to another, and the subscripts p, s the wave-type in each medium.

For slownesses between 4.6 and 8.8 s deg$^{-1}$, the P wave ($2 \to 1$) and reverberations in the uppermost layer ($1 \to 2$) impinges on the Moho with pre-critical angles of incidence, 20–40° for the P wave and 10–20° for the S wave (grey bands in Fig. 2). Over the range of slownesses, the transmission and reflection coefficients have different sensitivities to the elastic and density contrasts at the Moho interface. For instance, the pure P-wave reflection $R_{pp}^{2 \to 1}$ atop the Moho discontinuity contributing to the phase PPmp (with amplitude $a_{[\Phi=P Pmp]} \propto a^{2,0}$ on the RF), depends mainly on the contrast in P-wave elastic impedance (the product $\rho v_p$), and due to the partitioning of energy in SV, also on the shear-wave elastic impedance. In contrast, the transmission coefficient ($T_{ps}^{2 \to 1}$) for a converted phase contributing to the phase Pms (with amplitude $a_{[\Phi=Pms]} \propto a^{-1,1}$) depends mostly on $\Delta v_s$.

Fig. 2(a) shows that most of the amplitude of the upward propagating P wave that interacts with the Moho is transmitted as a P wave, and only $T_{ps}^{2 \to 1} \approx 15$ per cent of it is converted to a shear wave. Fig. 2(d) shows that the reflection at the free surface is quite efficient ($R_{pp,ps}^{1 \to air} \approx 75–100$ per cent), and changes the polarity of the pure P-wave reflected pulse. The reflections on top of the Moho discontinuity are less efficient ($R_{pp,ps,ss}^{1 \to 2} < 20$ per cent; Figs 2b and c), preserving the polarity for P-P, but changing the sign of the pulses for P-S and S-S modes of reflection.

The amplitude of the phases on the RFs scales with the product of these coefficients (eq. 4). With the same notation as above,

$$a_{[\Phi=Pms]} \propto T_{ps}^{2 \to 1} a_0,$$

$$a_{[\Phi=P Pms]} \propto R_{ps}^{1 \to 2} R_{pp}^{1 \to air} T_{pp}^{2 \to 1} a_0,$$

$$a_{[\Phi=P Sms]} \propto R_{ss}^{1 \to 2} R_{ps}^{1 \to air} T_{pp}^{2 \to 1} a_0,$$

$$a_{[\Phi=P Pmp]} \propto R_{pp}^{1 \to 2} R_{pp}^{1 \to air} T_{pp}^{2 \to 1} a_0,$$

where $a_0$ is the initial amplitude of the P incident wave in medium 2. The multiplication of coefficients over the incidence angles given by the grey bands in Fig. 2 results in Pms and PPms amplitudes that are positive on the radial component. In contrary, the PSms and PPmp phases have a negative polarity on this component. On the vertical component, the amplitudes of all these phases is negative. Due to symmetry relations in the reflection and transmission coefficients, the case of a downward negative velocity/impedance contrast gives opposite signs of coefficients. The amplitudes of the Pms and PPms phases tend also to decrease with decreasing incidence angle (Fig. 2), whereas the amplitudes of the PSms and PPmp phases increase. This is of importance when considering waves impinging on dipping interfaces (e.g. Cassidy 1992).

When the width of velocity gradients gets large and their shape more complex, then the use of reflection/transmission coefficients for first order interfaces is not anymore a good approximation, except for long-period seismic waves. We illustrate this in Fig. 3. The response can be reconstructed by subdividing the velocity gradient in thin homogeneous layers and then computing the effective response from plane-wave matrix algorithms such as designed by Kennett (1983) or Chapman (2004, 2003). Following Ma et al.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Variations of the reflection and transmission coefficients with the angle of incidence of the incoming wavefield. The coefficient $a/a_0$ plotted along the Y-axis correspond to the amplitude ratio of the transmitted/reflected waves with respect to the incident wave. We consider the configuration of a one-layer crust (medium 1) of elastic properties $v_{p1} = 6.5 \text{ km s}^{-1}$, $v_{s1} = 3.75 \text{ km s}^{-1}$, $\rho_1 = 2.85 \text{ g cm}^{-3}$ sandwiched between the air and a mantle half-space (medium 2) of properties $v_{p2} = 8.2 \text{ km s}^{-1}$, $v_{s2} = 4.55 \text{ km s}^{-1}$, $\rho_2 = 3.39 \text{ g cm}^{-3}$. We represent the coefficients for the cases of (a) upward (2 → 1) P-P and P-S transmissions, (b and c) downward (1 → 2) P-P, P-S, and S-S reflections and (d) the upward P-P and P-S reflections at the free surface (1 → air).</figcaption></figure>

(2012), we modified the algorithm from Chapman (2003) to calculate the generalized reflection-transmission coefficients from any stack of homogeneous layers. We show in Fig. 3 the comparison of the plane-wave frequency responses $T_{ps}^{2→1}$, $R_{ps}^{1→2}$ and $R_{pp}^{1→2}$, obtained for a first-order discontinuity at the base of the crust and a more gradual boundary with 10 km thickness (Fig. 3a). A gradual boundary acts as a low-pass filter for the coefficients (Figs 3b–d). As expected, the responses at frequencies below 0.03 Hz are similar to the response obtained for the first order discontinuity, meaning that the analysis of the signal at these frequencies does not permit to distinguish the broad gradient from the first order discontinuity. At higher frequencies the transmission and reflection coefficients quickly decrease in amplitude and show a typical low-amplitude ringing. This behavior is at the origin of the frequency dependence

of the visibility of deep mantle seismic discontinuities (e.g. Helffrich & Wood 1996). An interesting observation in this context is the different frequency cut-off of the responses of broad gradients depending on the considered mode (dotted red lines in Fig. 3). This feature has been demonstrated through numerical experiments by Bostock (1999). Transmitted conversions are less sensitive to transition width than reflected modes, with the first minimum in the transmission response occurring at roughly twice the characteristic frequency for P-P and P-S reflections. The numerical tests from Bostock (1999) showed that an appropriate rule of thumb for the detection of a broad velocity gradient is that the thickness of the gradient must be smaller than $\Delta z \approx \lambda/4$ for pure-mode reflections, and $\Delta z \approx \lambda_P/2$ for conversions, $\lambda$ being the mode wavelength. In
<figure><img src="image_7.png" /><figcaption>**Figure 3.** The low-pass effect of a 10-km-thick velocity gradient on the reflection and transmission coefficients of the different modes for a plane-wave of slowness $p = 0.05$ s km$^{-1}$. (a) The $v_p$ profiles used for computing the effective responses in transmission and reflection. The $v_s$ and $\rho$ structures are scaled to the $v_p$ structure. The first-order gradient structure is similar to the one used for computing the coefficients in Fig. 2. The gradual structure is composed of 10-km-thick homogeneous layers at the top and the base of the gradient, the latter being subdivided in 18 homogeneous layers with 0.5 km thickness. (b) Comparison between the upward P-S transmission responses for the 1st order interface (continuous line) and for the gradual structure (line highlighted with circles). (c) Comparison for the downward P-P reflection. (d) Comparison for the downward P-S reflection. We indicate in red the frequency of the first-minimum in coefficient amplitude.</figcaption></figure>

consequence, the joint use of the different modes, including the frequency dependence of their amplitudes, provides complementary constraints on $v_p$, $v_s$ and $\rho$ beneath the stations.

## 3 PROCESSING

The source function $S(t, p)$ is defined here as the superposition of the effects of variable source time function, frequency content, radiation pattern of the earthquakes, and wave-propagation away from the receiver. Extracting from the coda coherent information on the structure across large-aperture seismic networks and with different sources is often not possible unless we account for or equalize the effect of the source, that is apply the deconvolution operation in eqs (7) and (8).

### 3.1 Radial receiver functions ($RF_r$): removing source characteristics through water-level deconvolution

In practice, this division is often unstable because holes exist in the spectrum of the vertical component $Z(\omega)$, and some sort of regularization needs to be introduced to avoid the division by small or zero values. From eqs (8) and (9),

$$ 
\begin{aligned}
RF_r(\omega, p) &= \frac{R(\omega, p) Z^*(\omega, p)}{Z'(\omega, p)} G(\omega), \\
&= \frac{S^2}{Z'} E_r(\omega, p) E_z^*(\omega, p) G(\omega),
\end{aligned} 
\quad (15) 
$$

where $Z'(\omega, p)$ is the regularized vertical component, and $G(\omega) = \exp(-\frac{\omega}{4a^2})$ is a gaussian filter. The gaussian filter can be applied at any time on the data, with the parameter $a$ controlling the width of the gaussian.

The regularization can be applied for instance using a water-level parameter, $c$, as proposed by for example Clayton & Wiggins (1976), where:

$$ Z'(\omega, p) = \max \left( \|Z(\omega, p)\|^2, c \max_{\forall \omega} \|Z(\omega, p)\|^2 \right). \quad (16) $$

An example of the application of the water-level deconvolution is shown in Fig. 4 for one teleseismic event recorded at the station ME25 from the FAME Mendocino experiment in northern California. In practice, this operation flattens the spectrum in the bandpass of the Gaussian filter (Fig. 4b).

### 3.2 Vertical receiver functions ($RF_z$): the case of water-level deconvolution

If the above approach is efficient to recover an approximation of the radial response $R/Z \approx E_r$, as shown in Fig. 4(c), it does not allow to recover the vertical response $E_z$ by computing the vertical RF, $RF_z = Z/Z'$. The vertical RF is

$$ RF_z(\omega, p) = \frac{Z(\omega, p) Z^*(\omega, p)}{Z'(\omega, p)} G(\omega). \quad (17) $$

The vertical response is included into the deconvolved source function $Z'(\omega)$ (contrary to the assumption in eq. 12). The division of the
<figure><img src="image_8.png" /><figcaption>(a) spectral domain deconvolution</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) spectral whitening</figcaption></figure>

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

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

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

**Figure 4.** Similarities between spectral domain deconvolution and spectral whitening for computing receiver functions. (a) Vertical $Z(t)$ and radial $R(t)$ components for the September 9th 2008 teleseismic event recorded at station ME25 from the FAME Mendocino experiment in northern California. The frame in red corresponds to the time window of interest for imaging the crust and uppermost mantle structure. The waveforms have been low pass-filtered below 1 Hz with a similar filter as $G(\omega)$ in eqs (11) and (13). (b) Illustration of the water-level deconvolution in the spectral domain with (left) the original spectrum of the radial component $R(\omega)$ (in black) and the regularized spectrum of the vertical component $Z'(\omega)$ (in red), and (right) the result of their division in the spectral domain (the vertical and radial RFs). The spectrum $Z'(\omega)$ is used as an approximation for the earthquake source function and the structure away from the receiver side. (c) Result of the deconvolution in the time domain, with at the bottom the radial RF and at the top the vertical RF. As expected for vertical RF, because we divided $Z$ by a scaled version of itself ($Z'$), we do not recover the vertical impulse response. (d) Illustration of the spectral whitening with (left) the original spectrum of the vertical component $Z(\omega)$ (in black) and its regularized version $Z''(\omega)$ (in red), and (right) the result of the division of $Z$ by $Z''$ (the vertical RF). $Z''(\omega)$ represents a better approximation for the earthquake source function. The resulting transfer function is the approximate vertical response of the structure in the frequency domain. (e) Result of spectral whitening and correlations in the time domain. The auto-correlation at the top is an approximation for the vertical response $E_z(t)$, and the cross-correlation at the bottom is an approximation of the radial response $E_r(t)$. The autocorrelation is symmetric with respect to the zero-lag whereas the cross-correlation is not.

vertical component by a slightly modified version of itself results in the time domain in a single Dirac impulse (or Gaussian) function $\delta(t)$ (Fig. 4c). P-wave reverberations associated with horizontal interfaces arrive at the surface at the same angle of incidence as the direct P wave and are therefore removed from the RFs by source equalization (Langston 1979).

### 3.3 Vertical receiver functions ($\hat{R}F_z$) through spectral whitening and autocorrelation

The method introduced by Pham & Tkalčić (2017) to extract the PPmp phases from teleseismic P-wave records shares similarities with the spectral domain deconvolution described above. Similarly to eq. (17), the vertical response of the medium is synthesized with

$$ \hat{R}F_z(\omega, p) = \frac{Z(\omega, p) Z^*(\omega, p)}{Z''(\omega, p)} G(\omega). \quad (18) $$

The difference with eq. (17) arises from the choice of the regularized source function, $Z''(\omega, p)$, that is a smoothed version of the spectrum of the vertical component (Pham & Tkalčić 2017; Fig. 4d)

$$ Z''(\omega, p) = \frac{1}{2W} \int_{\omega-W}^{\omega+W} \|Z(\nu, p)\|^2 d\nu \quad (19) $$

in the continuous domain. In the discrete domain, each individual frequency sample $Z''_n$ of this smoothed spectrum is obtained with

$$ Z''_n = \frac{1}{2N+1} \sum_{\nu=n-N}^{n+N} \|Z_\nu\|^2, \quad (20) $$

where $Z_\nu$ are the individual frequency samples of the initial complex spectrum of the vertical component. If $\Delta\omega$ is the frequency sampling of the spectrum, $W = 2N\Delta\omega$ is the spectral whitening
width (in Hertz) over which the spectrum is smoothed. This formulation follows the running-absolute-mean normalization method in the time domain (Bensen et al. 2007; Pham & Tkalčić 2017).

With this regularized source spectrum, eq. (18) is called a spectral whitening (or balancing) operation because, as illustrated in Fig. 4(d), it amplifies the high-frequency content and reduces at the same time the low-frequency content of the teleseismic records (Pham & Tkalčić 2017). If we omit the complex conjugate $Z^*$ at the numerator of eq. (18) (the correlation), the amplitudes in the spectrum are equalized, but the phase in the original seismogram is preserved. This phase includes, in particular, the source mechanism (radiation pattern) of the earthquake. The multiplication of $Z$ by its conjugate in eq. (18) cancels out the phase and provides the vertical reflectivity of the structure $E_z$ (Fig. 4e, top). In comparison with the result of water-level deconvolution in Fig. 4(c), the trace obtained through spectral whitening and autocorrelation shows multiple peaks that can now be interpreted in terms of P-wave reflectivity (Claerbout 1968).

## 3.4 Radial receiver functions ($\hat{R}F_r$) through P-wave coda cross-correlation

Because spectral whitening and autocorrelation allow to extract the P-wave reflectivity (Pham & Tkalčić 2017), we can conclude that this processing is more efficient than classical deconvolution schemes to equalize the contribution from the source in the spectra of seismic records, and thereby isolate the vertical response of the structure. What about extracting the classical P-to-SV radial RFs using correlation? This approach has been discussed by Galetti & Curtis (2012) and implemented by Sun & Kennett (2016), and shown to give good results for the records of teleseismic earthquakes at stations located in Australia.

We compute here the radial RF by analogy with the water-level deconvolution with a simple spectral division of the radial component $R$ by the smoothed vertical component $Z''$:

$$ \hat{R}F_r(\omega, p) = \frac{R(\omega, p) Z^*(\omega, p)}{Z''(\omega, p)} G(\omega) \quad (21) $$

This operation gives a RF in agreement with the one obtained from water-level deconvolution (Figs 4c, e and 5a).

Among standard methods exists also the iterative time-domain deconvolution (Kikuchi & Kanamori 1982; Ligorria & Ammon 1999). This approach is based on the least-square minimization of the difference between the observed radial component and a predicted signal generated by the convolution of an iteratively updated spike-train with the vertical component (Ligorria & Ammon 1999). We show for the event recorded at the station ME25 in Fig. 5(a) a comparison between the radial RF obtained from iterative time-domain deconvolution (in blue), and the RFs obtained through regularization with water-level and spectral whitening (red and black). The three waveforms are quite similar, although it appears that the RFs obtained from iterative time-domain and water-level deconvolutions are the most similar to each other, especially between 15 and 20 s. This can be explained by two reasons.

First, if classical deconvolution approaches (Langston 1979; Ammon 1991; Ligorria & Ammon 1999) allow to recover an approximation of the radial seismic response $E_r$, they do not provide the P-wave reflection response $E_z$, even in projection on the radial component. To illustrate this, we show in Fig. 5(b) the result of applying these three deconvolution methods on the synthetics computed from a simple Earth's model composed of a flat horizontal crustal layer

on top of a homogeneous half-space. Radial RFs obtained from iterative time-domain and water-level deconvolutions (blue and red) are identical, but the RF obtained from whitening (black) shows, additionally, the PPmp arrival, which is the P-wave reflection in the layer (Fig. 1c). Therefore compared to classical deconvolution methods, spectral whitening allows to recover the projection of the P-wave energy onto the radial component. Similarly, we find that the vertical autocorrelation (Fig. 5b) records the projection of converted modes from P into SV (Pms, PPms and PSms) in addition to the P-wave reflection response (PPmp). An obvious improvement would be to turn the recording components along the main axis of vibration in L, Q, T or in P, SV, SH (e.g. Reading et al. 2003; Rondenay 2009). A second explanation for the observed differences in the RFs in Fig. 5(a) is the approximation for the far-field spectrum of the source function.

## 3.5 Approximations in the far-field spectrum of the source function

We sketch in Fig. 6 the sequence of linear filters that are applied to the initial source time function through the wave propagation from the source to the receiver. As written in eq. (1) and (5), the seismograms result in the time domain from the convolution of the source with the Earth's seismic response, and correspond to a multiplication of these two transfer functions in the frequency domain. If the event can be approximated by a simple point source, its far-field spectrum corresponds to a low-pass (Fig. 6, top) with a roll-off (the rate at which attenuation increases beyond the cut-off frequency) as a negative power of $\omega$, and the cut-off depending on the earthquake magnitude (e.g. Houston & Kanamori 1986; Brune 1970). Accounting for anelasticity, the wavefield is attenuated along its path. Its spectrum is therefore progressively shifted toward lower frequency, but remains flat over its corresponding bandpass (Fig. 6, middle). In contrast, the Earth modulates the seismic response by adding spikes and holes in the far-field spectrum of the teleseismic wavefield. The body-wave reflectivity has a 'sparse' spectrum, with spikes corresponding to the first and multiple-order reverberations recorded with increasing time in the seismic record (eq. 11). Under some conditions, the smoothed spectrum of the vertical component (in blue, bottom left in Fig. 6) can be a good approximation for the source function, and dividing the original component spectrum by this smoothed version allows to extract the response of the structure (e.g. Fig. 4d, right).

However, the response of the structure near the hypocenter is also sparse and is, at the receiver side, superimposed on the source function representing the wavefield incoming beneath the stratification. The spectrally whitened spectrum of the vertical component may not fully incorporate this response, and the deconvolution using this approximate source function may bias the recovery of the vertical and radial RFs. We investigated, by means of a simple synthetic test, the effect of the response of the structure near the source not being properly accounted for. The result is shown in Fig. 7. First, we designed a model with a low-velocity zone (LVZ) sandwiched in between a one-layer crust and a homogenous mantle half-space (Fig. 7a). This model is an approximation for the structure expected in the Cascadia forearc (e.g. Li & Nábělek 1999; Rondenay et al. 2001; Tauzin et al. 2017). Second, we built synthetic vertical and radial component seismograms with the same distribution of slownesses as recorded at the station ME25 of the FAME experiment (Fig. 7c). The responses are computed with a reflectivity algorithm
<figure><img src="image_13.png" /><figcaption>**Figure 5.** Comparison between deconvolution methods for the observation of a single teleseismic earthquake, and for synthetic data predicted by a simple structure. (a) Result of computing the RFs from the records of the 2008 September 9 teleseismic event shown in Fig. 4. We make a comparison between the classical water-level deconvolution (e.g. Clayton & Wiggins 1976; Langston 1979) (trace in red; eq. 11), iterative time-domain deconvolution (Ligorria & Ammon 1999) (trace in blue), and the spectral whitening (trace in black; eq. 15) methods. The three RFs are superimposed at the top. (b) A comparison of RFs (blue, black and red traces from the bottom to the top) and vertical and radial autocorrelations (black traces at the top) for synthetic data computed in a simple, flat, horizontal crustal layer (thickness $H = 35$ km, $v_p = 6.00$ km s⁻¹, $v_s = 3.5$ km s⁻¹, $\rho = 2690$ kg m⁻³), on top of a homogeneous half-space ($v_p = 8.0$ km s⁻¹, $v_s = 4.6$ km s⁻¹, density $\rho = 3330$ kg m⁻³). In addition, the vertical and radial responses (shown in red at the top) are computed using a reflectivity algorithm (Kennett 1983; Randall 1989) for a plane-wave of slowness $p = 0.07$ s km⁻¹.</figcaption></figure>

(Randall 1989). The source functions, convolved with these responses, are approximated by the observed vertical components at the station, and are therefore not sparse in the frequency domain. The resulting 365 synthetic vertical and radial component seismograms are then deconvolved using spectral whitening and compared with those obtained from water-level deconvolution (Fig. 7c).

The vertical responses (Fig. 7b) show a simple sequence of arrivals, the main ones being related to the P-wave reverberations at the top and the base of the LVZ, noted here PPtp and PPbp. The radial responses are more complicated because the top and the base of the LVZ give direct conversions Pts and Pbs, as well as two pairs of reverberations PPts, PPbs and PSts, PSbs (see e.g. Audet et al. 2009). The result of water-level deconvolution (Fig. 7c, left) is a good approximation for the theoretical radial response except near 10 s where the PPtp and PPbp waveforms are not recovered (explanation in previous section). In contrast, the average response obtained from spectral whitening (Fig. 7c, middle) shows a better recovery of the PPtp and PPbp pulses on the radial RF. However, individual RFs obtained from spectral whitening, after cross- and autocorrelations (Fig. 7c, middle and right), are contaminated by incoherent noise that blurs out the signal of major converted and reflected arrivals. This noise is the result of incompletely accounting for the sparse component in the source functions. This component is extreme here, because the observed vertical components at ME25 (that we have used as an approximation for the source functions) include, in addition to the source structure, the structure at the receiver side and a component of ambient seismic noise. This general behavior is sometimes observed on individual autocorrelations (Section 5). This does not mean that the autocorrelations cannot be used, because a simple linear stack is sufficient to extract the coherent structure in the data (Fig. 7c, top).

Finally, we conclude this section by referring to a simple strategy to tackle the problem of source contamination in P-wave coda autocorrelations. These strategies have been investigated before, notably by Li & Nábělek (1999), Hansen et al. (2013), and probably many others. They mainly consist in finding a smoothed average source function for each teleseismic event across large-aperture seismic arrays. This approach is restrictive in the sense that every single earthquake needs to be recorded at most of the seismic stations, which limits its use to a single or synchronous experiments. A ‘brute force’ approach such as employed here is efficient because stacking in the time (or depth) domain leads to destructive interference of the near-source structure contributions.

# 4 IMAGING PRINCIPLE

## 4.1 Time-to-depth conversion

Using the generalization of eq. (15) to 1-D stratified velocity models, it is possible to determine the delays between the transmitted P wave, and the phases in the coda (the P-S, P-P-S, P-S-S and P-P-P modes) for any depth of interaction in the model. There is no need for a physical interface in the model. A time-to-depth relationship is then constructed for each mode, allowing to convert the single time-domain seismic trace into four different depth traces associated each with a different mode.

In Fig. 8, we show after stacking the effect of applying the depth conversion to the time-domain radial and vertical synthetic seismic responses of Fig. 7. After this operation, the pulses associated with the top (Pts, PPTS, PSts and PPtp) and bottom (Pbs, PPbs, PSbs and PPbp) of the low-velocity zone all fall at depths of 25 and 35 km, as specified in the model. To get similar amplitudes, we flipped the
<figure><img src="image_14.png" /><figcaption>**Figure 6.** A sequence of linear filters is applied to the source (top left) through wave-propagation in the Earth (middle), which results in the seismogram recorded at the seismometer (bottom left). For simplicity, we omit the instrument response. The seismogram is the multiplication in the frequency domain of the source with the transfer functions from the elastic and anelastic structure. The response of the structure near the source, along the path away from the receiver, and beneath the receiver is ‘sparse’. The bottom left panel shows superimposed the initial spectrum of the source in red, the spectrum after attenuation through the Earth in thick black, the spectrum recorded at the seismometer in thin black, and the spectrally whitened spectrum in blue. The spectrally whitened spectrum of the vertical component only partially incorporates the response of the structure near the source. This biases the recovery of the responses for the receiver side (see Fig. 7).</figcaption></figure>

polarity of the traces associated with the P-S-S and P-P-P modes. In addition, the depth-conversion shrinks the signal at later time more significantly than at early time. To obtain a similar wavelength in the signals of the different modes, the time-to-depth conversion is applied on data filtered with different low-pass corners, at 1 s for the P-S mode, 3.5 s for the P-P-P mode and 5 s for the P-P-S and P-S-S modes.

## 4.2 Common reflection point stacking for vertical autocorrelations

We include an imaging step after the correlations to benefit from the dense coverage of 1-D or 2-D seismic arrays, and synthesize spatial variations of the Earth’s seismic response (e.g. Schuster *et al. 2004*). Our approach is based on a simple and efficient algorithm of back-projection used for stacking by common conversion point (CCP) single-mode (Zhu 2000; Wittlinger *et al. 2004*) and multimode P-to-SV RFs (Hetényi 2007; Tauzin *et al. 2016, 2017*; Audet & Schaeffer 2018). The present work extends the derivation of this imaging

principle to the vertical RFs, that is the P-P-P back-scattered mode obtained through P-wave coda autocorrelation. The amplitude is back-projected to the last reflection point of the P-wave in the stratification (see Fig. 1), so we refer to this operation as ‘common reflection point stacking’. In a subsequent paper, we will focus on the application of this imaging principle to the data from dense arrays in the Cascadia forearc, Mendocino Triple Junctions system, and San Andreas fault zone.

We discretize the medium with a stack of layers and, through 3-D ray tracing in this medium, compute the coordinates of the P- and S-rays beneath the array. Then, the seismic traces are time-to-depth converted using a 1-D or 3-D velocity model. Using the ray-coordinates, the amplitude of each depth sample on the correlations/RFs is back-projected (migrated) and stacked to the medium to obtain an image from each different mode (Wittlinger *et al. 2004*; Tauzin *et al. 2016*).

We illustrate such an imaging principle in Fig. 9 with synthetic data computed from the simple seismic model in Fig. 8(b), and the acquisition geometry associated with the Cascadia93 experiment in Oregon (Tauzin *et al. 2017*). The linear Cascadia93 seismic array
<figure><img src="image_15.png" /><figcaption>**Figure 7.** A simple synthetic experiment showing the (negative) impact of incorporating a 'sparse' component in the source functions. (a) We use a 1-D model of a low-velocity zone sandwiched between a one-layer crust and a half-space mantle. (b) We computed the vertical $E_z(t)$ and radial $E_r(t)$ responses for plane-waves incoming on the stratification with a distribution of slownesses similar to observed at the station ME25 of the FAME experiment. To compute synthetic vertical $Z(t)$ and radial $R(t)$ seismograms (eq. 1), we used as source functions the observed waveforms on the vertical components of station ME25. (c) We applied a water-level deconvolution (left; eqs 11 and 12) and spectral whitening followed by cross- and auto-correlation (middle and right; eqs 13 and 15). The recovered individual traces (bottom, in black) are sorted by order of increasing slowness and superimposed on the theoretical responses (traces in red). At the top, we plot the linear stack of all traces (black) and of the theoretical responses (red).</figcaption></figure>

in Oregon and theoretical geographical location of PPmp reflection points at 35 km depth are shown in Fig. 9(a). We show also the medium discretization, using a regular parametrization in the three spatial dimensions x, y and z (Fig. 9b). We define x as the direction

of the profile. With such a linear array, we only consider the 2-D variations of the structure along the x- and z-direction. We therefore elongate the parameterization in the y-direction ($\Delta y \to \infty$), and all the information obtained from the back-projection in the 3-D space
<figure><img src="image_16.png" /><figcaption>**Figure 8.** Effect of time-to-depth conversion on the theoretical vertical and radial responses computed from a low-velocity zone (LVZ) model. (a) Time-domain signals. The components are filtered with different low-pass filters (1, 3.5 and 5 s periods) depending on the mode used for imaging. We indicate with blue and red lines the seismic phases related to the top and bottom of the LVZ. (b) Traces obtained from the four different modes after the depth conversion. The traces have a similar vertical wavelength. For each mode, the interaction with the top and bottom of the LVZ is mapped to its correct depth. We have flipped in polarity the traces associated with the P-S-S and P-P-P modes. At the right-hand side, we represent the vertical shear-wave velocity used to compute the synthetics. The top and bottom of the LVZ are marked with blue and red lines.</figcaption></figure>

is projected on to the vertical plane along the seismic array (see Wittlinger *et al.* 2004, for the first application). An advantage of this approach is that by increasing the distance of projection $Δy$, we can artificially reduce aliasing effects due to a poorer ray coverage (i.e. when using sparse 2-D seismic arrays), and therefore improve the quality of the imaging if the structure is void of significant 3-D heterogeneity.

In Fig. 9(c), we show the result of the back-projection of the four different modes for the synthetic data associated with the flat LVZ model. In this example, the discretization has cells of $Δx = 2$ km and $Δz = 1$ km in the horizontal and vertical directions. If the seismic interfaces associated with the LVZ are mapped at their correct depths, each individual image for a given mode presents artifacts related to the other modes. For instance, the image obtained from the P-P-P mode is strongly contaminated by the P-P-S mode (PPbs in Fig. 9c), and more weakly by the P-S mode (Pts and Pbs on the vertical component in Fig. 8a). For this depth range, the P-P-S and P-S-S modes show the higher level of contamination. In our previous work (Tauzin *et al.* 2016), we showed how it is possible to partially mitigate this contamination problem by building a single ‘multi-mode’ image through phase-weighted stack (Schimmel & Paulssen 1997). We apply such an approach in the application Section 5. We note also that various approaches have been proposed to filter out the different contributions from the different modes

directly on the raw data/RFs (e.g. Guan & Niu 2017; Ainiwaer & Gurrola 2017).

Compared to the other modes, a benefit of applying an imaging principle to the P-P-P mode is the increase in lateral sampling of the structure beneath the array. We plot in Fig. 10 the lateral distance of sampling of the rays for the P-P-P (continuous lines) and P-S, P-P-S and P-S-S modes (dashed lines). The distance is two times larger for the P-P-P mode, and it increases with increasing slowness (decreasing epicentral distance) and with the depth of investigation. At 400 km depth for a slowness of 0.04 s km⁻¹, the lateral extents would be ~75 and 145 km, respectively for the P*-S and P-P-P modes. These extents further increase to ~155 and 350 km for a slowness of 0.08 s km⁻¹.

## 5 APPLICATION TO A DENSE BROADBAND SEISMIC NETWORK

We use broad-band records of teleseismic earthquakes at stations from a network located in northern California, the Flexible Array Mendocino Experiment (FAME). This 2-D temporary deployment with ~25 km interstation spacing has been used by Liu *et al.* (2012) to illuminate the lithospheric structure of the Gorda Plate and the Mendocino Triple Junction system (inset in Fig. 11). We
<figure><img src="image_17.png" /><figcaption>**Figure 9.** Interferometric seismic imaging from different modes of the seismic wavefield. (a) Locations of the Cascadia93 linear array (black triangles) and of PPmp reflection points at 35 km depth (red dots) beneath the array. The background color at depth describes the low-velocity zone (LVZ) model. The frame in black delimits the domain used for imaging the structure beneath the array. A small rectangular zone near 122°E longitude and 35 km depth is enlarged in (b). (b) Enlargement over the parameterization. The parameterization is elongated along the direction orthogonal to the profile, that is $\Delta y \gg \Delta x$ and $\Delta y \gg \Delta z$. Red dots correspond to the piercing points of the *P* waves at the depth of 35 km. (c) Result of the back-projection of different modes obtained from the synthetic data for the LVZ model. The acquisition geometry associated with the observed data at the Cascadia93 array is used. Seismic rays (black lines) for the *P*-wave (mode P-P-P) and *S* wave (modes P-S, P-P-S and P-S-S) are shown for a selected portion of the profile.</figcaption></figure>

use the original records of teleseismic earthquakes from our previous RF studies (Tauzin et al. 2016, 2017). These earthquakes have magnitudes larger than 5.5, and are located at epicentral distances between 30 and 90°. The waveforms have been preliminary filtered in the 0.013–1 Hz band. We then select the best-quality waveforms based on a signal-to-noise ratio criterion. We compute

this signal-to-noise ratio (snr) from the root-mean-square (rms) in a 100 s-wide window after the P-arrival for the signal, in a 100 s window before the P-arrival for the noise, setting the ratio threshold at 2.0. The records are rotated along the vertical and radial directions using the theoretical backazimuth from the station to the epicentre.
<figure><img src="image_18.png" /><figcaption>**Figure 10.** Lateral sampling around the station for different depths of investigation and the different modes (P-P-P, solid; P-S, P-P-S and P-S-S, dashed line). The curves are plotted for slownesses of 0.04 (blue) and 0.08 s km⁻¹ (red), corresponding to the bounding epicentral distances for the observation of the converted waves.</figcaption></figure>

## 5.1 Station ME25

We illustrate in Fig. 11 the processing on the records of a single seismic station from the FAME experiment, ME25. In Fig. 11(a), the waveforms are the broadband vertical components aligned on the theoretical P-wave arrival-time predicted by the earth model iasp91 (Kennett & Engdahl 1991), and sorted by increasing slowness. Due to the variable frequency content, source-time functions and radiation patterns, the P-wave codas are not directly comparable from one earthquake to another. In Fig. 11(b), we apply spectral whitening but do not account for the complex conjugate $Z^*$ at the numerator of eq. (18) (the correlation). The fast Fourier transform is computed here with a frequency step of 0.015 Hz. We compute the running average of the vertical component spectrum in a normalization window of $W = 0.15$ Hz, with $2N + 1 = 11$ frequency samples (eq. 20). We then divided the original spectrum by this running average to flatten out the spectrum and amplify the high frequencies in the coda. The result in the time domain shows that the phase is preserved, recognizable by the similar polarity of the direct P-arrival in the original and spectrally whitened records (Figs 11a and b).

The multiplication of $Z$ by its conjugate in eq. (18) cancels out the phase and provides the vertical reflectivity of the structure $E_z$ (Fig. 11c). The autocorrelations are computed in the time-domain and normalized by the amplitude of the peak at zero-lag. The maximum amplitude on these vertical RFs corresponds to the maximum of the correlation at zero-lag. Due to the coefficients of reflection $R_{pp}^{1 \to air}$, the free-surface reflection of the incoming teleseismic P-wave changes its polarity (see Section 2.6). The reflection $R_{pp}^{1 \to 2}$ on top of a layer with a higher elastic impedance ($\rho v_p$) conserves its polarity, so the overall P-P reflection on a discontinuity such as the Moho has a negative polarity (Section 2.6). Using eq. (15), we emphasize at this station two potential arrivals with amplitudes and moveout consistent with reflections from positive seismic discontinuities located at about 20 and 30 km depths (red lines in Fig. 11c).

In addition to these P-wave reflections observed at ~6.5 and 9.0 s on the vertical component RFs, we should also observe on the radial component two triplets of converted waves associated with the interaction with the interfaces at 20 and 30 km depths. The direct P-S conversions should be recorded at ~2.5 and 3.6 s (eq. 15). The P-P-S and P-S-S reverberations should be recorded at ~9.0, 12.5 and 11.5, 16.5 s, respectively. We show in Fig. S1 in supplementary material a similar section as in Fig. 11(c), but for the radial component.

Because of interference effects, the identification of individual P-P-S and P-S-S reverberations is quite difficult. However, the direct conversions at the two discontinuities are quite remarkable near ~3 s after the P wave. In Figs 12(a) and (b), the result of time-to-depth conversion indicates that the recovery of the sequence of conversions/reflections in the stratification below the station ME25 is consistent between the vertical and radial components, further validating our processing and inference on the structure.

In Fig. 11(d), we plot the vertical RFs as a function of backazimuth. As the structure is not perfectly horizontal beneath the station (dipping slab; see inset in Fig. 11), significant azimuthal variations in the waveforms is to expect, in particular when comparing the western and eastern backazimuths. The arrivals associated with the 20 and 30-km-deep interfaces are partially disrupted in the eastern backazimuths. These are potentially effects of the dipping structure. There are also distinct, sloping events appearing, for instance at times 3–10 s for backazimuths 29–127°, and time ~10 s for backazimuths 228–235° (red arrows). A detailed analysis of these features in the context of the dipping structure beneath the ME25 stations is out of the scope of this study. An important assumption in this seismic interferometry/RF principle is the presence of a planar horizontal stratification beneath the recording seismometer (Fig. 1). If this assumption seems valid at first order to explain the arrivals associated with the discontinuities at 20 and 30 km depth, it is likely that complex wave-propagation effects related to the 2-D/3-D structure are also observed on individual RFs.

## 5.2 Structure beneath the FAME Mendocino experiment

Despite the complexity in lithospheric structure under the FAME experiment, we show in Figs 12(c) and (d) that the simple back-projection of the radial and vertical RFs using a 1-D earth model gives a good first-order idea of the general architecture of the subduction zone. The seismic section in Fig. 12(c) is the result of combining the three converted modes in a single image through phase-weighted stack (Schimmel & Paulssen 1997). This broadband structure is identical to described in the study from Tauzin et al. (2017), and very similar to the one observed for the subduction of the Juan de Fuca Plate further north in Oregon, northwestern Washington, and beneath Vancouver (Bostock et al. 2002; Abers et al. 2009; Audet et al. 2009; Audet & Schaeffer 2018). In the western portion, the multimode profile shows a dipping low-velocity layer (LVZ) with a red interface at its base and a blue interface at its top, corresponding either to the Gorda oceanic crust and/or a residual/accreted layer of sediments (Rondenay et al. 2001). The continental Moho is observed in the eastern portion, dipping westward from ~30 to 40 km beneath the volcanic arc. The continental Moho then gets shallower near 122°W. A positive converter sits on top of the LVZ near 22 km depth at 122.5°W longitude (Tauzin et al. 2017).

P-wave reverberations recorded on broadband vertical RFs provide comparable information on the structure (Fig. 12d). The subhorizontal North America continental Moho and the dipping Gorda oceanic Moho are clearly identified in red at similar depths and with similar dip angles. However, the oceanic Moho seems doubled. A possibility is that the uppermost interface corresponds to the true reflection PPmp atop the Moho, and the deepest horizon is associated with the PPms phase, which is also observed in projection on synthetic vertical autocorrelations for simple structures (e.g. Fig. 5b). However, our synthetic experiment in Section 4.2 showed that the PPms phase (PPbs in Fig. 9c) is mapped significantly deeper.
<figure><img src="image_19.png" /><figcaption>**Figure 11.** Processing waveforms of teleseismic earthquakes recorded at station ME25 from the FAME Mendocino array. (a) Original vertical components after signal-to-noise ratio selection, low-pass filtering at 1 Hz, alignment according to the theoretical arrival time of the *P* wave in the reference *iasp91* velocity model (Kennett & Engdahl 1991), and sorting by increasing slowness of the *P* wave. (b) Same vertical components after spectral whitening. (c) Auto-correlation of the spectrally whitened vertical components sorted as a function of the *P*-wave slowness. Red curves are the expected delays for P-reflections atop 20- and 30-km deep discontinuities using the *iasp91* model. (d) Same as (c) except that the autocorrelations are sorted by backazimuth. The two red arrows point anomalous sloping arrivals possibly related to 3-D propagation effects. The inset at the top right hand side locates the FAME Mendocino broad-band experiment (black triangles) in northern California, and the ME25 station (red triangle). This station is located at 40.3°N latitude, 121.8°W longitude. The contours in black locate at depth the subducted Gorda Plate (McCrory *et al.* 2014). The blue line locates the profile associated with seismic sections shown in Figs 12(c) and (d). MTJ: Mendocino Triple Junction.</figcaption></figure>

Another possibility is that the dual interface results from the difference in sensitivity to velocity gradients of the reflection coefficient, $R_{pp}^{1→2}$, and the coefficients associated with modal changes, $T_{ps}^{2→1}$, $R_{ps}^{1→2}$ and $R_{ss}^{1→2}$, such as described in Section 2.6. The dipping

portion of the plate is also a likely problem for the imaging. The imaging principle is indeed based on ray-theory, uses a 1-D model, and does not account for the elastic scattering patterns. In addition, if the structure is not perfectly 2-D, such as
<figure><img src="image_20.png" /><figcaption>**Figure 12.** Seismic imaging from a dense broadband seismic network. (a) Time-to-depth converted broad-band vertical RFs for station ME25 of the FAME Mendocino experiment. (b) Depth converted radial RFs. The blue and red lines mark discontinuities at 20 and 30 km depth, respectively. (c) Section obtained from multimode CCP stacking of radial RFs for the whole FAME Mendocino broad-band seismic array. The interpretation of first-order structural discontinuities (i.e. Tauzin et al. 2017) is superimposed. (d) Result of common reflection point stacking the vertical RFs for this same network. We have flipped the signal in polarity. The first order structures are recovered, that is the dipping red interface at west and the subhorizontal continental Moho at east.</figcaption></figure>

suggested by the plate contours in the inset of Fig. 11, it is likely
that the image is contaminated by effects of projection of the struc-
ture in the direction perpendicular to the profile. A method rely-
ing on 3-D Kirchhoff pre-stack depth migration may help to im-
prove the images (Cheng et al. 2016, 2017; Hansen & Schmandt
2017).

6 DISCUSSION

6.1 Short-period single-component sensors

An advantage of applying spectral whitening and auto-correlation
to single component records is that we are not anymore
dependent on the limited geographical distribution of three-component stations to make detailed regional studies of the seismic structure. We can now make use of the records from single-component sensors. One-component instruments have been densely deployed in some regions for a very long duration. For instance, the Northern California Seismic Network (NCSN) shown in the map in Fig. 13(a) is a permanent network that has operated since 1967 (USGS, 1967). Because it mostly comprises short-period (vertical) single-component sensors, it has never been exploited to extract the structural response associated with near-vertically propagating seismic waves. Within the 30–90° range of epicentral distances, we recovered 435 735 raw waveforms from 487 stations over the period 1985–2016. Several hundred of broad-band stations are available in the same region, but with a very different spatial distribution, achieving for instance a less dense coverage of the San Andreas fault zone.

We can wonder if the narrow bandpass of short-period instruments at frequencies above ~1 Hz represents a limitation for RF estimates. To investigate this question, we show in Fig. 13 the result of processing the vertical-component records of a large (magnitude 7.6) earthquake. Fig. 13(b) shows a subsample of the high-quality earthquake waveforms, with snr of raw seismograms >2.0, analysed within the bandpass of the instruments and within the 3–10 Hz band. In Fig. 13(c), we represent in grey the smoothed spectra (eq. 14) of all the available vertical components, and an example of applying spectral whitening to one of the traces, in blue before, and black after processing.

Within the bandpass of the instruments (Fig. 13b, left), this event gives high-quality waveforms at 300 stations of the NCSN network, with a low noise level before the teleseismic P-wave arrival. However above 3 Hz (Fig. 13b, right), the records are largely dominated by the seismic ambient noise. This is confirmed in Fig. 13(c) where most (grey) spectra are narrowly focused within 1 decade in the 0.1–4 Hz bandpass, whereas above 4 Hz the spectra spread over 2–4 decades revealing variable ambient noise conditions or site effects. Due to noise characteristics, most teleseismic (i.e. almost vertical incidence) energy is observed below 4 Hz. This particular event shows that for short-period sensors, the usable teleseismic band is 0.2–4 Hz, the upper frequency limit being significantly higher than used for typical RF studies. We require here the Gaussian filter to cut frequencies beyond ~3 Hz and apply spectral whitening and autocorrelation (Fig. 13 c; blue and black spectra). The result is shown in the time-domain in Fig. 13(d) for a subsample of the data sorted by epicentral distances. The general characteristics of these vertical RFs are close to the ones obtained from the broadband station ME25 in Fig. 6, although higher frequency. Due to the analysis of data associated with a common source, it is possible that part of the reflectivity is related to the near-source structure. However, the autocorrelograms differ from a station to another, suggesting that the change in nature of reflectivity is also related to spatial variations of the structure below the receivers. In conclusion, our approach can be applied not only to process broad-band data, but it also applies to records from short-period single-component sensors up to the frequency of ~4 Hz.

A key assumption in the seismic interferometry principle applied here is that seismic energy comes with near-vertical incidence (Claerbout 1968). That is why we use the coda of teleseismic earthquakes to extract the receiver side response. We have shown that for a magnitude 7.6 earthquake, the 0.1–4 Hz frequency band corresponds to the teleseismic energy, whereas frequencies above 4 Hz are associated with the seismic ambient noise. This noise is often attributed to shallow anthropogenic or environmental sources

that involve propagation of surface waves (e.g. Nakamura 1989; Bonnefoy-Claudet et al. 2006).

Some authors make use of correlations of continuous records of seismic ambient noise in the 0.5–4 Hz frequency band to recover the body-wave reflectivity (e.g. Poli et al. 2012; Gorbatov et al. 2013; Kennett 2015; Becker & Knapmeyer-Endrun 2018). This naturally leads to the question of the respective contributions of teleseismic earthquakes and ambient noise in the apparently successful extraction of vertically propagating body-wave reflectivity.

## 6.2 Noisy earthquake codas

The seismic ambient noise has long been considered as a nuisance for studying the Earth through earthquake responses. It is only recently that seismic interferometry methods have been developed to turn the noise from distant sources into a transient signal between seismometers for structural Green's functions retrieval (e.g. Roux et al. 2005; Shapiro et al. 2005). An advantage of using noise records compared to earthquake records is that the deployment time of seismic arrays to recover the Earth's seismic response can significantly be reduced. However, the process responsible for the emergence of body waves in noise correlations is debated (e.g. Ruigrok et al. 2011; Poli et al. 2012; Lin & Tsai 2013; Lin et al. 2013). Theoretical and empirical studies show that to obtain body waves from noise sources on the free surface, these must be recorded at teleseismic distances and include a significant component of vertically propagating body waves in the source mechanism. Array analysis suggests that a component of the seismic noise field indeed consists in compressional body waves originating from large storms at sea and coupling of the ocean with the solid Earth (e.g. Lacoss et al. 1969; Vinnik 1973; Ruigrok et al. 2011). However, Lin et al. (2013) also found that the amplitudes of body waves retrieved from noise correlations are highly correlated with the amount of earthquake-excited energy present in the continuous noise records. Earthquake body-wave amplitudes seem to be an integral part of the ambient seismic noise field. This component seems prevalent at frequencies higher than 0.3 Hz (Lacoss et al. 1969; Ruigrok et al. 2011). This likely explains why studies using ambient noise often make use of the 0.5–4 Hz frequency band to recover crustal body-wave reflectivity (e.g. Poli et al. 2012; Gorbatov et al. 2013; Kennett 2015; Becker & Knapmeyer-Endrun 2018).

Here, we are interested in examining the effect of the noise on the retrieval of the response from noisy earthquake codas. We therefore investigate the special case where the energy in the teleseismic coda is present but cannot be distinguished from the noise. This can occur either by a decrease of the signal in the teleseismic band in Fig. 13(c), an increase of the noise level over the whole spectrum, or both. The first case occurs when examining teleseismic records of small magnitude earthquakes: with decreasing magnitude, the phase of the earthquake waveform remains the same but the amplitude decreases to the point that it does not emerge from the noise. The second case occurs when considering large earthquakes—in our case beyond magnitude 5.5—but with a high level of noise at the recording site.

We designed simple experiments to investigate the effect of the noise. Our objectives are multifold: (i) to extract a single-station response from vertical-component records of teleseismic earthquakes, (ii) to check if it is possible to recover this response from stacking waveforms with poor signal-to-noise ratios (snr) and (iii) to estimate the minimal number of data to stack to obtain stable responses.
<figure><img src="image_21.png" /><figcaption>**Figure 13.** Processing short-period vertical-component records from the NCSN network in California. (a) Map of the 300 NCSN short-period seismometers that have recorded high-quality waveforms for the magnitude 7.6, 2014 April 12, teleseismic earthquake. Red triangles locate Quaternary active volcanic fields. The inset locates the map on the globe. (b) A subsample of the raw waveforms (left), and filtered waveforms with a 4th-order butterworth filter with a 3–10 Hz bandpass (right). The vertical components are normalized to the maximum amplitude within the 300 s time window. The instrumental response is not removed here. (c) A spectral analysis of the records and an example of application of spectral whitening on the first seismic trace. The traces in grey correspond to the smoothed spectra of the 300 vertical components (eq. 14). Thick and thin red lines show the average and 1 standard deviation on the smoothed spectra distribution. The spectrum in blue corresponds to the raw waveform of the first trace in (b). We divide this spectrum by its regularized version, and multiply it by the gaussian filter in green, to obtain the spectrally whitened spectrum in black (eq. 13). Three bands of frequencies can be distinguished, delimited by vertical dashed lines: the seismometers have little sensitivity to the lowest frequencies ($10^{-2}$–$10^{-1}$ Hz); the smoothed spectra (grey) are the most similar to each other and show a peak in amplitude between 0.1 Hz and 3 Hz, corresponding to the band of teleseismic energy; the spectra more significantly differ in the high-frequency band, associated with variable ambient noise conditions and site effects. (d) High-frequency (up to ~3 Hz) vertical autocorrelations in the time domain sorted by increasing epicentral distance.</figcaption></figure>

We identified within the NCSN network the station having recorded the largest number of teleseismic earthquakes with magnitude larger than 5.5 (GCMT catalogue): the station NSM, whose location is shown in magenta in the map of Fig. 13(a), has recorded 1166 events over the period 1997–2016, with variable snr (Fig. 14a). We estimated the snr similarly as in Section 5, from the ratio of rms of the signal in the coda over the noise before the P-arrival. The distribution of snr has a maximum between 1 and 1.1, and is skewed toward higher values (Fig. 14a). We show in Fig. 14(b) the

corresponding waveforms. Visually, good quality waveforms are obtained for snr larger than 2.0. However, 94 per cent of the data have a snr between 0.8 and 1.4 (Fig. 14a), for which it is difficult to observe body waves (Fig. 14b). We find a positive correlation (+0.29) between the magnitude of the earthquake and the snr, but the small value suggests that, in California, the quality of the waveforms mainly depend on the ambient noise conditions rather than the earthquake magnitude.
<figure><img src="image_22.png" /><figcaption>**Figure 14.** Effect of the noise on the retrieval of the receiver-side responses from earthquake codas. (a) A histogram of signal-to-noise ratios (snr) for the 1166 short-period vertical-component records at the station NSM from the NSCN network in California. 94 per cent of the data are within 0.8–1.4 snr. (b) A subsample of the raw waveforms sorted by increasing order of snr. We emphasize the range of snr for 94 per cent of the data. (c) Result of applying spectral whitening and autocorrelation on the waveforms shown in (b). The autocorrelograms are stacked by bins and ordered as a function of the snr in the initial waveforms. (d) Stack for the ensemble of data. The amplitude is normalized relative to the peak at the zero-lag. Red plusses and blue minuses emphasize the potential arrivals associated with reflections at sharp interfaces below the NSM station. The early part of the autocorrelograms, up to ~3 s, is dominated by the negative sidelobe from the peak at the zero lag.</figcaption></figure>

Building a diagram of the spectra similarly as in Fig. 13(c) (Fig. S2), we found as expected a cross-over in the spectra of high- and low-quality waveforms. This cross-over is near 2–3 Hz, with high quality waveforms (snr above 1.4) having high-energy in the teleseismic band and low-noise level above 2–4 Hz. Low-quality waveforms (snr below 0.8) have a small surviving peak in the tele-seismic band, but the energy of the noise above 3 Hz is as strong as the energy in the teleseismic band. In general, even for a low snr, the earthquake waveforms remain significantly enriched in high fre-quencies in the teleseismic band (0.5–1.0 Hz) compared to typical spectra for the microseismic noise (e.g. Ruigrok et al. 2011).

We then applied spectral whitening and autocorrelation on the 1166 vertical component waveforms recorded at the station NSM. We show in Fig. 14(c) the result of stacking the corresponding vertical RFs in non-overlapping bins with fixed width (0.05) of snr. Given the distribution in Fig. 14(a), the number of stacked data is highly variable, and is the largest in bins of snr between 0.8 and 1.4 (up to 287 waveforms). The behavior of the responses stabilizes with the number of stacked data, independently of the initial quality of the

data (Fig. 14c). In the 0.8-1.4 snr range, clear coherent arrivals are observed on the autocorrelograms. These arrivals are still visible in the phase of higher quality waveforms (snr > 2.0), but blurred in the middle of larger amplitude signals with no apparent coherency from an autocorrelogram to another. This behavior is in agreement with the one obtained from our synthetic experiment in Fig. 7, suggesting that the stacking operation is critical to improve the ‘signal-to-noise ratio’ of auto-correlograms. We show in Fig. 14(d) the average vertical RF for station NSM, resulting from the stack of all auto-correlograms. We also emphasize high-frequency arrivals that may correspond to reflections. By minimizing the least square difference between this averaged autocorrelogram and the stacks of a variable number of autocorrelograms in the distribution, we found that the response remains stable after 10 to 20 stacked autocorrelograms.

This result has significant implications. The typical implementa-
tion of RF studies is to select the best quality waveforms prior to
deconvolution. In this context, it is clear that 94 per cent of the data
(within the 0.8 and 1.4 snr range) would have been considered as bad
quality data, and not used for deconvolution and stacking. We show
<figure><img src="image_23.png" /><figcaption>**Figure 15.** Vertical receiver functions for low-magnitude earthquakes recorded at the NSM short-period sensor. (a) Number of vertical component records (blue histogram) and logarithm of the cumulative number of data (black dots) per bins of magnitude in the USGS catalogue. The fit to the cumulative number of data (red straight line) suggests that the catalogue is complete for such magnitudes. (b) Stacked high-frequency autocorrelograms from low-magnitude event-based data. We emphasize with blue and red lines arrivals that can possibly be associated with reflections atop negative and positive discontinuities.</figcaption></figure>

here that stacking the numerous data in this snr range gives however coherent signals in the autocorrelograms. This raises a question of quantity versus quality, and a possibility that using lower magnitude earthquakes would improve our imaging capabilities.

Indeed the gain in data would be considerable. Globally, for a de-
crease of 1 unit in magnitude, the number of available earthquakes
multiplies roughly by 10. A limitation in using low-magnitude earth-
quakes is in the completeness of global catalogues. An examination
of the USGS catalogue indicates that magnitudes above 4.5 are usu-
ally well documented (Fig. 15a). Therefore, in addition to our pre-
liminary 1166 waveforms (magnitude larger than 5.5; GCMT cata-
logue), we could obtain 15 620 supplementary vertical-component
records of earthquakes with magnitude between 4.5 and 5.4 (USGS
catalogue) at the station NSM. 3759 records are available for earth-
quakes with a magnitude of 4.5 (Fig. 15a). In comparison, only
142 earthquakes would be used in a typical RF study because they
present a sufficient magnitude (>5.5) and a good snr (>2.0).

We show in Fig. 15(b) the result of stacking autocorrelograms
by non-overlapping bins of magnitude between 4.5 and 5.4. The
stacked autocorrelograms appear stable regardless of the magni-
tude. Likewise, the phase of these traces is comparable with the
one obtained in the stack for earthquakes with magnitude larger
than 5.5. Based on this observation, we conclude that the use of
low-magnitude earthquakes for reconstructing structural responses
from seismic interferometry seems feasible and promising. How-
ever, when processing such data, care must be taken with regards
to the effect of instrumental noise. Indeed, we realized that over its

period of activity (1997–2016), a significant monochromatic instru-
mental noise at specific frequencies of 1, 2 and 3 Hz was present
at the station NSM. We show in Fig. S3 the result of stacking these
data without discarding spectra with amplitude larger than 1 decade
from the background level at these frequencies. Large amplitude
signals are observed at regular time intervals in the late portion
of the autocorrelations, mainly associated with the 1 Hz compo-
nent. Our empirical data selection succeeds to partly remove this
component (Fig. 15b).

We have shown that stacking autocorrelograms from noisy earth-
quake codas can be used to recover an approximate vertical re-
flectivity within the framework of the seismic interferometry prin-
ciple from Claerbout (1968). The autocorrelograms seem stable
although slightly enriched in low-frequencies for magnitudes larger
than 5.5 (Fig. 15b). The recovered amplitudes however largely dif-
fer. The origin of this difference might have a simple explana-
tion: the records of low-magnitude teleseismic earthquakes con-
sists in a mixture of a vertically propagating (teleseismic) wave-
field and waves propagating with a large horizontal component,
such as surface waves (e.g. Bonnefoy-Claudet et al. 2006). This
deviation from near-vertical incidence departs from the original as-
sumption in the seismic interferometry principle (Claerbout 1968;
Wapenaar 2003; Galetti & Curtis 2012). This likely introduces
a bias in the recovery of the body-wave reflectivity in a simi-
lar manner as the inhomogeneous distributions of surface noise
sources affects the retrieval of surface waves (e.g. Paul et al.
2005).
## 6.3 Seismic ambient noise

We conclude this section by examining at the same station, NSM, the nature of the responses recovered through seismic ambient noise autocorrelations. We requested continuous waveforms for the first 6 months of the year 2004. This year was chosen because it is free from the problematic 1, 2 and 3 Hz monochromatic instrumental noise. Following Kennett (2015), 6 months of record is usually adequate to obtain a stable autocorrelogram. The 6-month record was subsequently cut into 4-hr windows. This is shorter than the 6-hr windows used by Gorbatov *et al.* (2013), Kennett (2015) and Kennett *et al.* (2015), but longer than the 1–3 hr used by Becker & Knapmeyer-Endrun (2018). For computational purposes, we down-sampled the waveforms from 100 to 10 points per second (Kennett 2015). We then applied the processing from Pham & Tkalčić (2017)—that is eqs (13) and (14)—to each 4-hr time segments, making sure that the window width for spectral whitening (2*N* + 1, with *N* the number of frequency samples) is the same as in processing the event-based data. The Gaussian filter is also the same, cutting frequencies above ~3 Hz. We found a significant artificial long-period component in the recovered responses that we cut out with a Butterworth high-pass filter with 0.33 Hz corner frequency. Consequently, there is a difference in the frequency content and time sampling in the results from event- and noise-based data. We do not thoroughly discuss the alternatives to this processing, but rather focus on the comparison with earthquake data. The studies referred above give examples of different processing.

We show in Fig. 16(a) the daily autocorrelograms, obtained by stacking the 4-hr time window autocorrelograms in 2004. The result is stable across the whole year. In Fig. 16(b), we make a comparison of the stacked noise-based autocorrelations (black trace) with the autocorrelations obtained from large magnitude earthquake data with high signal-to-noise ratio (red trace), and the autocorrelations obtained from smaller magnitude earthquakes (blue trace). The noise-based trace has large side-lobs near the peak at zero-lag because of the narrow filtering and shows a more aliased appearance due to a coarser time sampling. The waveform of the noise-based autocorrelogram share some similarities with the one from low-magnitude earthquakes, but largely differ from the one for large magnitude earthquakes. Possibly, the noise component and its directivity affect the retrieval of the body-wave reflectivity.

## 7 CONCLUSIONS

This paper aimed to review the concepts underlying the well-known RFs method, and points to an analogy with seismic interferometry in a pedagogic way. By using a data processing method that involves spectral whitening and auto/cross-correlation, we have shown that the recovered RFs (or receiver-side responses) are similar to the ones obtained from more classical deconvolution methods. Apart from providing the radial RFs, this approach also provides the vertical RFs including the multiple reflected *P*-waves in the stratification beneath the receivers. The key assumption in seismic interferometry is that the seismic energy is recorded with near-vertical incidence (Claerbout 1968). A remaining limitation is the contamination by the source side structure due to imperfect destructive interference when stacking the empirical responses. We have shown an example of 'interferometric seismic imaging' to illuminate the crustal and mantle structure from the reconstructed vertical and horizontal seismic responses. We also have shown that the processing can be applied to records from short-period single-component sensors, up to the frequency of ~4 Hz, with a benefit in terms of resolution and

<figure><img src="image_24.png" /><figcaption>**Figure 16.** Vertical receiver functions from the seismic ambient noise and comparison with earthquake data. (a) Daily vertical autocorrelograms at the station NSM over a 6-month period. Blue colors represent positive amplitudes and red colors represent negative amplitudes. We could not compute autocorrelograms during 4 d (vertical red traces) due to an absence of contiguous 4-hr long time windows on these days. (b) A comparison between the autocorrelograms obtained from seismic ambient noise (black), and earthquake data with magnitudes larger than 5.5 (red) or below 5.5 (blue). The time step and bandpass for the noise study differ slightly from the ones for the event-based studies.</figcaption></figure>

coverage of the crustal structure. Finally, utilizing the same interferometry principle, we have established a link between event- and noise-based seismic interferometry. We make a case that approximate receiver-side reflectivity can be retrieved from low-magnitude earthquakes. We also take into consideration autocorrelations from the continuous seismic ambient noise wavefield. We speculate that differences that arise among different methods are related to the directivity of the seismic wavefield.

## ACKNOWLEDGEMENTS

Waveform data of the National California Seismic Network were accessed through the Northern California Earthquake Data Center (NCEDC), doi:10.7932/NCEDC. The IRIS Data Management Center was used for access to waveforms of the FAME Mendocino Experiment. IRIS Data Services are funded through the Seismological Facilities for the Advancement of Geoscience and EarthScope (SAGE) Proposal of the National Science Foundation under Cooperative Agreement EAR-1261681. We thank the universities of Oregon, California at Berkeley, Rice and the Earthscope program of the National Science Foundation for acquiring and providing the FAME waveform data. B.T. was funded by a Délégation CNRS and
Congé pour Recherches et Conversion Thématique from the Université de Lyon to visit the Research School of Earth Sciences (RSES), The Australian National University (ANU). We thank the CNRS, the Université de Lyon, and ANU, for providing the funding and research environment. We thank Brian L.N. Kennett and Malcolm Cambridge for constructive discussions about the concepts underlying RF analysis and seismic interferometry. We thank Brian L.N. Kennett for his insightful advices on plane-wave matrix algorithms for the computation of the effective response of stratified media.

## REFERENCES

Abers, G.A., MacKenzie, L.S., Rondenay, S., Zhang, Z., Wech, A.G. & Creager, K.C., 2009. Imaging the source region of cascadia tremor and intermediate-depth earthquakes, *Geology*, **37**(12), 1119–1122.

Abe, S., Kurashimo, E., Sato, H., Hirata, N., Iwasaki, T. & Kawanaka, T., 2007. Interferometric seismic imaging of crustal structure using scattered teleseismic waves, *Geophys. Res. Lett.*, **34**(19).

Ainiwaer, A. & Gurrola, H., 2017. Wavefield iterative deconvolution to remove multiples and produce phase specific ps receiver functions, *Geophys. J. Int.*, **212**(3), 2016–2030.

Aki, K. & Richards, P., 1980. *Quantitative Seismology*, W. H. Freeman.

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

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

Audet, P. & Schaeffer, A.J., 2018. Fluid pressure and shear zone development over the locked to slow slip region in cascadia, *Sci. Adv.*, **4**(3), eaar2982.

Becker, G. & Knapmeyer-Endrun, B., 2018. Crustal thickness across the trans-european suture zone from ambient noise autocorrelations, *Geophys. J. Int.*, **212**(2), 1237–1254.

Bensen, G., Ritzwoller, M., Barmin, M., Levshin, A., Lin, F., Moschetti, M., Shapiro, N. & Yang, Y., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, *Geophys. J. Int.*, **169**(3), 1239–1260.

Bonnefoy-Claudet, S., Cotton, F. & Bard, P.-Y., 2006. The nature of noise waveform and its applications for site effects studies: A literature review, *Earth-Sci. Rev.*, **79**(3-4), 205–227.

Bostock, M., 1999. Seismic waves converted from velocity gradient anomalies in the Earth's upper mantle, *Geophys. J. Int.*, **138**, 747–756.

Bostock, M., Hyndman, R., Rondenay, S. & Peacock, S., 2002. An inverted continental moho and serpentinization of the forearc mantle, *Nature*, **417**(6888), 536–538.

Brune, J.N., 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes, *J. geophys. Res.*, **75**(26), 4997–5009.

Cassidy, J., 1992. Numerical experiments in broadband receiver function analysis, *Bull. seism. Soc. Am.*, **82**(3), 1453–1474.

Chapman, C., 2003. Yet another elastic plane-wave, layer-matrix algorithm, *Geophys. J. Int.*, **154**(1), 212–223.

Chapman, C., 2004. in *Fundamentals of Seismic Wave Propagation*, Cambridge University Press.

Cheng, C., Bodin, T. & Allen, R.M., 2016. Three-dimensional pre-stack depth migration of receiver functions with the fast marching method: a kirchhoff approach, *Geophys. J. Int.*, **205**(2), 819–829.

Cheng, C., Bodin, T., Tauzin, B. & Allen, R.M., 2017. Cascadia subduction slab heterogeneity revealed by three-dimensional receiver function kirchhoff migration, *Geophys. Res. Lett.*, **44**(2), 694–701.

Claerbout, J.F., 1968. Synthesis of a layered medium from its acoustic transmission response, *Geophysics*, **33**(2), 264–269.

Clayton, R.W. & Wiggins, R.A., 1976. Source shape estimation and deconvolution of teleseismic bodywaves, *Geophys. J. Int.*, **47**(1), 151–177.

Crotwell, H.P., Owens, T.J. & Ritsema, J., 1999. The taup toolkit: Flexible seismic travel-time and ray-path utilities, *Seism. Res. Lett.*, **70**(2), 154–160.

Galetti, E. & Curtis, A., 2012. Generalised receiver functions and seismic interferometry, *Tectonophysics*, **532**, 1–26.

Gorbatov, A., Saygin, E. & Kennett, B., 2013. Crustal properties from seismic station autocorrelograms, *Geophys. J. Int.*, **192**(2), 861–870.

Guan, Z. & Niu, F., 2017. An investigation on slowness-weighted ccp stacking and its application to receiver function imaging, *Geophys. Res. Lett.*, **44**, 6030–6038.

Hansen, S. & Schmandt, B., 2017. P and s wave receiver function imaging of subduction with scattering kernels, *Geochem., Geophys., Geosyst.*, **18**(12), 4487–4502.

Hansen, S.M., Dueker, K.G., Stachnik, J.C., Aster, R.C. & Karlstrom, K.E., 2013. A rootless rockies? support and lithospheric structure of the col- orado rocky mountains inferred from crest and ta seismic data, *Geochem., Geophys., Geosyst.* **14**(8), 2670–2695.

Helffrich, G. & Wood, J., 1996. 410 km discontinuity sharpness and the form of the olivine α-β phase diagram: Resolution of apparent seismic contradictions, *Geophys. J. Int.*, **126**, F7–F12.

Hetényi, G., 2007. Evolution of deformation of the Himalayan prism: from imaging to modelling, *PhD thesis*, Université Paris Sud-Paris XI.

Houston, H. & Kanamori, H., 1986. Source spectra of great earthquakes: teleseismic constraints on rupture process and strong motion, *Bull. seism. Soc. Am.*, **76**(1), 19–42.

Kennett, B., 1983. in *Seismic Wave Propagation in Stratified Media*, Cambridge University Press.

Kennett, B., 2015. Lithosphere–asthenosphere p-wave reflectivity across australia, *Earth planet. Sci. Lett.*, **431**, 225–235.

Kennett, B. & Engdahl, E., 1991. Travel times for global earthquake location and phase identification, *Geophys. J. Int.*, **105**, 429–465.

Kennett, B. & Pham, T.-S., 2018. Evolution of the correlation wavefield extracted from seismic event coda, *Phys. Earth planet. Inter.*, **282**, 100–109.

Kennett, B., Saygin, E. & Salmon, M., 2015. Stacking autocorrelograms to map moho depth with high spatial resolution in southeastern Australia, *Geophys. Res. Lett.*, **42**(18), 7490–7497.

Kikuchi, M. & Kanamori, H., 1982. Inversion of complex body waves, *Bull. seism. Soc. Am.*, **72**(2), 491–506.

Lacoss, R.T., Kelly, E.J. & Toksöz, M.N., 1969. Estimation of seismic noise structure using arrays, *Geophysics*, **34**(1), 21–38.

Langston, C., 1979. Structure under mont Rainier, Washington, inferred from teleseismic body waves, *J. geophys. Res.*, **84**.

Lay, T. & Wallace, T.C., 1995. *Modern Global Seismology*, Vol. 58, Elsevier.

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

Lin, F.-C. & Tsai, V.C., 2013. Seismic interferometry with antipodal station pairs, *Geophys. Res. Lett.*, **40**(17), 4609–4613.

Lin, F.-C., Tsai, V.C., Schmandt, B., Duputel, Z. & Zhan, Z., 2013. Extracting seismic core phases with array interferometry, *Geophys. Res. Lett.*, **40**(6), 1049–1053.

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

Li, X.-Q. & Nábělek, J.L., 1999. Deconvolution of teleseismic body waves for enhancing structure beneath a seismometer array, *Bull. seism. Soc. Am.*, **89**(1), 190–201.

Ma, Y., Wang, R. & Zhou, H., 2012. A note on the equivalence of three major propagator algorithms for computational stability and efficiency, *Earthq. Sci.*, **25**(1), 55–64.

McCrory, P.A., Hyndman, R.D. & Blair, J.L., 2014. Relationship between the cascadia fore-arc mantle wedge, nonvolcanic tremor, and the downdip limit of seismogenic rupture, *Geochem., Geophys., Geosyst.*, **15**(4), 1071–1095.

Nakamura, Y., 1989. A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface, *QR Railway Tech. Res. Inst.*, **30**(1), 25–33.

Oren, C. & Nowack, R., 2017. Seismic body-wave interferometry using noise auto-correlations for crustal structure, *Geophys. J. Int.*, **208**(1), 321–332.

Paul, A., Campillo, M., Margerin, L., Larose, E. & Derode, A., 2005. Empirical synthesis of time-asymmetrical Green's functions from the correlation of coda waves, *J. geophys. Res.: Solid Earth*, **110**(B8).
Pham, T.-S. & Tkalčić, H., 2017. On the feasibility and use of teleseismic p-wave coda autocorrelation for mapping shallow seismic discontinuities, *J. geophys. Res.: Solid Earth*, **122**, 3776–3791.

Poli, P., Pedersen, H. & Campillo, M., 2012. Emergence of body waves from cross-correlation of short period seismic noise, *Geophys. J. Int.*, **188**(2), 549–558.

Randall, G., 1989. Efficient calculation of differential seismograms for lithospheric receiver functions, *Geophys. J. Int.*, **99**(3), 469–481.

Reading, A., Kennett, B. & Sambridge, M., 2003. Improved inversion for seismic structure using transformed, s-wavevector receiver functions: removing the effect of the free surface, *Geophys. Res. Lett.*, **30**(19).

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

Rondenay, S., Bostock, M. & Shragge, J., 2001. Multiparameter two-dimensional inversion of scattered teleseismic body waves 3. Application to the cascadia 1993 data set, *J. geophys. Res.: Solid Earth*, **106**(B12), 30795–30807.

Roux, P., Sabra, K.G., Gerstoft, P., Kuperman, W. & Fehler, M.C., 2005. P-waves from cross-correlation of seismic noise, *Geophys. Res. Lett.*, **32**(19).

Ruigrok, E., Campman, X. & Wapenaar, K., 2011. Extraction of p-wave reflections from microseisms, *Compt. Rendus Geosci.*, **343**(8-9), 512–525.

Ruigrok, E. & Wapenaar, K., 2012. Global-phase seismic interferometry unveils p-wave reflectivity below the himalayas and tibet, *Geophys. Res. Lett.*, **39**(11).

Schimmel, M. & Paulssen, H., 1997. Noise reduction and detection of weak, coherent signals through phase-weighted stacks, *Geophys. J. Int.*, **130**(2), 497–505.

Schuster, G., Yu, J., Sheng, J. & Rickett, J., 2004. Interferometric/daylight seismic imaging, *Geophys. J. Int.*, **157**(2), 838–852.

Shapiro, N.M., Campillo, M., Stehly, L. & Ritzwoller, M.H., 2005. High-resolution surface-wave tomography from ambient seismic noise, *Science*, **307**(5715), 1615–1618.

Sun, W. & Kennett, B., 2016. Receiver structure from teleseisms: autocorrelation and cross correlation, *Geophys. Res. Lett.*, **43**(12), 6234–6242.

Sun, W. & Kennett, B., 2017. Mid-lithosphere discontinuities beneath the western and central north china craton, *Geophys. Res. Lett.*, **44**(3), 1302–1310.

Tauzin, B., Bodin, T., Debayle, E., Perrillat, J.-P. & Reynard, B., 2016. Multi-mode conversion imaging of the subducted gordà and juan de fuca plates below the north american continent, *Earth planet. Sci. Lett.*, **440**, 135–146.

Tauzin, B., Reynard, B., Perrillat, J.-P., Debayle, E. & Bodin, T., 2017. Deep crustal fracture zones control fluid escape and the seismic cycle in the cascadia subduction zone, *Earth planet. Sci. Lett.*, **460**, 1–11.

Tibuleac, I.M. & von Seggern, D., 2012. Crust-mantle boundary reflectors in nevada from ambient seismic noise autocorrelations, *Geophys. J. Int.*, **189**(1), 493–500.

Tseng, T.-L., Chen, W.-P. & Nowack, R.L., 2009. Northward thinning of tibetan crust revealed by virtual seismic profiles, *Geophys. Res. Lett.*, **36**(24).

USGS Menlo Park, 1967. International Federation of Digital Seismograph Networks. Other/Seismic Network, 10.7914/SN/NC.

Vinnik, L., 1973. Sources of microseismic P waves, *Pure appl. Geophys.*, **103**(1), 282–289.

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

Wapenaar, K., 2003. Synthesis of an inhomogeneous medium from its acoustic transmission response, *Geophysics*, **68**(5), 1756–1759.

Wittlinger, G. et al., 2004. Teleseismic imaging of subducting lithosphere and moho offsets beneath western Tibet, *Earth planet. Sci. Lett.*, **221**, 117–130.

Yu, C.-Q., Chen, W.-P., Ning, J.-Y., Tao, K., Tseng, T.-L., Fang, X.-D., Chen, Y.J. & van der Hilst, R.D., 2012. Thick crust beneath the ordos plateau: Implications for instability of the North China craton, *Earth planet. Sci. Lett.*, **357**, 366–375.

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**(3), 1932–1941.

Zhu, L., 2000. Crustal structure across the San Andreas fault, southern California, from teleseismic converted waves, *Earth planet. Sci. Lett.*, **179**, 183–190.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

**Figure S1.** (a) Vertical and (b) radial RFs calculated with spectral whitening and correlation of the records of good signal-to-noise ratio teleseismic P-wave codas at station ME25. Red curves in (a) are the theoretical delays for P-reflections atop 20 km and 30 km deep discontinuities in iasp91. Blue, magenta and green curves are traveltime predictions for the direct P-S conversions, and P-P-S and P-S-S reverberations. The arrivals near 3 s after the P-wave, associated with the direct conversions at the two discontinuities are remarkable.

**Figure S2.** Spectral analysis of the records of magnitude >5.5 teleseismic earthquakes at station NSM from the NCSN network in California. The smoothed spectra are color coded as a function of the signal-to-noise ratio (snr), with high snr shown in red and low snr shown in blue. The two thick lines in black and red correspond to spectra with snr of 0.8 and 1.4, respectively. 94 per cent of the data are within this range of snr. The teleseismic band of frequency corresponds to data with the highest quality (red), with spectra showing a peak up to 6 Hz before flattening due to the seismic ambient noise level. There is an intersection of the spectra with low (black) and high (red) snr at 2–3 Hz (dashed black line). Even for low snr, the earthquake waveforms remain significantly enriched in high frequencies in the teleseismic band (0.5–1.0 Hz) compared to typical spectra for the microseismic noise.

**Figure S3.** Vertical receiver functions for low-magnitude earthquakes recorded at the NSM short-period sensor, without removing the records contaminated by the 1, 2 and 3 Hz monochromatic components. (a) The number of vertical component records (blue histogram) and logarithm of the cumulative number of data (black dots) per bins of magnitude in the USGS catalogue. The fit to the cumulative number of data (red straight line) suggests that the catalogue is complete for such magnitudes. (b) Stacked high-frequency autocorrelograms from low-magnitude event-based data. Blue and red lines emphasize arrivals that can possibly be associated with reflections atop negative and positive discontinuities.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.