https://orcid.org/0000-0002-4358-1189
SUMMARY
Teleseismic receiver functions are widely used to map the depth and topography of various major discontinuities in the Earth’s mantle. To determine what precisely contributes to the receiver functions, we applied the adjoint method of full waveform inversion to calculate their sensitivity kernels. These kernels illustrate the extend to which model parameters may influence the waveforms. We calculated synthetic data for a realistic event measured at a realistic receiver array, whereby we focused on the waveforms of the P410s and P660s phases, that convert a P to an S wave at the 410- and 660-discontinuity, respectively. We calculated both the volumetric sensitivity kernels for density, P- and S-wave speeds, as well as boundary kernels that illustrate receiver functions’ sensitivity to topography on the discontinuity. In the boundary kernels, we observe that receiver functions are highly sensitive to a discontinuity’s topography, in particular to an area surrounding the conversion point with a radius comparable to the Fresnel zone. However, the volumetric kernels illustrate a sensitivity to model parameters in large areas of the mantle. This includes sensitivity to the Fresnel zone of the converted wave far before the conversion, as well as sensitivity to scatterers of other phases. We therefore conclude that receiver functions are sensitive to the topography of discontinuities. However, effects of an incorrect velocity model, even far from the conversion point, may erroneously be projected onto the topography of the discontinuity. Therefore, a simultaneous inversion of topography and velocity parameters is required to image topography with high accuracy.
1 INTRODUCTION
It is well known that the Earth’s interior consists of various layered regions separated by large jumps of seismic wave speeds, indicative of a phase transition or a major change in composition. Mapping the local depth of these discontinuities provides insight into the temperature and compositional state of the Earth’s interior, the locations of sinking slabs or a mantle upwelling, giving a glimpse into the evolution of our planet. Examples of such discontinuities are the Moho, that separates the crust from the mantle and the core–mantle boundary, which signals the separation of the mantle from the liquid iron core. More subtle, but still well observed, discontinuities exist in the middle of the mantle, located roughly at depths of 410 and 660km (Niazi & Anderson 1965; Engdahl & Flinn 1969; Whitcomb & Anderson 1970; Shearer 1991, 2000). These boundaries define the transition zone that separates the upper from the lower mantle. The jumps in wave speeds and density observed in this zone are believed to be mainly related to the phase transition from olivine to wadsleyite for the 410-discontinuity and from ringwoodite to bridgmanite and ferropericlase for the 660-discontinuity (Ringwood 1975; Akaogi et al. 1989; Katsura & Ito 1989; Shearer 1991; Ita & Stixrude 2017). The depth of the boundaries can be used as a proxy for local temperature and/or composition, and can therefore provide important information on mantle dynamics.
Various seismological methods used to observe these boundaries exploit the properties of the interaction of the wavefield with the discontinuity. When a travelling wave encounters a discontinuity, part of its energy can reflect or convert to another wave type. For the deeper mantle discontinuities, reflections are widely used in studies of SS and PP precursors (Whitcomb & Anderson 1970; Estabrook & Kind 1996; Deuss 2009; Koroni & Trampert 2021), but wave conversions are often used as well (Vinnik 1977; Langston 1979; Paulssen 1985; Lawrence & Shearer 2006; Rondenay 2009; Gao & Liu 2014; Van Stiphout et al. 2019).
In this work, we focus on one of these latter methods, namely receiver functions. Receiver functions use P waves that convert to S waves or vice versa. When a wave reaches a discontinuity, energy can be split, such that an incoming P wave will result in a transmitted P and SV waves. The weak signal of the converted SV wave can be more easily observed using a receiver function. In this study, we define a receiver function as a deconvolution of one spatial component of the incoming wave with another. We use P-to-SV converted waves measured in the radial-vertical plane, although other configurations are possible (Rondenay 2009). For P-to-S conversions, the radial component is deconvolved with the vertical:
Here |$F^{-1}\left[\hat{u}_i\right]$| denotes the inverse Fourier transform of |$\hat{u}_i(\mathbf {x}_r,\omega )$| (where i = R, V, T denotes the radial, vertical or transverse component and ω represents the frequency) and we use the fact that a deconvolution in the time domain corresponds to a division in the frequency domain. We follow Langston (1979), and consider a seismogram as a convolution of multiple signals:
where S(xs, t) is the source time function at the source location xs, Pi(xs, x, t) is a propagator term that describes the signal acquired between the source, xs, and x, an arbitrary location below the discontinuity. Ei(x, xr, t) is the near-receiver Earth response containing the signal from x upwards to the receiver location xr, and I(t) is the instrument response. In this notation, RF(xr, t) is a time-series measured at the receiver location xr, but implicitly depends on the structure between xs and xr. The far-earth contribution, Pi, is generally not considered (Langston 1979), as the propagation for the V and R components is thought to be similar, resulting in the assumption that PR(xs, x, t) ≈ PV(xs, x, t). This reduces the receiver function to:
with the assumption in eq. (3), RF(xr, t) now only implicitly depends on the structure between x and xr.
Many methods have been employed to project receiver function observations to locations on the discontinuity in the subsurface, such as the ray-theoretical backprojection of receiver functions stacked by slowness (Kind et al. 2002) or backazimuth (Kosarev et al. 1993). Other approaches, inspired by techniques used in reflection seismics such as a normal moveout correction (Chen & Niu 2013) or common conversion point stacking (Dueker & Sheehan 1997), have been applied as well. Since the start of this millennium, migration of receiver functions or (waveforms of) scattered waves are used to image the discontinuities of the subsurface (Ryberg & Weber 2000; Bostock et al. 2001; Poppeliers & Pavlis 2003; Wilson & Aster 2005; Liu & Levander 2013; Cheng et al. 2016). In addition, methods using scattering kernels to image receiver function observations have been developed, whereby the forward wavefield is calculated either using ray theory (Hansen & Schmandt 2017) or spectral elements (Harmon et al. 2022) assuming a single incoming wave. Nearly all of the aforementioned methods assume, at least in part, ray-theoretical wave propagation and an incoming plane wave. To reduce the amount of simplifying assumptions, various techniques employ full waveform inversion methods where waveforms are modelled and compared to real observations. Hybrid techniques have been used by Monteiller et al. (2013) and Tong et al. (2014) to reduce the high computational cost of full waveform modelling. Only part of the wave propagation (usually near the receiver, i.e. Ei(x, xr, t)) is then calculated to a high frequency using a complex, 3-D model of the region. This significantly reduces the computational cost of complete full waveform calculations. However, a plane incoming wavefront is still sometimes assumed (Tong et al. 2014). More recent developments consider a more complex, multiphase incoming wavefield either by inserting multiple phases explicitly (Wang et al. 2021a), or using AxiSEM (Nissan-Meyer et al. 2014) for global wave propagation outside the high-resolution region (Beller et al. 2018; Pienkowska et al. 2020). However, the hybrid modelling always assumes that no perturbations occur outside of the considered high-resolution area.
In previous work (De Jong et al. 2022), we have shown that this last assumption needs to be treated with care. We found that a receiver function waveform is sensitive to a substantial part of medium the wavefield travels through. We also noted the limitations of ray theory. Here, we will investigate the waveform sensitivity kernels of receiver functions using the full physics of a realistic teleseismic receiver function applied to the converted phases of the mantle transition zone: P660s and P410s. We determine what precisely contributes to their waveform and to what extend they are actually sensitive to the boundary topography, as well as what their sensitivity to P- and S-wave velocities and impedance throughout the mantle is. With the complete physics of a realistic problem, we investigate the validity of the commonly used assumptions and at what cost they can be relaxed.
2 METHODS
2.1 Adjoint method and sensitivity kernels
We use the adjoint method with full waveform calculations to investigate teleseismic receiver functions sensitivity to 3-D structural parameters throughout the mantle and to boundary topography (Tarantola 1984; Tromp et al. 2005; Fichtner et al. 2006). To this end, we calculate the kernels based on a least-squares waveform misfit:
In previous work (De Jong et al. 2022), we have derived the corresponding adjoint source:
where F−1 denotes the inverse Fourier transform, |$\Delta \hat{\mathrm{ RF}}$| the difference between synthetic and observed receiver functions at xr in the frequency domain, and w(t) the time window of observation. † signifies the parameter being in the adjoint space, which propagates backward in time. Hence, it depends on T − t rather than t, where T corresponds to the final time stamp. The adjoint wavefield interacts with the forward field (u(x, t)) to give the sensitivity kernels (Tromp et al. 2005):
where u(xr, t) denotes the forward field and u†(xr, T − t) the time-reversed adjoint wavefield, D(x) constitutes the deviatoric strain tensor and ρ, κ, μ represent the structural parameters of density, bulk modulus and shear modulus, respectively. In the following, we will show the sensitivity to a linear combination of these kernels, where α refers to the P-wave speed, β to the S-wave speed and ρ′ is equivalent to impedance (Zhu et al. 2009):
These kernels represent the derivatives of the misfit function with respect to relative perturbations in the volumetric structural parameters.
Beside the volumetric parameters in the mantle, we also consider the sensitivity to changes in boundary topography (the 410- and 660-discontinuity). Boundary sensitivity kernels for a solid–solid boundary are given by Tromp et al. (2005):
where c and ϵ represent the elastic and strain tensor, respectively, |$\hat{\mathbf {n}}$| the surface normal and ∂n the normal derivative given by |$\hat{\mathbf {n}}\cdot \nabla$|. The total misfit perturbation can then be expressed as:
and will guide the optimization workflow to infer structural volumetric and boundary topography parameters.
2.2 Practical implementations
Deconvolutions have inherent instability when the denominator in the frequency domain becomes too small. To stabilize this division for both receiver functions and their adjoint sources, we use the deconvolution method as specified by Langston (1979), where both the numerator and denominator are multiplied with the complex conjugate of the denominator and a small value, ϵ, is added to the denominator as a water level. An additional Gaussian low-pass filter is applied, such that no high frequencies are introduced by the deconvolution following (Langston 1979):
where we set ωa = 0.05 Hz. Our receiver functions are thus estimated as:
The adjoint source is also stabilized using a water-level, ϵ:
In both eqs (11) and (12), ϵ is set to 1 per cent of the maximum value of the denominator.
2.3 Synthetic models
To investigate the sensitivity of teleseismic receiver functions to 3-D structural mantle parameters and topography of the major transition zone discontinuities, we calculated the full synthetic wavefield using SPECFEM3D-GLOBE (Komatitsch & Tromp 2002a, b; Liu & Tromp 2008). The parameters in SPECFEM3D-GLOBE are set so that all wavefields are calculated to a maximum frequency of approximately 0.2 Hz. For details regarding our computational set-up, please see the Supporting Information. We filtered all waveforms using a fourth-order high-pass Butterworth filter at 30 s and additionally used a source-time function with a 6.7 s half-time, which acts as a low-pass filter on the seismograms. The source-time function is implemented as a Gaussian filter in SPECFEM3D-GLOBE that slowly decays to 0.2 Hz. The complete frequencies content of our waveforms thus ranges between 0.01 and 0.2 Hz, with most energy between 0.033 and 0.067 Hz. An additional Gaussian low-pass filter is applied to the receiver functions (eq. 10). The low-pass filter in the receiver functions is set to a slightly lower corner frequency than the source-time filter to eliminate higher frequencies that might be introduced by the deconvolution. We have used in total four synthetic models: (i) PREM (Dziewonski & Anderson 1981), (ii) PREM with added topography (Meier et al. 2009) on the 660- and the 410-discontinuity (‘Topo’), (iii) a model with no topography on the discontinuities but model S20RTS (Ritsema et al. 2000) with scaled dln Vp and dln ρ throughout the mantle (‘Velo’) and (iv) a model with both topography on the discontinuity and S20RTS in the mantle (‘ToVe’). Attenuation was not included in any of the models. Its effect was investigated in a separate test and was found to be negligible. Fig. 1 shows the topography on the 660- and 410-discontinuity. The figure also shows the source and receiver locations. All data are synthetic, but realistic. We use a real teleseismic event that occurred in 2010 on the Russia-China border at a depth of 578 km with a magnitude of 6.9. Our 227 stations are based on the station distribution of the Alaska Region Network (AK) to emulate a realistic data availability.
![Topography models from Meier et al. (2009) of the (a) and (c) 410- and (b) and (d) 660-discontinuity. (a) and (b) show the whole discontinuity, together with the source (beach-ball) and two stations (blue triangles), ATKA in the west and BESE in the east. In (c) and (d), we show the area of the discontinuity beneath the receivers [depicted in (a) and (b) by the black box], and the other receivers in black.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/gji/235/1/10.1093_gji_ggad274/1/m_ggad274fig1.jpeg?Expires=1747915227&Signature=W8llTBIofSQtGLTJdx~hvU9VZiBgg-J-1KEXfK3QPtZwcQcB8xg8RlZZLtGEBfRtaOnwcet-8Ag5cq5Gnd1iEZUodl6uFpGsT62WoAidME8ei35ae~EvgM4Jy7RdY0sD~C~ugEfCvfHhJ9H2GBud1MXGKfimzX3Yqe6SgJez9KC1ZlvFHmQdIST9Zt2M77~nQGpRUbT8voYx6RVsTQ8HBk8Q-a67vHSX4UrOnYT8E7EqpG0k6WVFOTpHSYgc-uTqMm2k6~pjh9tfKtuFAqt3zmexMUkW1wnsStE6pMeRT3dM5VKD7ejST1MxnWLa5h9ygEo0Sn~KGVYncZFE-2TDZw__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Topography models from Meier et al. (2009) of the (a) and (c) 410- and (b) and (d) 660-discontinuity. (a) and (b) show the whole discontinuity, together with the source (beach-ball) and two stations (blue triangles), ATKA in the west and BESE in the east. In (c) and (d), we show the area of the discontinuity beneath the receivers [depicted in (a) and (b) by the black box], and the other receivers in black.
3 RESULTS
3.1 Single source–receiver kernels: sensitivity to mantle parameters
We start by showing various single source–receiver sensitivity kernels of the mantle P- and S-wave velocities, and impedance. Initially, we use two stations (shown in Figs 1a and b), located at epicentral distances of 38° (ATKA) and 57° (BESE) from the source. Fig. 2 shows the synthetic seismograms, receiver functions and adjoint sources for the reference (PREM) and the topography-only model (Topo). The arrival times of various phases (e.g. P, P410s, PcP and P660s) are calculated using TauP (Crotwell et al. 1999) and depicted by the vertical lines in the traces. Around the arrival of the P410s- and the P660s-phases, a 15 s wide time window has been selected that is used for the respective adjoint sources. Note that this window is large enough to analyse information from waveforms with frequencies between 5 and 15 s. We also experimented with longer time windows and found that the kernels looked essentially the same, but with effects from other phases that arrive within the larger time windows.
From top to bottom: seismograms (uV(xr, t), green and uR(xr, t), magenta) for receivers at epicentral distances of 38° (ATKA) and 57° (BESE). Next: receiver functions (RF(xr, t)) for PREM (solid) and the Topo model (dashed). The differences between these are hard to detect by eye. The middle row shows the difference between receiver functions of PREM and Topo (ΔRF(xr, t)). In the fourth row, the full adjoint source (|$f^{\dagger }_i(\mathbf {x}_r,t)$|) and at the bottom, the windowed adjoint sources, where green and magenta again correspond to the vertical and radial components are shown. Two 15 s long time-windows are shown (red, dotted lines, around the P410s and P660s arrivals). In the adjoint calculations, only one window at a time was applied, depending on whether the sensitivity of the P410s- or the P660s-phase was calculated.
We use the adjoint sources given in eq. (12) and shown on the bottom two panels of Fig. 2, for the P660s-phase time window, to calculate the sensitivity of receiver functions’ P660s waveforms to the mantle parameters. Fig. 3 shows this sensitivity for the two single source–receiver pairs to α (P-wave speed), β (S-wave speed) and ρ′ (impedance). The kernels in Fig. 3 illustrate several features, namely (i) little overall sensitivity to impedance (ρ′) aside from some near-receiver surface reflections and source effects, (ii) significant sensitivity to the P660s Fresnel zone in the α-kernel prior to conversion and after conversion in the β-kernel and (iii) strong sensitivity to the direct P-wave scatterers, particularly in the α-kernels, but also faintly in the β- and ρ′-kernels (i.e. the outer ‘halo’). Additionally, sensitivity to near-surface reverberations can be observed in the α- and β-kernels, predominantly for the 38° source receiver pair (Fig. 3, top). In contrast to our previous 2-D, explosive source kernels (De Jong et al. 2022), the β-kernel demonstrates significant sensitivity to the direct S wave near the source, which arises from the fact that a realistic moment tensor is used rather than an explosive source. The 57|$^{\circ }\, \alpha$|-kernel (Fig. 3, bottom), shows some sensitivity to the scatterers of the PcP-phase, which arrives 24 s before P660s, while the 38|$^{\circ }\, \alpha$|-kernel (top) seems to indicate weak sensitivity to the upper mantle and the transition zone, possibly related to a surface reflection that reflected downwards on the 410-discontinuity.
![Receiver function sensitivity kernels of P660s for the mantle parameters α (Vp, left), β (Vs, middle) and ρ′ (impedance, right), for the two source–receiver pairs of Fig. 2 (Topo model with PREM as reference). The units of the volumetric kernels are [s−1 m−3]. The upper sensitivity kernels are for a source–receiver pair with an epicentral distance of 38°; and the bottom kernels have an epicentral distance of 57°. Ray paths are depicted in black (P), red (P660s) and blue (PcP).](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/gji/235/1/10.1093_gji_ggad274/1/m_ggad274fig3.jpeg?Expires=1747915227&Signature=VRO2gUrDA05DLQLHmrtNWTVi4zQUhAWugVvJuzNjGBp~37j-Lk0HwY548OnsUKA1S4GiKR0PQnYCxgSVX39GkiJZu0JBS6h5XIoe0dm0wcJhzD2W2bs8M6PpgmXtPBNCEPXiCwSxcHPSbewNKtIW~8vXNf~GsOOn9VPVAOzpGy4Ec6BMcHPVFoeK9hQDG7q2dfNIvxr1sYoQ-2fEvmnO2sOrtNKsGV7ypvhyEicT7XxpitgpsdiziPjMOb1vR~LF1JQLYzCL3gOw6zDClhBegze1HNLlq2jRqJG5fcUpXIF0RlVlAemWPjnoiwyPbuK74YM1sMLuHlPp43DxCFImEw__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Receiver function sensitivity kernels of P660s for the mantle parameters α (Vp, left), β (Vs, middle) and ρ′ (impedance, right), for the two source–receiver pairs of Fig. 2 (Topo model with PREM as reference). The units of the volumetric kernels are [s−1 m−3]. The upper sensitivity kernels are for a source–receiver pair with an epicentral distance of 38°; and the bottom kernels have an epicentral distance of 57°. Ray paths are depicted in black (P), red (P660s) and blue (PcP).
The mantle sensitivity kernels for the P410s phase, shown in Fig. 4, are similar to the P660s kernels of Fig. 3. In particular, for the 38° source receiver pair (top), the same sensitivity to the scatterers of the direct P and S waves is observed, as well as to the Fresnel zone of the P410s before (α-kernel) and after (β-kernel) conversion. The bottom kernels (57°), however, predominantly show α-sensitivity to the Fresnel zone of the PcP-phase (blue area surrounding the PcP-ray path). Fig. 2 shows that this phase arrives 3 s after P410s, well within the selected time window. There is relatively little sensitivity to the Fresnel zone of the P410s-phase in the 57°α-kernel. However, there is sensitivity to scattering in the near-receiver region of the β-kernel corresponding to scatterers of the P-to-S converted wave.
Same as Fig. 3, but here the P410s arrival time window is used for the adjoint source. The ray paths are depicted in black (P), purple (P410s) and blue (PcP).
3.2 Boundary topography kernels
Receiver functions are generally used to image local topography on discontinuities. To this end, we also calculated the boundary kernels as specified in the method section. These kernels depict a receiver function sensitivity to topography on a discontinuity and are therefore of particular interest. Fig. 5(a) depicts two single source–receiver boundary sensitivity kernels (660-discontinuity) using the P660s-adjoint sources shown in Fig. 2. Both 660-boundary kernels demonstrate strong sensitivity to the area surrounding the ray-theoretical conversion point, with a diameter of roughly 500 km. The topography of the model (Fig. 1) differs in sign beneath the two stations, a feature observed in the opposite sign of the sensitivity value on the kernel at the ray-theoretical conversion point locations, one indicating an elevation of the discontinuity at the conversion point (Fig. 5a, left) and the other a depression (Fig. 5a, right). Like in the mantle kernels, the boundary kernels also display sensitivity to the scatterers of the direct P wave (the outer ‘halo’) and S waves (surrounding the source area). The 38°-kernel (Fig. 5a, left) illustrates sensitivity to the area of the discontinuity between source and receiver, albeit relatively weak.
![Sensitivity to topography of the receiver functions generated by the model Topo with PREM as a reference on (a) the 660-discontinuity for the two P660s-adjoint sources shown in Fig. 2 (at 38°, left and at 57°, right) and on (b) the 410-discontinuity. The units for the boundary kernels are [s−1 m−2]. The segment of the discontinuities shown is 120° wide (between longitudes of 120° in the east to −120° in the west) and 50° tall (latitudes of 35° to 85°). The source is located on the left and the receiver on the right-hand side of the segment.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/gji/235/1/10.1093_gji_ggad274/1/m_ggad274fig5.jpeg?Expires=1747915227&Signature=x~IA0ZA6O~7obOcfjI8i5RuwGY-E92POTWmmYmx~22S0-EWqhYUG39vCOF3OHQlWgtktqN-urDQx2xOo0g7dkDjampINPs8u6LtXZNg-1CO31deUi87Y5aZdtzuEWH4~AKDSLMRAY2Dk6ztTqlvkeyjknJpFy~jNuMbehZH32PMoVy4SEb5-9JnqYjdbiTe3w55Ij0MjtDM1RsbAjGh9rs~3aCwJIOoxfxBX-J1GB6l7w6~UDs9wqL0zlFA27rehAO~YXrKCrWkjOAlIA9keJQNZQfVOlQ-NUnQv98QB6V4OL0USa0wNBWbDlgHB9NAzy-cZH4doFwP4maWkHRBKgg__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Sensitivity to topography of the receiver functions generated by the model Topo with PREM as a reference on (a) the 660-discontinuity for the two P660s-adjoint sources shown in Fig. 2 (at 38°, left and at 57°, right) and on (b) the 410-discontinuity. The units for the boundary kernels are [s−1 m−2]. The segment of the discontinuities shown is 120° wide (between longitudes of 120° in the east to −120° in the west) and 50° tall (latitudes of 35° to 85°). The source is located on the left and the receiver on the right-hand side of the segment.
The kernels for the 410-discontinuity, using the P410s-adjoint source of Fig. 5(b) demonstrate many of the same features as the 660-kernels: strong sensitivity to the conversion point area (with a diameter of roughly 350 km), sensitivity to the scatterers of the direct P wave and to scatterers surrounding the source area. In the left-kernel (38°), we again observe some sensitivity to the area between the source and receiver.
In the Supporting Information, we provide additional volumetric and boundary kernels for single source–receiver pairs systematically spanning a range of epicentral distances.
3.3 Boundary topography kernels, multiple receivers
We also calculated the sensitivity for the complete receiver array. The resulting kernels contain the sensitivity of the entire wavefield measured at all receiver array in the given time window (15 s around the ray-theoretical Ps-arrival). In Fig. 6, we show the sensitivity to boundary topography of the same deep event used previously, using all stations shown in Figs 1(c) and (d). As in the single source–receiver kernels of Fig. 5(a), the sign of the sensitivity depends on the local topography. On the 410-discontinuity, we predominantly observe negative sensitivity values beneath the receiver array (indicating elevation) with some positive values (depression) towards the east, in accordance with the local topography as shown in Fig. 1(c). The 660-sensitivity kernel mostly suggest elevation in the west and depression in the east, generally in agreement with the topography shown in Fig. 1(d). This will be further elaborated on Section 4.
![Boundary sensitivity kernels [s−1 m−2] of receiver functions generated by the model Topo at the 410-discontinuity with P410s (left) and the 660-discontinuity with the P660s (right) for the 578 km deep event using the complete receiver array. The kernels are shown on the same segment as in Fig. 5.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/gji/235/1/10.1093_gji_ggad274/1/m_ggad274fig6.jpeg?Expires=1747915227&Signature=30TSgvXRcVfHxSaIPSVSVhbD1SyEcC~Fw6gZEL0q5YQkqQUEku~uKxRhg6lBQcRXRiCNRXgIKe2I3J~4b-NlK~4UdvpxDAlnYD3CGyMU7RlbTjjmMNsv-pecyATo8iToOWtcEdf5RvWTvj6Snq1tgvGEOllrRZLRzNEOqKs80I2bdXKMcvNuoVUtp-aIiviECoOn4Jph9fdh-EVKRh2LIeORWZ1TubZpTvUjq683J3ZIS~K7ZJ6YNHfU3CIK-NDVkQP3cI7ViIjkDn2XgTwzFCfiHHrq~2JvwNHMRSPpeSWseksphq53SrsWNTn8xBRB4ACOKgNOBTLaqGxkWNIYYw__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Boundary sensitivity kernels [s−1 m−2] of receiver functions generated by the model Topo at the 410-discontinuity with P410s (left) and the 660-discontinuity with the P660s (right) for the 578 km deep event using the complete receiver array. The kernels are shown on the same segment as in Fig. 5.
This kernel also demonstrates that the amplitude of the sensitivity at a certain location does not just depend on the misfit between the data and reference models, but also on the array set-up. For the 660-discontinuity kernel, we see a strong sensitivity in an area where the topography is relative small, but where the data coverage is large. For this particular example, that indicates that the array set-up may lead to stronger sensitivity to a small change in this area than to larger changes in outer, less well covered areas of the discontinuity. Additionally, the sensitivity to source area and to the scatterers of the direct P and S waves remain visible.
3.4 Additional models
As mentioned before, we also calculated the synthetic wavefields for two additional models, namely one model with no topography on the discontinuities, but with smoothly varying mantle velocities from S20RTS (Ritsema et al. 2000), and another with both topography on the 410- and 660-discontinuities and with S20RTS in the mantle. In Fig. 7, we show the corresponding traces, receiver functions called ToVe (referring to the model with both topography and velocity perturbations), and Velo (velocity-only). The bottom three graphs show ΔRF and the (windowed) adjoint sources for a case where we consider the synthetic data model (ToVe-Pr) with PREM as a reference and one where we use the velocity-only model as a reference (ToVe-Ve), see Table 1 for an overview of the naming convention. The single source, multiple receiver 660-boundary kernels for the four adjoint calculations are shown in Fig. 8, together with the sum of the ToVe-Ve and Velo-Pr kernels. This will be further elaborated on in Section 4.
Same event and stations as in Fig. 2. Top: seismograms (ui(xr, t)). Second and middle row: (RF(xr, t)) of the models: PREM, Velo and ToVe and the differences between receiver functions, ΔRF(xr, t) (ToVe w.r.t. Velo and ToVe w.r.t. PREM). Bottom rows: the full and windowed adjoint sources (|$f^{\dagger }_i(\mathbf {x}_r,t)$|) for the P660s-arrival.
![Boundary sensitivity kernels [s−1 m−2] on the 660-discontinuity using the P660s arrival, shown on a 2-D projection of the receiver area (same region as Figs 1c and d). We show the kernels for the Topo-Pr (as in Fig. 6), Velo-Pr, ToVe-Pr and ToVe-Ve adjoint runs. The bottom-right kernel is the sum of ToVe-Ve and Velo-Pr.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/gji/235/1/10.1093_gji_ggad274/1/m_ggad274fig8.jpeg?Expires=1747915227&Signature=oxBAeX7cai02DBqPpMM5g975QbJH4NtbM3l3LI6ue2A0~Myg1tq3XVa95pVsScZn4C5xk9GXu6JrO4FH61cxuST971urtJJ9CS3F8HyyvkBfZa0Vum4a8GPPbvRMuV94rMHFwW~xFq-o1iqfAQIpS0wZCGtbz4-bCiRvPSUHOKIb9Safn5ApLj1ASJsxZuRidNE0wPq90StEqkDlxidPbqrAA9frT9iTBCKn7OAjohjCxKLcqmcTrECdGAvzKZsRQcrI2vsrP7Gd1FWO9380dag9S3vo6uPzrJfikpAqdO37vqBGLH~nvMHl-MJ75gskit3~UpKamvkWTfg4rq0htA__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Boundary sensitivity kernels [s−1 m−2] on the 660-discontinuity using the P660s arrival, shown on a 2-D projection of the receiver area (same region as Figs 1c and d). We show the kernels for the Topo-Pr (as in Fig. 6), Velo-Pr, ToVe-Pr and ToVe-Ve adjoint runs. The bottom-right kernel is the sum of ToVe-Ve and Velo-Pr.
The models and runs used in this research.
| Model name | Description | Reference model (if applicable) |
|---|---|---|
| PREM | 1-D variation, no topography, no velocity | n/a |
| Topo | Topography on discontinuities | n/a |
| Velo | 3-D mantle velocity model (S20RTS) | n/a |
| ToVe | Topography and 3-D velocity (S20RTS) | n/a |
| Topo-Pr | Topo adjoint run | PREM |
| Velo-Pr | Velo adjoint run | PREM |
| ToVe-Pr | ToVe adjoint run | PREM |
| ToVe-Ve | ToVe adjoint run | Velo |
| Model name | Description | Reference model (if applicable) |
|---|---|---|
| PREM | 1-D variation, no topography, no velocity | n/a |
| Topo | Topography on discontinuities | n/a |
| Velo | 3-D mantle velocity model (S20RTS) | n/a |
| ToVe | Topography and 3-D velocity (S20RTS) | n/a |
| Topo-Pr | Topo adjoint run | PREM |
| Velo-Pr | Velo adjoint run | PREM |
| ToVe-Pr | ToVe adjoint run | PREM |
| ToVe-Ve | ToVe adjoint run | Velo |
The models and runs used in this research.
| Model name | Description | Reference model (if applicable) |
|---|---|---|
| PREM | 1-D variation, no topography, no velocity | n/a |
| Topo | Topography on discontinuities | n/a |
| Velo | 3-D mantle velocity model (S20RTS) | n/a |
| ToVe | Topography and 3-D velocity (S20RTS) | n/a |
| Topo-Pr | Topo adjoint run | PREM |
| Velo-Pr | Velo adjoint run | PREM |
| ToVe-Pr | ToVe adjoint run | PREM |
| ToVe-Ve | ToVe adjoint run | Velo |
| Model name | Description | Reference model (if applicable) |
|---|---|---|
| PREM | 1-D variation, no topography, no velocity | n/a |
| Topo | Topography on discontinuities | n/a |
| Velo | 3-D mantle velocity model (S20RTS) | n/a |
| ToVe | Topography and 3-D velocity (S20RTS) | n/a |
| Topo-Pr | Topo adjoint run | PREM |
| Velo-Pr | Velo adjoint run | PREM |
| ToVe-Pr | ToVe adjoint run | PREM |
| ToVe-Ve | ToVe adjoint run | Velo |
4 DISCUSSION
These sensitivity kernels provide us with valuable insight into what contributes to the observations made using teleseismic receiver function waveforms.
4.1 General implications
The mantle kernels, shown in Figs 3 and 4 illustrate the sensitivity to mantle parameters and they clearly indicate that a receiver function sensitivity is not restricted to the discontinuities and near-receiver area alone, but to significant parts of the mantle. Sensitivity to the S-wave speed is largely restricted to the upper mantle, in particular to the Fresnel zone of the converted S-wave part of the P660s/P410s phases and to the scatterers of the converted wave. That said, S-wave sensitivity to the scatterers of the direct S wave, near the source, and the direct P wave cannot be ignored entirely. However the most widespread sensitivity is observed to the P-wave speed (α). In these kernels, we see sensitivity to the Fresnel zone prior to conversion, to the direct P wave, as well as several reverberations near the surface. The assumption of PR(xs, x, t) ≈ PV(xs, x, t) which implies that the contribution of the far earth cancels out by the deconvolution is not supported by our sensitivity kernels. We observe little sensitivity to the impedance, and what we observe appears mostly restricted to reflective phases near the surface.
Another important feature of the sensitivity to mantle parameters, demonstrated primarily by the α-kernels, is the significant sensitivity to other phases. The bottom α-kernel in Fig. 4 shows a significantly stronger sensitivity to the PcP-Fresnel zone than to the P410s, which happens to arrive within the time window. But even when the other phase does not arrive in our considered time window, its scatterers might still contribute to the receiver function’s shape, as can be observed in the bottom α-kernel of Fig. 3. There, we observe sensitivity to the PcP-scatterers even though the PcP-phase arrives significantly before our time window. Additionally, all kernels clearly show significant sensitivity to scatterers of the strong direct P and S waves.
In short, our results demonstrate that perturbations of the reference mantle velocity model can partially explain the observations made by receiver functions. These velocity perturbations may be located before or after conversion, may be related to P- or S-wave speed, and can occur within the Fresnel zone of the converted phase as well as outside of it. In other words, one needs to have an accurate velocity model before a receiver function waveform can be used to infer topography of the considered discontinuity.
On the other hand, the boundary kernels in Figs 5 and 6 show that receiver functions have significant sensitivity directly to discontinuity topography as well, in particular near the conversion point. The sensitivity to boundary topography, shown in Fig. 5, is strongest in the region surrounding the ray-theoretical conversion point, although there is some sensitivity to the source area as well. The sign of the sensitivity value that we observe in the kernels of Fig. 5 and 6 mostly corresponds to elevations and depressions in the input topography model indicating that the boundary sensitivity kernels, actually ‘see’ the topography on the respective discontinuities.
So far, we can draw two main conclusions: (i) the receiver functions are sensitive to the velocity field through large parts of the mantle and, apart from the Ps waves, may be influenced by many other phases and their scatterers, and (ii) receiver functions show clear sensitivity to the discontinuity topography near the conversion point. Although we have demonstrated that the assumptions of ray theory and a known incoming wave front commonly used in receiver function studies are not supported by our results, we also show that receiver functions are sensitive to the boundary topography information that they are generally used for. This sensitivity is highly localized around the conversion point area, especially compared to the boundary sensitivity of PP- or SS-underside reflections (Koroni et al. 2019; Koroni & Trampert 2016), that display a wide-spread X-shaped sensitivity due to their minimax nature.
The kernels in Fig. 8 allow us to further investigate the interaction of the sensitivity between the boundary topography and the mantle velocity model. In this figure, we compare several kernels corresponding to varying data and reference models. In the top left, we show the kernel of the adjoint run used thus far (Topo-Pr), with on the right the kernel of the adjoint run using only velocity perturbations (Velo-Pr). The top-right kernel of this figure shows the sensitivity kernel of the receiver functions generated with a model containing both topography and velocity perturbations and using PREM as a reference model (ToVe-Pr). It clearly differs from the Topo-Pr kernel, particularly in the north, as a consequence of the velocity perturbations. However, when using the S20RTS-velocity model as a reference, as done for the kernel on the bottom-left (ToVe-Ve), we get a result that is almost identical to the Topo-Pr sensitivity kernel, implying that:
This relation shows that when the correct velocity model is used, regardless of what it is, the boundary kernel will show the sensitivity caused by topography perturbations that is required to update the boundary topography model. In this hypothetical case, we could invert for topography only and use the boundary kernel to update our model, while leaving the elastic parameters untouched. The opposite is also true, however. The ToVe-Pr boundary kernel (top right in Fig. 8) shows that when the incorrect velocity model is used, an additional misfit change is projected onto the discontinuity sensitivity. The top- and bottom-right kernels in Fig. 8 illustrate the following relation:
This implies that the topography will be updated incorrectly when the elastic parameters are ignored in an inversion with a reference velocity model as a reference that deviates significantly from the actual velocity structure. Because the receiver function sensitivity kernels cannot a priori ‘see’ which perturbations arise from velocity and which arise from topography, it only ‘sees’ where they could have originated. As such, it will project the velocity perturbation onto the discontinuity and vice versa, which will lead to an inaccurate result if not accounted for. Here, we show that when the observations are solely attributed to topography, you need to be sufficiently certain of the velocity model in the entire mantle or your results will be incorrect. If you suspect your reference model to be incorrect, a simultaneous inversion of both topography and velocity is advised.
4.2 Topography sensitivity and mapping techniques
We further investigated how our results relate to commonly used methods of mapping receiver function observations onto discontinuity topography. It is clear that mapping the observation to a single, ray-theoretical conversion point is both a strong approximation, and impractical when receiver functions are stacked, as they often are. Therefore, the observations are generally mapped to some area beneath a single receiver (Langston 1979) or when the array is dense enough, to regions such as common conversion point bins (Dueker & Sheehan 1997), which vary in size, but often range somewhere between 50–200 km (Dueker & Sheehan 1997; Niu et al. 2004; Van Stiphout et al. 2019).
Lekic et al. (2011) developed an approach that maps receiver function observations into bins using a weighting factor that depends on the distance to the ray-theoretical conversion point. The cubic weighting function decreases from 1 at the ray-theoretical conversion point to 0 at the Fresnel zone half-width taken at half the dominant S-wave wavelength. In Fig. 9, we compare the spread of Lekic’s weighting function (Fig. 9a) to our 660-discontinuity kernel (Fig. 9b). It is important to note that while sensitivity kernels and weighting functions are physically two completely different quantities, both express what part of the discontinuity is assumed to contribute to the observations. It is in that light that we are comparing them. Comparing the size of the weighted area in Fig. 9(a) to the sensitivity displayed in Fig. 9(b), we observe that the weighted area roughly corresponds to the inner circle in the kernel (blue area). This is better visualized in Fig. 9(e), where we compare the weight to the sensitivity values across the source receiver great-circle line. There we observe that the weighted area proposed by Lekic et al. (2011) corresponds to the area surrounding the conversion point that the receiver function is most sensitive to.
![Weighting function area compared to our sensitivity kernels [s−1 m−2] at 660km conversion. (a) and (c) weighting functions from Lekic et al. (2011) for (a) a single source–receiver pair (BESE) and (c) the full receiver array. (b) and (d) boundary sensitivity kernel for (b) the single source–receiver pair, and (d) the full receiver array. (e) Sensitivity (right y-axis) and weighting function (left y-axis) along the source–receiver great circle.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/gji/235/1/10.1093_gji_ggad274/1/m_ggad274fig9.jpeg?Expires=1747915227&Signature=ynwODYFl7eSlMjMu6QUDuBE~QGNgIIcOVa3QOzQpRCjecvGjGeirEROAu-QliuYcdpRdCUVgAqSgOTL7bREdy832UkeTncRF9PgEcf9r5cSCJD9yGdAzTEeIsglzajDjicNpzWYbT0xCg~xN9YAHaLhJrju008-uJ6aH3APpt9pbIE10SQueZ1M9X5NmrxTOSpvo8Dg5Ca~l77fdoZnDBzj1RKEI46Y-l7RF5sRYs48U2GGQQrR~JWtHSmrW1cHMynJwi9q4TdWPsWjsCcyT~0iOTOnzCYvLntWjNtWCOBjGZbxUQr-o3uosJUIftzVBh~9Cg4edIWTeX93eRsZDPg__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Weighting function area compared to our sensitivity kernels [s−1 m−2] at 660km conversion. (a) and (c) weighting functions from Lekic et al. (2011) for (a) a single source–receiver pair (BESE) and (c) the full receiver array. (b) and (d) boundary sensitivity kernel for (b) the single source–receiver pair, and (d) the full receiver array. (e) Sensitivity (right y-axis) and weighting function (left y-axis) along the source–receiver great circle.
In Figs 9(c) and (d), we compare the contributing area to the sensitivity kernel for multiple receivers. The most important observation here is that in the weighting method the area with the highest contribution to the observations is the region with the most data coverage, that is, the highest receiver density. However, in the sensitivity kernels, the amplitude of the sensitivity depends both on receiver density as well as on the local magnitude of the difference between the ‘observed data’ and ‘synthetic data’. This results in the relatively high sensitivity towards the edges of the array, where the receiver coverage is relatively low but the difference between the reference and data model is the highest.
When comparing our sensitivity kernels to migration (Ryberg & Weber 2000; Bostock et al. 2001), scattering kernels (Hansen & Schmandt 2017; Harmon et al. 2022) and the hybrid full waveform propagation method (Monteiller et al. 2013; Tong et al. 2014), we note three important observations: (i) our kernels show a direct sensitivity to velocity, in particular to the Fresnel zone of the converted phase, in contrast with the scattering kernels shown in the work of Hansen & Schmandt (2017), (ii) our kernels demonstrate sensitivity to velocity throughout the entire mantle and (iii) they show possible sensitivity to other, interfering phases. These last two observations are particularly relevant for the hybrid waveform modelling method. Compared to full waveform modelling and inversion, the hybrid method has the advantage of modelling high-resolution waveforms relatively fast, but it comes at the cost of information loss. Considering the full global propagation outside the region (Beller et al. 2018; Pienkowska et al. 2020), or explicitly inserting multiple phases (Wang et al. 2021a) is an important step to reduce the amount of information loss. However, the sensitivity of receiver functions to the velocity structures in the far field is not insignificant, and therefore might still influence the observations in significant ways if not properly taken into account.
4.3 Moment tensor
An important contributor to the details of the receiver function sensitivity kernels has not been addressed so far: the moment tensor. The moment tensor has a significant influence on the wavefield, its propagation and the relative amplitude of various arrivals at the receiver, and thus, inevitably, on the receiver functions and their sensitivity. The event we use in our example is realistic, but also purposefully chosen as it generates a high P-excitation towards the receivers. Fig. 10 presents two additional 660-boundary kernels, where we use events at the same location but with different moment tensors. The kernel on the left shows the sensitivity for a moment tensor that does not generate strong P-wave excitation towards the receiver, and thus fails to excite clear Ps-conversions. This results in a generally low sensitivity (a factor of 30 smaller) to boundary topography beneath the receivers for the waveforms in the chosen time window. We note that the sensitivity near the source is not reduced as that is related to scatterers of other phases (direct S-wave predominantly). The second kernel (on the right) does have a high P excitation and its sensitivity is more similar to the original kernel on Fig. 6, both in amplitude and sign.
![Boundary sensitivity [s−1 m−2] for events at the same location as before, but with two different moment tensors.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/gji/235/1/10.1093_gji_ggad274/1/m_ggad274fig10.jpeg?Expires=1747915227&Signature=g~AtwU0W2I8vsLsTsaC0EdwRT6QICVGGsuLWX-MKhkmdSRFG4DfdjDprlKFvxPZZlDO2CkvzZWPd-AyPrewR9t3F4l3Sd3wfCRgTurrkhLk8ZwWIZAOtP~0BjkN6NalUOAxoR3TWFFFdB-iLeFvHUXprp0sOnNBNhGfwPn1m3cn9CIB5RPec2VWdzJ3h-Mfw5dvMnUC1R6PAiFc3If5EE3gYvX8twuL4Bo-KAll1nIAdK0DiRBpDj7h37NYmJeWvIuEW4pfR7gC1RUX3ZpfNOJ0RnZAK~ETS7aXOPvy4AV5HXgnwdb~GjUgA-xyUaza~I2VChmbB8gfC6wQNW3CWQA__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Boundary sensitivity [s−1 m−2] for events at the same location as before, but with two different moment tensors.
4.4 Additional events
In many receiver function investigations, more than one event is used. The receiver functions of these multiple events are binned, either by station (Kosarev et al. 1993) or by common conversion point (Dueker & Sheehan 1997), and stacked to extract the relatively weak signal of the converted phase. We calculated kernels for three additional events in order to investigate how multiple events might affect a full waveform approach. We consider the sensitivity kernels for various events individually, and then directly stack the kernels, not the receiver functions. These additional earthquakes are also based on real events. They were selected to have a similar magnitude, but vary with depth, azimuth and moment tensor, see Table 2 and Fig. 11(f).
![(a)–(d) Boundary sensitivity kernels [s−1 m−2] for the four events in Table 2. (e) The sum of the four kernels. (f) Topography map (as in Fig. 1b) with the moment tensors plotted at the earthquake locations.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/gji/235/1/10.1093_gji_ggad274/1/m_ggad274fig11.jpeg?Expires=1747915227&Signature=FzigyLIq34tSohxDC4r7OQX8nsgw2Z2fVmTnYVstoFD63nBPsxyvmFi4llqfCPBpCLrczSiTJ2hHr2H~ZigHW8NkGjMiyDuoyVfXXcUJt~DwyAMAFcQWe-RgZceFufcyy1X2ZoMiOvZAqc1e7nUTSvOyIyNzvLfA2~4J9tnzeFAMsEY3b6jdRK1zz9IBTKXp~ol07AGLgSKxLS4EtT83V5VRJDL2m1qwCJVaAut~oofIcAv-jxWGE3CrcJeW726egw~eXA9MHyCmxk1x1axCuc6PNCmwBpoe5c~DBKS0t9EEMfmXb~4AqD2K8dXaFB7RBRrBdKakZ6Qw7-EWGws8mA__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
(a)–(d) Boundary sensitivity kernels [s−1 m−2] for the four events in Table 2. (e) The sum of the four kernels. (f) Topography map (as in Fig. 1b) with the moment tensors plotted at the earthquake locations.
Characteristics of four real events that our examples are based on. The first event was used as an example in the rest of this paper. The moment tensors for event II and III were used for the kernels in Fig. 10
| Event | Place | Exact location | Depth (km) | Epicentral Distance | Magnitude | Date |
|---|---|---|---|---|---|---|
| I | Russia-China | 42.5°, 130.7° | 578 | 47.7° ± 20° | 6.9 | 18-02-2010 |
| II | Mariana Islands | 13.5°, 144.5° | 140 | 65.8° ± 20° | 6.7 | 17-09-2014 |
| III | Hindu Kush | 36.5°, 70.5° | 214 | 75.7° ± 20° | 7.0 | 09-08-1993 |
| IV | Romania | 45.9°, 26.8° | 74 | 72.1° ± 20° | 7.0 | 30-05-1990 |
| Event | Place | Exact location | Depth (km) | Epicentral Distance | Magnitude | Date |
|---|---|---|---|---|---|---|
| I | Russia-China | 42.5°, 130.7° | 578 | 47.7° ± 20° | 6.9 | 18-02-2010 |
| II | Mariana Islands | 13.5°, 144.5° | 140 | 65.8° ± 20° | 6.7 | 17-09-2014 |
| III | Hindu Kush | 36.5°, 70.5° | 214 | 75.7° ± 20° | 7.0 | 09-08-1993 |
| IV | Romania | 45.9°, 26.8° | 74 | 72.1° ± 20° | 7.0 | 30-05-1990 |
Characteristics of four real events that our examples are based on. The first event was used as an example in the rest of this paper. The moment tensors for event II and III were used for the kernels in Fig. 10
| Event | Place | Exact location | Depth (km) | Epicentral Distance | Magnitude | Date |
|---|---|---|---|---|---|---|
| I | Russia-China | 42.5°, 130.7° | 578 | 47.7° ± 20° | 6.9 | 18-02-2010 |
| II | Mariana Islands | 13.5°, 144.5° | 140 | 65.8° ± 20° | 6.7 | 17-09-2014 |
| III | Hindu Kush | 36.5°, 70.5° | 214 | 75.7° ± 20° | 7.0 | 09-08-1993 |
| IV | Romania | 45.9°, 26.8° | 74 | 72.1° ± 20° | 7.0 | 30-05-1990 |
| Event | Place | Exact location | Depth (km) | Epicentral Distance | Magnitude | Date |
|---|---|---|---|---|---|---|
| I | Russia-China | 42.5°, 130.7° | 578 | 47.7° ± 20° | 6.9 | 18-02-2010 |
| II | Mariana Islands | 13.5°, 144.5° | 140 | 65.8° ± 20° | 6.7 | 17-09-2014 |
| III | Hindu Kush | 36.5°, 70.5° | 214 | 75.7° ± 20° | 7.0 | 09-08-1993 |
| IV | Romania | 45.9°, 26.8° | 74 | 72.1° ± 20° | 7.0 | 30-05-1990 |
In Figs 11(a)–(d), we show the boundary sensitivity kernels for the four individual events. We note a wide variety of sensitivity. As already mentioned in the previous section, not all moment tensors generate a strong P-excitation, and consequently strong P-to-S conversions. For events II and IV, this results in sensitivity kernels that seem weaker and less coherent, in particular beneath the receivers. The other two events (I and III) demonstrate significant sensitivity to the discontinuity under the receivers. Their sensitivity is also about an order of magnitude larger than the sensitivity of events II and IV. In the summation kernel, we therefore predominantly observe the contribution from events I and III. We only see the contribution from events II and IV near their respective source locations. The sensitivity to the region beneath the middle of the receiver array, where the local topography is the smallest, the kernels for evenst I and III differ significantly. This relates in part to sensitivity to the arrivals of other phases in the kernels of event III (see Supporting Information) and to the small topography perturbation in that area. However, in the region beneath the south and eastern edge of the receiver array, where the topography varies most with respect to PREM, the kernels are more similar.
From this brief investigation into the contribution of multiple events, we tentatively infer the following: (i) using multiple events is likely to better constrain the summed sensitivity kernels to regions where the observed data are most affected by the perturbations. (ii) The relative strength of sensitivity to the discontinuity topography differs widely between events due to variations in magnitude, epicentral distance, azimuth, depth and the moment tensor. The great advantage of using full waveform inversion sensitivity kernels is that all these effects are automatically considered, and therefore the source–receiver pairs with the strongest sensitivity are automatically weighted the most. Stacking, however, introduces significant uncertainty as stacked receiver functions are generally assumed to be the sum of more or less equal contributions from many individual receiver functions, while it is more likely that some events will contribute significantly more than others.
5 CONCLUSION
Receiver functions are generally used under the assumption that only the depth of the discontinuities and the near-receiver velocity model are contributing significantly to the observations. Assuming that the near-receiver velocity structure is reasonably well known, the observed arrivals of converted phases can accurately be backprojected to determine the local depth of a discontinuity. In addition, both a ray theoretical approximation is made and a incoming plane wave of a certain slowness is often assumed.
We investigated the validity of these assumptions by applying the adjoint method of full waveform inversion to calculate sensitivity kernels of receiver function waveforms for synthetic examples inspired by real events. We calculated both the sensitivity to topography on the boundaries and sensitivity to structural parameters in the mantle. The boundary kernels demonstrate that receiver functions are indeed strongly sensitive to the topography of the discontinuity in the conversion point area. Although sensitivity to the source region and along the wider discontinuity is observed, the strongest sensitivity is located near the receiver surrounding the conversion point, matching the weighting function of Lekic et al. (2011). When calculating sensitivity for multiple source–receiver pairs simultaneously, we note that the sign and strength of the local sensitivity depend on local topography of the discontinuity, as well as on density of the receiver array.
Additionally, we demonstrate that the calculated boundary sensitivity is independent of the velocity model, only if the velocity model is perfectly known. In that hypothetical case, the kernels only show the sensitivity to the topography perturbations and can be used to invert for topography alone. However, when the reference model is incorrect, the misfit change induced by velocity perturbations will be projected onto the discontinuity and generate biased topography results.
This is further highlighted by the mantle sensitivity kernels that show widespread sensitivity of receiver functions to the velocity parameters of the mantle. These kernels show that it is not just the near-receiver Earth structure that contributes to the observations. Instead they demonstrate sensitivity to the whole Fresnel zone of the converted P-to-S phases, including the path prior to conversion. In addition, they indicate sensitivity to scatterers of other phases such as the direct P and S waves, or the PcP-arrivals. The strength of these phases’ sensitivity may exceed the sensitivity to the P-to-S phase, especially when they arrive in the same time window. But even when the phases arrive much earlier, their scatterers might still contribute significantly to the observations.
We therefore propose that receiver functions are best incorporated into a full inversion scheme which inverts for both topography and velocity structure simultaneously. Our sensitivity kernels show that receiver functions indeed contain valuable information about the topography and local depth of discontinuities. However, without a perfectly accurate velocity model throughout the mantle, it is difficult to attribute their information content solely to the local topography on the discontinuity (Ammon et al. 1990). Nevertheless, when applied in a simultaneous inversion which incorporates other data to update the mantle velocity parameters, receiver function adjoint sources can be used to calculate boundary sensitivity kernels and update topography.
Alternatively, the sensitivity kernels could be used as a tool to assess which phases and mantle regions might be contributing to the waveform in the considered time window and to what extend, or whether the event’s moment tensor is even capable of generating receiver functions with a sensitivity to the discontinuity. All of this knowledge could prove valuable in other, less computationally expensive methods using receiver functions (Ryberg & Weber 2000; Bostock et al. 2001; Monteiller et al. 2013; Tong et al. 2014).
SUPPORTING INFORMATION
Figure S1. The topography of the 660-discontinuity with the theoretical station array and explosive source used for the kernels given in Figs 2–4.
Figure S2. Boundary sensitivity kernels [s−1 m−2] for the receiver array given in Table 3. Right: the kernel for the model with the topography shown in Fig. 1. Left: the kernel for the model with the flat, shifted discontinuity.
Figure S3. Volumetric sensitivity kernels [s−1 m−3] corresponding to the boundary kernels shown in Fig. 2. Top: the kernel for the topography model, and bottom: for the shifted 660-discontinuity.
Figure S4. (a)–(d) Volumetric [s−1 m−3] and (e)–(h) boundary [s−1 m−2] sensitivity kernels for four individual source–receiver pairs using the topography model. Note: the colour scale of the various kernels varies.
Figure S5. 660 boundary sensitivity kernels [s−1 m−2] for (a) and (b) event I and (c) and (d) event III. In (a) and (c), we show the kernels using all source–receiver pairs, while in (b) and (d) the pairs that might be contaminated by other phases were removed.
Table S1. Settings used for the simulations discussed in the main text.
Table S2. Settings used for the lower resolution simulations shown in the Supporting Information.
Table S3. Topography perturbation at the conversion points. Negative values indicate elevation.
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 paper.
ACKNOWLEDGEMENTS
This work has been financed by the research programme DeepNL of the Dutch Research Council (NWO) under project number DeepNL.2018.033. We used the open source software SPECFEM3D GLOBE (available at https://geodynamics.org/resources/specfem3d), for which we would like to thank the developers. Additionally, we would like to thank the reviewers, Daniel Peter and Qinya Liu, for their kind and constructive comments, as well as Thomas Cullison for his invaluable help with the SPECFEM3D workflow.
DATA AVAILABILITY
The synthetic data generated for this study, as well as the SPECFEM3D source files and codes are available upon request from the authors.
