# Modelling P wave reflections on MTZ discontinuities from distant oceanic sources

P Pawlowski, H Pedersen, B Tauzin, Pierre Boué, Lisa Tomasetto

► To cite this version:

P Pawlowski, H Pedersen, B Tauzin, Pierre Boué, Lisa Tomasetto. Modelling P wave reflections on MTZ discontinuities from distant oceanic sources. Geophysical Journal International, 2025, 242 (3), 10.1093/gji/ggaf275. http://dx.doi.org/10.5594/gji.12684

HAL Id: hal-05217498

https://hal.science/hal-05217498v1

Submitted on 22 Aug 2025

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
# Modelling *P* wave reflections on MTZ discontinuities from distant oceanic sources

P. Pawlowski<sup>1</sup>, H. A. Pedersen<sup>1</sup>, P. Boué,<sup>1</sup> B. Tauzin<sup>2,3</sup> and L. Tomasetto<sup>1</sup>

Accepted 2025 July 14. Received 2025 July 10; in original form 2025 May 8

## SUMMARY

We investigate the modelling of *P*-wave reflections on the mantle transition zone (MTZ) discontinuities (Pv410p* and Pv660p*) using ambient seismic noise generated by distant oceanic sources. Using ray theory and waveform simulations, we assess biases in arrival times and amplitude ratios when interpreting noise correlations as Green's functions. Our results show that source distribution and the b-caustic effect strongly influence signal recovery. Simulations based on realistic oceanic models (WAVEWATCH III) demonstrate that appropriate source conditions significantly reduce biases. This approach enables reliable imaging of the MTZ, particularly in regions like the greater Alpine area with favourable microseismic source distribution.

**Key words:** Seismic discontinuities; Seismic interferometry; Seismic noise; Wave propagation.

## 1 INTRODUCTION

The mantle transition zone (MTZ), often referred to as the boundary between the upper and lower mantle, is a region located between approximately 410 and 660 km depth in the Earth's mantle, where mineral phase changes occur due to increasing pressure. The heterogeneous composition of the MTZ has a major impact on mantle convection, making it a key element in mantle dynamics (for a review, see Goes *et al.* 2022). The use of long-period seismograms has provided global scale images of the MTZ over the last few decades using specific phases as SS and PP precursors from underside reflections (e.g. Shearer 1991), or *P*-to-*SV* converted phases receiver functions (e.g. Chevrot *et al.* 1999; Lawrence & Shearer 2006; Tauzin *et al.* 2008).

A decade ago, Poli *et al.* (2012) showed that topside reflections on the MTZ, a.k.a Pv410p and Pv660p phases (TauP phase name convention Crotwell *et al.* 1999), can be detected from the correlation of ambient noise in cratonic regions (stable continental crust). Poli *et al.* (2012) is a direct continuation of other studies that aimed to exploit microseismic noise for imaging the major interfaces of our planet, such as the Moho in the work of Ruigrok *et al.* (2011). We can also note that the microseismic noise can be used directly through the receiver function method for the detection and mapping of upper-mantle interfaces, as successfully recently demonstrated by Kato & Nishida (2023) for the MTZ beneath Japan.

Following Poli *et al.* (2012), other studies confirmed the possibility to use correlation functions in the frequency band corresponding

to the secondary microseismic peak to image the MTZ, especially in the particularly favourable cratonic region of North China (Feng *et al.* 2017, 2021). In more complex geological structures, such as collision zones, it seems less trivial to make similar observations and standard interferometric methods need to be improved, in particular through time and space selective approaches (Pedersen & Colombi 2018). Lu *et al.* (2021) showed that low H/V (horizontal over vertical) component ratios of continuous seismic noise data relative to yearly average are indicative of days with a significant amount of subvertically propagating *P* waves, in comparison to surface waves. They further demonstrated that low H/V occurred on days with significant storm activities in the Southern Hemisphere as compared to storm activities in the North Atlantic. Pedersen *et al.* (2022) showed that selecting days with low H/V significantly improved the extraction of Pv410p and Pv660p from noise correlations, consistent with the assumption that the correlation extract seismic waves from distant ocean sources. Based on this assumption, Lu *et al.* (2023) successfully mapped MTZ in the greater Alpine region using Pv410p and Pv660p extracted from seismic noise using data only from Northern Hemisphere summer months, when H/V is overall lower than during other times of year.

One of the main aims of MTZ studies is to constrain lateral variations in mantle temperatures (due to e.g. slabs and plumes). Additionally, lateral variations in mantle composition, including the effect of hydrous phases, also contribute to lateral variations of these reflectors. In turn, these variations may also influence mantle dynamics (Anderson 1987; Christensen 1995). Methods similar

<footnote n="1"><sup>1</sup>Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTERre, F-38000 Grenoble, France. E-mail: piel.pawlowski@univ-grenoble-alpes.fr</footnote>

<footnote n="2"><sup>2</sup>Mantle8, F-31100 Toulouse, France</footnote>

<footnote n="3"><sup>3</sup>Université Lyon 1, Ecole Normale Supérieure de Lyon, CNRS, LGL-TPE, F-69622 Villeurbanne, France</footnote>
to the one initially proposed by Poli *et al.* (2012) and developed
since then by, for example, Lu *et al.* (2023), seem particularly well
suited to tackling these problems, especially in combination with
the existence of large observational networks (USAarray, Alparray,
Adriarray,...). The combination of dense networks and these new
techniques is therefore very promising for the study of mantle dis-
continuities.

The previously cited observations are interpreted as if the ex-
tracted waves are parts of the Green’s function. Extracting Pv410p
and Pv660p using days with dominance of specific noise sources
(Pedersen *et al.* 2022) means moving significantly away from the
conditions necessary for the use of correlations as a reasonable
proxy for the equivalent Green’s function (Wapenaar 2004), that is,
fully distributed noise sources and/or a diffusive wavefield. Accu-
rate modelling of the mantle reflected waves using realistic sources
is, therefore, needed to understand the processes underlying the
wave reconstruction and the possible errors that are committed by
interpreting the extracted waves as part of the Green’s function
between pairs of stations.

Taking into account the noise sources in seismic tomography
with noise correlation (widely used for surface wave applications)
has been explored in the context of surface wave tomography, and
mostly at large scale (Fichtner et al. 2017). For approaches aimed at
exploiting body waves, the consideration of the sources is currently
limited to the implementation of a strategy to simulate the correla-
tions themselves, based on surface sources (typically ocean forcing,
Ermert et al. 2016; Tomasetto et al. 2024, 2025) or the passing of
freight trains on the smallest scale (e.g. Sager et al. 2021), with the
possibility in this case of also quantifying the sensitivity kernels to
the structure.

Moving away from the interpretation of correlations, assuming the extracted mantle reflected waves are part of the local Green's function, requires a detailed understanding of the underlying interferences at play. In particular, a significant pulse observed in the correlations should not be a priori interpreted as its Green's function counterpart (assuming it exists). This situation is particularly true in the case of coda correlations of large long-period earthquakes (e.g. Kennett & Pham 2018; Pham et al. 2018). Even though we are not specifically dealing with late coda propagation, the underlying principle remains the same. Following the terminology introduced by Wang & Tkalčić (2020), we refer to the prominent signals observed in correlation functions as features (denoted by a*), a term used to distinguish them from physical seismic phases. Each feature results from the constructive interference of multiple real seismic phases, whose individual contributions are referred to as constituents (i.e. correlation cross-terms). A key challenge in seismic interferometry is the presence of spurious signals, or spurious features, that do not correspond to any actual phase in the Green's function (e.g. Snieder et al. 2006; Li et al. 2020a). Understanding the nature of these spurious signals requires a careful examination of how different constituents combine within a given feature (e.g. Boué & Tomasetto 2024). By analysing this process in controlled cases, both in simple scenarios and in more realistic distributions, we aim to develop a framework for interpreting the visible correlation features. This approach is particularly relevant for studies such as those by Pedersen et al. (2022) and Lu et al. (2023), where strong signals appear in seismic correlations.

In this article, we propose to describe Pv410p* and Pv660p*
features and its constituents in the case of distant ocean sources and
explore the time and amplitude bias of the noise correlations with
respect to the local Green’s function. We first study the properties
of the Pv410p and Pv660p waves and their constituents in Pv410p*

and Pv660p* from a ray approach. We then model Pv410p* and
Pv660p* using both a vertical point force and extended source.
Finally, we model correlations and compare the observed times and
410/660 amplitude ratios using a realistic source models with the
objective to use modelling of ambient noise cross-correlations to
interpret observations in the greater Alpine area.

## 2 CONSTITUENTS OF P-WAVE REFLECTIONS ON THE MTZ IN THE CORRELATION—PV410P* & PV660P*

This section sets the scene for understanding the results of the full
waveform modelling of Section 3.1. We first present the definitions
and methods we use, the target Green’s function and ray-based
analysis of the possible constituents of the *P*-wave reflections on
the MTZ in the correlation functions in the case of distant sources.

2.1 Definitions and methods

Throughout this work we use the wave phase conventions of
Crotwell *et al.* (1999). The *P* waves reflected on the top and bottom
mantle transition discontinuities are therefore named Pv410p and
Pv660p, while an incoming core phase followed by rebound and
reflection would be named, for example, PKIKPPv660p. In some
cases and to simplify, we’ll refer to either or both PKPab and PKPbc
phases as PKP, distinguishing them only where necessary. Note that
the mantle discontinuities are at 410 and 660 km depth in the global
AK135F seismological model (Kennett *et al.* 1995), which is why
we refer to the mantle discontinuity reflected *P* waves as Pv410p
and Pv660p, even though the discontinuity depth is model depen-
dent. We will use AK135F throughout this work, but we tested our
procedure and conclusions on other models (e.g. Kennett 1991).
We use the term b-caustic for the particular focusing effect, due to
the Earth’s velocity structure, that occurs at the intersection of the
slower (PKPab) and faster (PKPbc) core seismic branches, where
the derivative of the ray parameter with respect to the epicentral dis-
tance is infinite. This well-known condition in global seismology
corresponds to a point at approximately 145° epicentral distance.

We will be studying (see Section 3.1) the bias embedded in
the *P* waves reflected on the mantle discontinuities which are ex-
tracted through noise correlations. Because we work with subver-
tical waves, the relevant Green's function is *G*<sub>zz</sub>, that is, a vertical
point force and the vertical component. Therefore, the correlations
will be calculated between vertical components. We will henceforth
drop the *zz* index on the cross-correlations, which we will name *CC*.

When we extract the reflected waves through cross-correlations
we do not obtain the exact Green’s function, nor the exact Pv410p
and Pv660p phases. The pseudo Pv410p and Pv660p obtained were
labelled Pv410p* and Pv660p* by Wang & Tkalčić (2020) for coda
extracted waves, and we will follow this notation.

Throughout this article, we will use *x* for the interstation distance and Δ for the distance between the vertical point source and the nearest station. Note that the furthest station will be located at a distance of Δ + *x* from the source if the two stations are on the same great circle as the source. We will, for extended symmetrical sources, use Δ as the distance to the centre of the extended source patch to the nearest station.

We use the TauP functionalities of the Obspy python package
(Beyreuther et al. 2010) for the ray tracing at local and global scale.
The TauP functionalities in Obspy started out as part of the Java
TauP Toolkit by Crotwell et al. (1999).
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Cartesian ray representation of Pv410p (blue) and Pv660p (green) phases for vertical source (black arrow) to station (blue/green), the two being separated by the distance of $\Delta$. (b) Time-distance section on the vertical component of the recording of vertical source propagation. Blue/green lines correspond to the arrival time of Pv410p and Pv660p waves, respectively.</figcaption></figure>

We use AxiSEM (Nissen-Meyer *et al.* 2014) for all full wave-form simulations. Coherent with models of P-wave generation on the sea floor, and with the fact that we are studying subvertically propagating P waves, we simulate vertical displacements created by a vertical force. For extended sources or realistic source distributions, we use closely located vertical forces (see Section 3.2). The target frequency range is that of the second microseismic peak, which is the main frequency band for generating globally propagating P waves at global scale (Roux *et al.* 2005). We keep the default sampling rate of 4 Hz used in WMSAN for this study. We therefore carry out the simulations using a Gaussian source function centred on 7s period, that is, on the dominant period of the secondary microseismic peak period. We will for simplicity use the term $G_{zz}$ for these filtered Green's functions.

At a local scale we will use the simulations, that is, $G_{zz}$, directly. For the global scale we will use synthetic cross-correlations using WMSAN (Tomasetto *et al.* 2025). Using this framework, we compute synthetic cross-correlation functions from modelled teleseismic P-wave propagation and modulated by any source distribution at the surface. WMSAN offers two possibilities for defining the surface source geometry: users may either supply a source distribution of their choosing or take advantage of a state-of-the-art secondary microseismic source model derived from WAVEWATCH III hinders (Tolman 2009), which provide azimuthal pressure spectra of the sea surface. These source models include frequency-dependent vertical force at the seabed, derived from oceanic wave-wave interactions in the 0.10–0.34 Hz band (Ardhuin *et al.* 2011), and local bathymetric amplification, which is specific depending on the type of seismic wave considered (Gualtieri *et al.* 2014). For seismic propagation at teleseismic distances, synthetic Green's functions are computed using AxiSEM (Nissen-Meyer *et al.* 2014) for a vertical point force within a global radial Earth model, AK135F (Kennett *et al.* 1995), and filtered in the relevant frequency band. Possible deflection and complex scattering of teleseismic waves due to the heterogeneous 3-D mantle, presumably weak, is not taken into account in this study, but could be explored by complexifying this Green's function data base. Synthetic cross-correlation functions are then obtained using the representation theorem by combining Green's

functions with the source spectral density at each ocean grid point to simulate the wavefield recorded between seismic station pairs, following the formulation of Sager *et al.* (2021). As radial Earth models do not account for lateral heterogeneities, scattering does not act to reduce the amplitude of surface waves, leading to an unrealistic representation of these waves compared to body waves in this frequency band (Li *et al.* 2020b; Tomasetto *et al.* 2024). In other words, strong scattering of surface waves in the frequency range favours the emergence of weaker body waves far from the source. To focus on the seismic phases of interest and remove unrealistically coherent surface wave cross-terms in the synthetic correlations (and thus mimicking scattering attenuation), WMSAN provides a simple way to window specific seismic phases in the Green's function data base. We here select a time window ranging from PKIKP to PKPabPv660p arrival times. As we here focus on time and amplitude biases in the extracted cross-correlations as compared to $G_{zz}$ and, as the synthetic seismograms are, practically, noiseless, we will not apply techniques aimed at improving the signal-to-noise ratio, for which an ample literature is available. In practice we do not need to stack the correlations, nor apply trace normalization prior to correlation (e.g. whitening, one-bit correlations).

## 2.2 Target $G_{zz}$ at short epicentral distances

Fig. 1 shows the geometry of Pv410p and Pv660p rays (Fig. 1a) and $G_{zz}$ waveforms as a seismic section, focusing on Pv410p and Pv660p (Fig. 1b). The seismic section focuses on short distances, that is, 0° to 1.5°. The seismograms are dominated by surface waves that arrive at times up to 100 s within this distance window. We also observe Pv410p and Pv660p arrivals at 100 and 150 s. Note that Pv410p is masked by surface waves at distances over 1.5°, coherently with the observation of Pedersen *et al.* (2022) on noise correlations in the greater Alpine region.

The amplitude difference between Pv410p and Pv660p stems from the difference between the reflection coefficients at the two discontinuities, the transmission coefficient across the 410 discontinuity and the geometric attenuation, due to the longer path of
<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) Ray representation of long distance *P*-wave phases (PKPab and PKPbc, in orange) and their extent to reflections at 660 km depth (PKPabPv660p and PKPbcPv660p, in green) for a source (star) located at the surface and a pair of receivers (triangles). The orange station is located where the orange and the green rays arrive at the surface, that is, after the reflection on the 660 discontinuity. (b) Zoom on PKPab (orange) and PKPabPv660p (green) near both receivers. (c) Zoom on PKPab (orange) and PKPabPv410p (blue). (d) Zoom on PKIKP (orange) and PKIKPPv660p (green), in the case where the first reflection at the surface of the green ray is located at the orange station.</figcaption></figure>

the phase Pv660p as compared to the path of Pv410p. This amplitude ratio will vary with distance, dominantly due to variations in reflection coefficient with incidence angle.

Previous studies have used as an implicit hypothesis that the Pv410p* and Pv660p* arrivals extracted from noise correlations correspond to a frequency filtered version of $G_{zz}$. We demonstrate that this is not the case. Section 3.1 contains more realistic scenarios for which we will compare differences in traveltimes and relative amplitudes with those observed in the simulation of Fig. 1, using distance sources.

## 2.3 Cross-correlation terms for global propagation: ray theory approach

Similar waveforms of different waves will correlate well, and the cross-correlation will have a maximum peak (or maximum envelope, in case of a phase shift) at the traveltime difference between the two waves. These are the so-called *cross* — *terms*. In this section we explore through ray theory the geometrical aspects (slowness and traveltime differences) of the cross-terms that can contribute to the features Pv410p* and Pv660p*, that is, we explore the major constituents of the two features.

Our target study region is Europe for which the study of Pedersen et al. (2022) showed that sources at very low latitudes are the ones that contribute efficiently to Pv410p* and Pv660p*. As Pv410p and Pv660p are subvertical, we need to recover cross-terms with a very small slowness (0-5 s deg⁻¹). Such waves are related to the PKP and PKIKP branches, that is, beyond the Earth’s outer core shadow zone. Fig. 2 shows some examples of constituents which would contribute to Pv410p* and Pv660p*. For simplicity we cite

below only the 410 reflected waves, but ‘410’ can in all cases be replaced by ‘660’.

As shown in Fig. 2, the constituents of Pv410p* will be waves that have had the extra reflection at the 410 discontinuity. This would, for example, be the case of the correlation of $CC_{PKIKPPv410p,PKIKP}$, that is, a PKIKP wave at station 1 (nearest the source) with a PKIKP410vp wave at station 2 (equivalent to Fig. 2b for the case Pv660p*). Depending on the distance *x* between the two stations and the distance *Δ* between the source and station 1, contributing cross-terms, or constituents, can also come from the two PKP branches PKPab and PKPbc, namely $CC_{PKPabPv410p,PKPab}$, $CC_{PKPbcPv410p,PKPbc}$, or even cross-branch contributions $CC_{PKPbcPv410p,PKPab}$ and $CC_{PKPabPv410p,PKPbc}$. The cross-correlation will contain the sum of all of these cross-terms. As both Pv410p and Pv660p waves arrive subvertically, with small slownesses (<5 s deg⁻¹ in our study), we expect no other cross-terms, such as multiple *P*-wave reflections that have higher slowness (>5 s deg⁻¹). Note that the case of source side cross-terms are discussed in Section 5.

For a given interstation distance *x*, the ideal situation is (in the case of the PKIKP phase) when the source is located at a distance *Δ* where PKIKP and PKIKPPv660p have exactly the same slowness. In that case, the reflection point on the surface of the PKIKPPv660p phase matches exactly the arrival point of the PKIKP phase (Fig. 2d). However, even when this is not the case, the cross-correlation will still yield a relevant cross-term but it does not have the same reflection point, nor the same traveltime as those of Pv410p between stations 1 and 2 (see Fig. 2b). Those are the root causes for the bias that we explore in Sections 3 and 4.

The concepts of Fig. 2 can be further explored in terms of slowness and arrival times. Fig. 3 shows the arrival times and slowness
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Slowness (upper panel) and time (lower panel) as a function of source–receiver distance (Δ) and/or interstation distance (x) for multiple phases in AK135F. Each branch with a specific colour correspond to a specific phase (e.g. PKPab, Pv410p,...). In each panel horizontal grey lines correspond to Pv410p (dashed) and Pv660p (dotted) slownesses at x = 1.5°. Grey vertical lines correspond to the intersection of the previous horizontal lines and the branches of the other phases, giving the corresponding distances for each crossed branch. In the bottom panel, horizontal grey lines are the projection on the arrival time axis of the previously obtained distances for the given branches. The difference between the arrival times of two grey lines with the same line style is read as dt and corresponds exactly to the arrival time of the simple reflection (e.g. Pv410p) between stations spaced x° apart.</figcaption></figure>

as a function of distance. It is organized such that we have the small
interstation distances *x* on the left and the source to station dis-
tance Δ to the right. Fig. 3 shows (grey lines) the example of two
stations with an interstation distance of *x* = 1.5°. In that case the
slowness of Pv410p is 2.5 s deg⁻¹. The ideal source distance and
branch is then PKPbc, with station 1 at 149° distance and station 2 at
149°+1.5° distance. The difference in arrival time of PKPbcPv410p
(149°+1.5°) and PKPbc (149°) will yield exactly the target arrival
time of Pv410p. However other phases will contribute as well and
so will the target phase from other non-optimal distances in the case
of extended sources, such as storms.

Before moving to full waveform modelling, we first show
the cross-correlations on synthetic seismograms based on the
arrival times in Fig. 3. We constructed synthetic seismograms
by placing a Gaussian function centred on each of the arrival
times of the target phases discussed previously. We then corre-
lated the seismograms for two stations for different interstation
distance x. These cross-correlations will have local maxima at

times corresponding to traveltime differences between different
waves.

The outcome of this simplified approach is shown in Fig. 4 for
*x* = 0° (Fig. 4a), 1°(Fig. 4b), 2° (Fig. 4c) and 3° (Fig. 4d). We carried
out this for different source-station distances (Δ). By construction
we observe the different cross-terms appear, with each cross-term
having some bias as compared to the target Pv410p arrival time (blue
dashed lines) for the given *x*. Note that the cross-terms between the
PKPPv410p and PKP phases from the ab and bc branches (e.g.
*CC*<sub>PKPabPv410p,PKPbc</sub>) are ahead of the Green's function arrival
times (here *t*<sub>Pv410p</sub>).

For real field data we will be stacking correlations using mul-
tiple station pairs. The extracted features Pv410p* and Pv660p*
will therefore be the sum of correlations with different interstation
distances *x* and different source-station distances Δ. The cross-
branch cross-terms (between PKPab and PKPbc branches) vary
rapidly with Δ, so they can be expected to decrease in amplitude
for extended sources. On the contrary, the cross-terms within the
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Synthetic correlation function as a function of epicentral distance for a point source, using as reference Green's function for the calculation of the cross-correlation in WMSAN (Tomasetto et al. 2025), Gaussian pulses located at PKPab, PKPbc, PKPabPv410p and PKPbcPv410p arrival times for an interstation distance of (a) $x = 0°$, (b) $x = 1°$, (c) $x = 2°$ and (d) $x = 3°$. Dark blue dashed lines represents Pv410p($x$) arrival time in the Green's function of the P-wave reflection on the 410 discontinuity for $x$ spaced stations. Light blue lines is the theoretical traveltime difference (ray theory) between PKPab and PKPabPv410p branches, and between PKPbc and PKPbcPv410p. In orange is shown the theoretical traveltime between PKPab and PKPbcPv410p, and PKPbc and PKPabPv410p cross-terms. The formulation used in the figure is chosen to show that we observe all the cross-terms between two recorded wavefields, one at the first station (at distance $\Delta$) and one at the second station (at distance $\Delta + x$).</figcaption></figure>

same branch (e.g. PKPabPv410P and PKPab) are likely to be robust also for extended sources. These terms are almost always in advance with respect to the arrival time of the target Pv410p and Pv660p. The major complication of this simple view is that the effect of extended sources may be significantly influenced by the amplitudes of the waves. For example, we can assume that the high amplitudes around the triplication source-station distance may play a particularly strong role in the cross-terms. This is why bias cannot be estimated from ray geometry alone; we need to use full waveform modelling and explore the effect of extended sources to be able to take into account the relative amplitude of the different seismic waves.

### 3 BIAS BETWEEN $G_{zz}$ AND NOISE CORRELATIONS USING SIMPLE MODELS OF DISTANT SOURCES

In this section, we explore differences in arrival times and amplitudes between the features Pv410p* and Pv660p* in the cross-correlation and the Pv410p and Pv660p in the Green's function $G_{zz}$. We will call these differences *bias*, as they correspond to a bias in interpretation when Pv410p* and Pv660p* as extracted from seismic noise are considered to be Pv410p and Pv660p.

For each interstation distance $x$ we calculate the theoretical arrival times of Pv410p and Pv660p, which we will call $t_{Pv410p}(x)$
and $t_{Pv660p}(x)$. As a reminder, the output of the modelling are cor-
relations, for a given source position and a set of stations. For
this source position and each station pair we calculate the arrival
time $t_{Pv410p*}(\Delta, x)$ of Pv410p* by the time of the maximum en-
velope within a window of 25 s around $t_{Pv410p}(x)$. We next calcu-
late the bias as $t_{Pv410p*}(\Delta, x) - t_{Pv410p}(x)$. We do the same thing
for Pv660p*, then calculate the bias on differential traveltimes as
$(t_{Pv660p*}(\Delta, x) - t_{Pv410p*}(\Delta, x)) - (t_{Pv660p}(x) - t_{Pv410p}(x))$.

A positive bias on Pv410p* or Pv660p* means that the maximum
envelope arrives later than the theoretical arrival time, and a positive
differential traveltime bias means that the differential traveltime
between the two reflections is larger than for the Green's function.

Note that we use the maximum envelope rather than the max-
imum of the correlation itself. Even for the simulations, we have
a complex phase of the cross-terms, due to the fact that we corre-
late waves with and without reflection at the free surface. On the
contrary, multiscattered waves in the Earth would contribute to the
Green's function. The sum of the two would potentially result in
a phase shift. Therefore, by using the envelope we avoid these ef-
fects and in practice mimic observations on field data where the
time of the maximum envelope is more stable than the time of
peaks and troughs. Using the envelope has the additional advantage
of avoiding potential $\pi/2$ phase offset issues related to stationary
point integration when comparing correlation to its Green's function
counterpart.

We also explore amplitude biasing. Already, we know that the
difference in the geometrical spreading is almost negligible be-
tween, for example, PKPab and PKPabPv410p, while it is approxi-
mately 660/410 in $G_{zz}$. We therefore calculate the amplitude ratios
$Ar_{CC}(\Delta, x) = A_{Pv410p*}(\Delta, x)/A_{Pv660p*}(\Delta, x)$. We do the same for
$G_{zz}$ by calculating $Ar_{GF}(x) = A_{Pv410p}(x)/A_{Pv660p}(x)$ on the syn-
thetic seismograms output from the Axisim simulation. Finally, we
correct for the difference in geometrical spreading. Once corrected
for geometrical spreading, a ratio $Ar_{CC}/Ar_{GFcorrected}$ of one corre-
sponds to no bias.

According to Fig. 4, we expect a dominance of negative time
bias. We also expect to observe an important difference in terms
of amplitudes in the cross-correlation for pairs located near to the
b-caustic and pairs that are far away from it. Finally, the relative
dominance of the PKPab, PKPbc and PKIKP phases will depend
on $\Delta$, and their slowness determines how close $x$ is to an optimal
situation (see Fig. 3). We therefore calculate biases for distances $x$
of less than 3° and for source-station distances $\Delta$ between 140° and
175°.

We first use a distant vertical force as source to explore what
emerges in the CC. We next study the time differences using an
extended, therefore more realistic, source to represent oceanic mi-
croseismic events.

**3.1 Bias for distant point source**

Within this subsection we give a first insight into the effect $\Delta$ and $x$ on bias, using a vertical point source at teleseismic distance. As we are using full-waveform modelling, the effect of relative amplitude between different phases will be respected, in particular the effect of the b-caustic.

Fig. 5 shows the correlation amplitude for the 410 cross-term,
and the bias as defined above. The maximum amplitude (Fig. 5a)
of the cross-correlation for the feature Pv410p* is reached for $\Delta$
between 145° and 150°, corresponding stations located at distances
around the b-caustic. For the direct PKP, in AK135, the b-caustic is

at 145° from the source, while the b-caustic with rebound is at 147°
(for PKPPv410p) and 148.7° (for PKPPv660p) from the source.
This means that we only have an effective cross-correlation at high
amplitude once the furthest station with the rebound (in the case of
PKPPv410p) is at distance $\Delta + x > 145°$, that is, $\Delta > 145° - x$.
Therefore, the smaller the value of $x$, the higher the value $\Delta$ for
which we observe high amplitude of cross-correlation related to the
cross-term in the correlation. Note that the results for Pv660p* (not
shown) are qualitatively the same, shifted towards higher values
of $\Delta$.

The time biases on Pv410p and Pv660p are, as expected, mainly
negative. They vary rapidly with both $\Delta$ and $x$, in a way not easily
predictable. Overall, for $\Delta < 160°$, the bias is mainly between be-
tween -2 and 0 s, for intermediate values of $x$. The high biases at
for large $x$ are not of so high relevance, as the interstation distance
is generally limited to 1.5° or 2° for field data. However, for small
values of $x$ we observe that biases can be high, with a strong $\Delta$
dependence. These high biases for small $x$ is significantly more
prominent for Pv410p*. This can be explained by Fig. 3, where we
see that very small values of slowness, as needed for Pv660p*, are
related to the PKIKP phase only. The biases increase systemati-
cally for $\Delta > 160°$. This is the distance where we loose the PKPbc
branch (see Fig. 3), and what remains are the competing effects of
the PKPab branch and PKIKP. Therefore, we move from three con-
tributing phases to two, which may be why we observed particularly
high biases and variability with $\Delta$.

The differential traveltimes show that overall the biases have
no systematic shift, but rather varies around, for $\Delta < 160°$. This
is a first indication that traveltime differences will have less bias
than absolute traveltimes. We do, however, again observed highly
variable and strong biases for $\Delta > 160°$. The amplitude ratio
$Ar_{CC}/Ar_{GFcorrected}$ is also more variable for $\Delta > 160°$, where it
can take values above 5 or below 0.2. Overall, the amplitude biases
are in the range 0.4-2.5, and only some combinations of $x$ and $\Delta$
give values between 0.75 and 1.25.

The positive outcome is that, once correlations are stacked over
different distances to storms and across different interstation dis-
tances, we can hope that bias will be reduced or, at least, the influ-
ence of fast varying cross-terms will decrease. Fig. 5 does, however,
indicate that there may be some conditions for small biases: that in-
terstation distances are bigger than some minimum value and that
storm distributions include sources at less than 160° distance, the
latter because those sources are likely to dominate in amplitude over
the far sources, with the result of avoiding large bias.

**3.2 Extended source: approaching realistic source geometries**

We now use a spatially extended patch to mimic ocean sources,
therefore including the interaction of the contributions from differ-
ent cross-terms. To do this, we use a 5° wide Gaussian patch (i.e.
$\sigma = 2.5°$), composed of closely spaced vertical point sources on the
surface. Because the strongest biases (see Fig. 5) vary rapidly with
$\Delta$, we expect an extended source to reduce unwanted contributions.
As a reminder, $\Delta$ designates the distance between the centre of the
patch and the nearest station.

In Fig. 6 we show cross-correlations of a point source (upper
panel) and of the extended source (lower panel) for $x = 2°$ and
for 140° < $\Delta$ < 175° (lower panel). Recall that, $x = 2°$ means that
the stationary point for Pv410p* is at the caustic (see Fig. 3). The
cancelling of cross-terms for the extended source as compared to
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Time and amplitude bias of Pv410p* and Pv660p* for a point source. Panel (a): Normalized maximum amplitude of cross-correlation functions, in the time window of 90–115s, as a function of interstation distance, *x*, and epicentral distance, Δ. Panel (b): Time difference between the lag time associated with the maximum amplitude from upper panel and *t*<sub>Pv410p(x)</sub> as a function of *x* and Δ in colourmap of 1 s levels. Panel (c): Time difference between the lag time associated with the maximum amplitude from upper panel and *t*<sub>Pv660p(x)</sub>. Panel (d): Time difference of panel (b) and (c) calculated as (*t*<sub>Pv660p*(Δ, *x*) − *t*<sub>Pv410p*(Δ, *x*) − (*t*<sub>Pv660p</sub>(*x*) − *t*<sub>Pv410p</sub>(*x*))). Panel (e): Amplitude ratio bias *Ar*<sub>CC</sub>/*Ar*<sub>GFcorrected</sub>.</figcaption></figure>

the point source is spectacular. Clearly, the b-caustic will stabi-
lize measurements also when the Gaussian source patch is located
at significant distance from the caustic: even a small amplitude
source at the b-caustic will dominate over sources in any other
location.

Considering that our simulations cover a radially layered Earth,
we would potentially see strong cross-terms also from reflections
on the mantle discontinuities on the source side. It turns out that
these cross-terms become negligible because the multiple stations
on the receiver side play a similar role as the extended source
to delete unwanted cross-terms. For more detail, we refer to Ap-
pendix A. Another issue is that sources located at intermediate
distances could also produce ambiguous cross-terms. This issue is
handled by the processing of field data, by choosing days with low
H/V ratio. As discussed in Section 4, the average distribution of

noise sources over days with low H/V has insignificant contribution
by sources outside PKP main phases. Therefore, many cross-terms
(such as PKPbcPv410pPKPbc, PKPbcPKPbc, from sources at ap-
proximately 70° distance) are *de facto* eliminated from the noise
correlations.

We will now explore how the stabilization due to the extended
source influences the biases, following the same procedure as in
Section 3.1. The outcome is shown in Fig. 7.

We first observe that the high amplitude area spreads over a
wide range of distances Δ. We also observe that cross-terms,
with associated large biases, that appeared for a point source at
165° < Δ < 170° are less energetic (one order of magnitude). This
confirms what we observed for the correlations: the b-caustic will
dominate even if only very small amplitude parts of the Gaussian
patch overlaps with the b-caustic point.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Synthetic correlation function as a function of epicentral distance for a point source and a Gaussian patch source of 5° width, with an interstation distance *x* = 2°. Upper panel: point source cross-correlation waveforms as a function of time and distance between source and closest station (Δ). Lower panel: Extended (Gaussian) patch source cross-correlation waveforms as a function of time and distance between the centre of the source patch and closest station (Δ). Dark blue dashed lines represents Pv410p(*x*) arrival time in the Green's function of the *P* wave reflection on the 410 discontinuity for *x* spaced stations. Dark green dashed lines represents Pv660p(*x*) arrival time in the Green's function of the *P* wave reflection on the 660 discontinuity for *x* spaced stations. Light blue lines are the cross-terms between PKPab and PKPabPv410p branches, the cross-terms between PKPbc and PKPbcPv410p and the cross-terms between PKIKP and PKIKPPv410p. Light green lines are the cross-terms between PKPab and PKPabPv660p branches, the cross-terms between PKPbc and PKPbcPv660p and the cross-terms between PKIKP and PKIKPPv660p. Orange lines are the cross-branches cross-terms for each discontinuity: the cross-terms between PKPab and PKPbcPv410p, the cross-terms between PKPbc and PKPabPv410p, the cross-terms between PKPab and PKPbcPv660p and the cross-terms between PKPbc and PKPabPv660p.</figcaption></figure>

The time biases have largely stabilized and decreased, both for small and big values of *x*. The spurious, very delayed, waves (large positive biases) have vanished both for Pv410p* and for Pv660p*. The biases for both Pv410p* and for Pv660p* are now, generally, smaller than 2 s. They remain negative, coherent with the systematic effects expected from the ray analysis of Section 2. The biases remain small for Δ up to 170° for Pv410p*, while they are small for Δ up to 162° for Pv660p*. However, the biases have decreased, and biases are small (less than 2 s) extending to 172° also for Pv660p* if *x* is in the range 0.5°–1.2°. This is likely due to the subvertical PKIKP phase dominating rather than various spurious phases, which have decreased in amplitude.

The decrease and stabilization of time biases is particularly clear if we consider the differential traveltime between Pv660p* and Pv410p*. If we consider analysis for field data where *x* < 2°,

the time differences are within ±2 s, except in the combination where Δ > 162° and *x* > 1°.

The bias in amplitude ratio is spectacularly reduced, going from an almost random and strongly oscillating behaviour for the point source to a smooth pattern for the patch, the latter with biases less than 25 per cent for broad ranges of Δ and *x*. However, we observe that the amplitude ratio $Ar_{CC}/Ar_{GFcorrected}$ is very high for several combinations of Δ and *x*, meaning that the amplitude of Pv660p* is significantly underestimated as compared to Pv410p*. For Δ > 160°, bias is mainly linked to the Pv660p*. For the high bias patch centred at Δ = 145° we see the effect of the b-caustic with reflection arriving at higher values for Pv660p* than for Pv410p* (see Section 3.1). On the other hand, the correlations will capture major changes. If the amplitude ratio in field data changes more than a factor 2.5, it is likely to be a local effect rather than a source effect.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Time and amplitude bias of Pv410p* and Pv660p* for a spatially extended Gaussian patch source. Panel (a): Normalized maximum amplitude of cross-correlation functions, in the time window of 90–115s, as a function of interstation distance, *x*, and epicentral distance, Δ. Panel (b): Time difference between the lag time associated with the maximum amplitude from upper panel and *t*<sub>Pv410p(x)</sub> as a function of *x* and Δ in colourmap of 1 s levels. Panel (c): Time difference between the lag time associated with the maximum amplitude from upper panel and *t*<sub>Pv660p(x)</sub>. Panel (d): Time difference of panel (b) and (c) calculated as (*t*<sub>Pv660p*</sub>(Δ, *x*) − *t*<sub>Pv410p*</sub>(Δ, *x*) − (*t*<sub>Pv660p</sub>(*x*) − *t*<sub>Pv410p</sub>(*x*)). Panel (e): Amplitude ratio bias *Ar*<sub>CC</sub>/*Ar*<sub>GFcorrected</sub>.</figcaption></figure>

4 ON THE USE OF OCEANIC MODELS
TO CHARACTERIZE AMBIENT NOISE
CORRELATION OBSERVATIONS

The aim of this section is to understand and quantify potential bias in
field data. We therefore model bias for our future target, which is the
MTZ beneath the greater Alpine region. For this type of simulation
we need to mimic realistic oceanic sources and data processing.

We will assume that we optimize the stacks using days with
low H/V ratio (Pedersen et al. 2022). As Switzerland is placed
centrally as compared to the greater Alpine region, we use a source
distribution based on days where Lu et al. (2021) observed low H/V
ratio in Switzerland. We selected days with H/V < 1, using H/V
from (Lu et al. 2021). On these days, the subvertically propagating
P waves are expected to be significant over the surface waves (Lu
et al. 2021).

For each day of these low H/V days (640 d between 2011 and
2019), we use the oceanic forcing model WAVEWATCH III (Tol-
man 2009), which models the strength of the seismic P-wave source
in the second microseismic peak, taking into account wave activ-
ity and their interaction with the sea bottom. This model has been
shown to provide a good proxy to sources as observed in seismic
noise (e.g. Ardhuin et al. 2011; Stutzmann et al. 2012; Gualtieri
et al. 2014). Each model is expressed as an amplitude map (see
an example of daily source amplitude in Fig. 8b). We use WM-
SAN (Tomasetto et al. 2025) to simulate daily cross-correlations
between pairs of stations centred at latitude 46°48' and longitude
8°13' (central part of Switzerland). For each day we use as source
distribution the source amplitudes from WAVEWATCH III, with
a Gaussian spectrum centred on the secondary microseismic peak
(~ 6 s).
<figure><img src="image_9.png" /><figcaption>**Figure 8.** (a) Profile of seismic stations in Switzerland (red triangles). The blue arrow corresponds to the interstation distance *x*. (b) Map of ocean noise sources for a given day with an H/V ratio <1. (c) Normalized stack of normalized ocean noise source maps used in our simulation. The red triangle shows the location of the station profile (a) and the orange curves are the points 145° away from this profile, corresponding to the b-caustic distance. Source maps are outputs of WMSAN (Tomasetto *et al.* 2025).</figcaption></figure>

In order to improve the robustness of the measurements, it is common usage in noise correlations processing workflow to stack, for each station pair, the cross-correlation over all selected days to obtain a mean cross-correlation (Bensen *et al.* 2007). In our case, the signal-to-noise ratio of Pv410p* and Pv660p* is enhanced by normalizing each mean cross-correlations, then stacking the cross-correlations within a given study-area and for which station distances is below a given threshold, here we take 1.5° (Pedersen *et al.* 2022). We mimic the spatial stacking by using station pairs with different orientations as compared to the central point, oriented at 0°, 45°,... 335° (see Fig. 8a).

It is not straightforward to illustrate the source distribution weighted by the correlation normalization factor. Above that, we have the problem of the b-caustic. Note that sources close to the b-caustic will create higher cross-correlation amplitudes than sources located elsewhere, and the daily correlations are normalized prior to stacking them. Fig. 8 shows a proxy for the weighted source distribution. All these daily source distribution maps are then normalized and stacked to obtain the map in Fig. 8(c). Note that the normalized source distribution has significant amplitudes at distances Δ smaller and bigger than the b-caustic (orange curve in Fig. 8c). Also, there are significant source amplitudes almost everywhere within the circle defined by Δ < 145°.

The result of this simulation is presented in Fig. 9, where we show bias as defined in Section 3. To mimic the spatial stack done

on field data, hence to maximize the SNR at a given reflection point, we also stack the cross-correlations over all azimuths (see above) and calculate the bias for this stack (solid black line).

We observe that the time bias for Pv410p* is within the range -2.5 to 0.5 s for individual cross-correlations (light grey lines, covering different azimuths). We observe that if we stack the cross-correlations over azimuth and, then, calculate bias (bold black line), we observe that the bias is within the range -1.5 to 0.5 s, increasing with *x*. Finally, if we also stack the correlations over distance (applying Normal Move Out correction before stack) before calculating the bias (which is independent of *x* but which we illustrate with the orange line), we observe a bias of -1.1 s. The time bias for Pv660p* behaves qualitatively the same way, but the bias for the stacked correlations is almost constant being -1.25 s if the correlations are stacked over distance. This value does not change if correlations are also stacked over distance. The bias of the differential traveltime varies between approximately -2 and 2 s, while the bias of the azimuth stack decreases with *x* from +0.5 s to -2 s. The bias on differential traveltimes is -0.15 s once stacked over distance as well.

Overall this means that the absolute depth would be slightly underestimated both for the 410 and 660 discontinuities for the imposed, but realistic, source distribution. In real data, we prefer to use the relative depth rather than absolute depths to avoid any effect of shallower crust or upper-mantle properties. A broad distribution
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Time and amplitude observations on ambient noise cross-correlation simulations with oceanic wave forcing. Each grey curve on all panels correspond to one pair of station result (always including the middle one) from Fig. 8(a). Black curves correspond to each observation made on the stack of all cross-correlations. Orange curves correspond to a Normal Move Out corrected observation of the stack. (a) Time delay between $t_{Pv410p*}$ and $t_{Pv410p_{GF}}$ in seconds. (b) Time delay between $t_{Pv660p*}$ and $t_{Pv660p_{GF}}$ in seconds. (c) Time delay between $t_{Pv660p*} - t_{Pv410p*}$ and $t_{Pv660p_{GF}} - t_{Pv410p_{GF}}$ in seconds. (d) Amplitude ratio $Ar_{CC}/Ar_{GF}$ corrected.</figcaption></figure>

of distance and azimuth ranges of interstation profiles significantly stabilizes the measurements, and the differential traveltime bias becomes negligible.

The effect of the stack is even more pronounced on the amplitude bias. Note that no bias corresponds to the value 1. While individual directions may have significant bias (values of 0.5 for some azimuths), the stack over directions reduce the bias to be within the range 0.9 to 1.2, that is, within -10 to +20 per cent of the correct amplitude ratio between the two reflections. The stack over distance gives the value 0.9, that is, we overestimate the 660 reflection amplitude as compared to the 410 amplitude by approximately 10 per cent as compared to the Green's function. This means that, for the target area, major variations (>> 10 per cent) in amplitude ratios of reflections from the two discontinuities are likely to be caused by temperature and/or compositional effects, rather than being an artefact due to source distribution.

We can conclude that by using stacks from several backazimuths, and several days of low H/V ratio, we give a good estimate of both delays and relative amplitudes on the MTZ in this realistic model, assuming the greater Alpine area being target study zone.

Strong caution must be taken as to whether our conclusions are generally applicable everywhere on Earth. We have set our tests on realistic sources within a precise geographic framework, which is that we wish to study the mantle transition zone in the greater Alpine area. We have seen that a source distribution with 'sufficient' sources at large distances from the study area is key for good imaging. Also, we have seen that the presence of sources at distances around the b-caustic play a major role in reducing bias on

the reflections on mantle discontinuities. The greater Alpine area is clearly favourable for such imaging while other locations might not be.

As we rely on the use of teleseismic P waves for deep Earth imaging, we can associate this P-wave production with microseismic events. An event is defined as an oceanic seismic forcing large enough to produce teleseismic P waves (Zhang *et al.* 2023). Zhang *et al.* (2022) has created a catalogue of these microseismic events between 1994 and 2020. Other P-wave event catalogues exist, based for example on the Japan Hi-Net network (Nishida & Takagi 2022). To illustrate the location dependent illumination from microseismic events, Fig. 10 shows, at each point on the Earth, the number of events that are located close to the b-caustic distance (145° ± 5°) from the point, using the microseismic events in 2020 from Zhang *et al.* (2022). We see that the Northern Hemisphere is well covered for MTZ imaging purposes using PKP phases generated by these secondary microseisms events. The north pole is particularly well illuminated. We also observe that the strongly localized north Atlantic and north Pacific events result in two dominant circular patterns at 145° of distances from these events in the Southern Hemisphere. East Antarctica is also well covered. However, large part of the Earth are also poorly covered, such as for example the north-western part of South America and large parts of the Indian and Pacific Oceans. This map should not be used to directly identify the feasibility of mantle discontinuity imaging with seismic noise but rather to encourage future studies to integrate modelling efforts and data-driven analysis of the noise field, for example H/V ratio and beamforming analysis and their variation with time.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Map of number of events located at b-caustic distance from any point on Earth. Each point correspond to the number of secondary microseisms events in 2020 from Zhang *et al.* (2022) catalogue, exceeding a force $2 \times 10^{12}$N (minimum force of meaningful events), located at $145^\circ \pm 5^\circ$ (b-caustic distance) from the chosen point. Colourbar corresponds to the number of events, from purple (0 event) to yellow (~2500 events).</figcaption></figure>

We conducted this synthetic experiment for a realistic source distribution and a pseudo-realistic station geometry located in our targeted greater Alpine region. For this configuration, we find that the combination of a favourable distribution of extended sources around the antipode, as well as the effect of azimuthal averaging and interstation distance, tends to constrain biases within reasonable limits and encourages comparison with real data. In particular, the azimuthal averaging has a significant effect in rendering the resultant effective source area much smoother and likely analogous to an extensive region as discussed in Section 3.2, despite source complexity, the latter being due to irregular bathymetry. It is likely, though not proven by this work, that the interior of the orange circle in Fig. 8(c) is homogenized through rotational symmetry, hence more effectively filling the Fresnel zone.

## 5 CONCLUSION

In this paper we explore major constituents properties of Pv410p* and Pv660p* features. By reconstructing the cross-correlation step by step from point source to oceanic forcing models, we can constrain better the potential bias between the observation and the Green's function. We show that the b-caustic plays a major role in cross-correlation amplitudes and bias results. Time bias are relatively small in any case of spatially extended (hence realistic) source. As we correlate signals emitted from localized noise sources, we show that Green's function amplitudes are not retrieved, but that this can be corrected for by taking into account the difference in geometrical spreading. In addition, we used a microseismic event catalogue to gain insight on the interpretation of seismic noise data for imaging the mantle discontinuities in the specific case of the greater Alpine area. Care should be taken when applying these results to other locations, but our workflow can be easily used elsewhere. It can potentially also be used for other body waves that probe the deep Earth. Further work would include taking into account how 3-D structure and attenuation impact the reconstitution of the target phases.

## ACKNOWLEDGMENTS

The authors would like to thank the Editor Christoph Senschönfelder, the assistant editor Fern Storey, Thanh-Son Pham and one anonymous reviewer for their time and excellent suggestions that helped to improve the manuscript. ISTerre is part of Labex OSUG@2020 (ANR10 LABX56).

## DATA AVAILABILITY

The ocean models used in this study come from WAVEWATCH III (Tolman 2009) and the secondary microseisms events catalog is from Zhang *et al.* (2022). To calculate synthetic noise correlations, we used the WMSAN python package from Tomasetto *et al.* (2025).

## REFERENCES

* Anderson, D.L., 1987. Thermally induced phase changes, lateral heterogeneity of the mantle, continental roots, and deep slab anomalies, *J. geophys. Res.: Solid Earth*, **92**(B13), 13 968–13 980.

* Ardhuin, F., Stutzmann, E., Schimmel, M. & Mangeney, A., 2011. Ocean wave sources of seismic noise, *J. geophys. Res.*, **116**(C9), C09004.

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

* Beyreuther, M., Barsch, R., Krischer, L., Megies, T., Behr, Y. & Wassermann, J., 2010. ObsPy: A Python Toolbox for Seismology, *Seismol. Res. Lett.*, **81**(3), 530–533.

* Boué, P. & Tomasetto, L., 2024. Opportune detections of global P-wave propagation from microseisms interferometry, *C. R. Géosci.*, **356**(S4), 85–100.

* Chevrot, S., Vinnik, L. & Montagner, J.-P., 1999. Global-scale analysis of the mantle Pds phases, *J. geophys. Res.: Solid Earth*, **104**(B9), 20 203–20 219.

* Christensen, U., 1995. Effects of phase transitions on mantle convection, *Annu. Rev. Earth Planet. Sci.*, **23**, 65–88.
Crotwell, H.P., Owens, T.J. & Ritsema, J., 1999. The taup toolkit: flexible seismic travel-time and ray-path utilities, *Seismol. Res. Lett.*, **70**(2), 154–160.

Ermert, L., Villaseñor, A. & Fichtner, A., 2016. Cross-correlation imaging of ambient noise sources, *Geophys. J. Int.*, **204**(1), 347–364.

Feng, J., Yao, H., Poli, P., Fang, L., Wu, Y. & Zhang, P., 2017. Depth variations of 410 km and 660 km discontinuities in eastern North China Craton revealed by ambient noise interferometry, *Geophys. Res. Lett.*, **44**(16), 8328–8335.

Feng, J., Yao, H., Wang, Y., Poli, P. & Mao, Z., 2021. Segregated oceanic crust trapped at the bottom mantle transition zone revealed from ambient noise interferometry, *Nat. Commun.*, **12**(1), 2531.

Fichtner, A., Stehly, L., Ermert, L. & Boehm, C., 2017. Generalized interferometry-I: theory for interstation correlations, *Geophys. J. Int.*, **208**(2), 603–638.

Goes, S., Yu, C., Ballmer, M.D., Yan, J. & Van Der Hilst, R.D., 2022. Compositional heterogeneity in the mantle transition zone, *Nat. Rev. Earth Environ.*, **3**(8), 533–550.

Gualtieri, L., Stutzmann, E., Farra, V., Capdeville, Y., Schimmel, M., Ardhuin, F. & Morelli, A., 2014. Modelling the ocean site effect on seismic noise body waves, *Geophys. J. Int.*, **197**(2), 1096–1106.

Kato, S. & Nishida, K., 2023. Extraction of mantle discontinuities from teleseismic body-wave microseisms, *Geophys. Res. Lett.*, **50**(18), e2023GL105017.

Kennet, B.L.N., 1991. IASPEI 1991 Seismological Tables, *Terra Nova*, **3**(2), 122–122.

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

Kennett, B.L.N., Engdahl, E.R. & Buland, R., 1995. Constraints on seismic velocities in the Earth from traveltimes, *Geophys. J. Int.*, **122**(1), 108–124.

Lawrence, J.F. & Shearer, P.M., 2006. A global study of transition zone thickness using receiver functions, *J geophys. Res.: Solid Earth*, **111**(B6).

Li, L., Boué, P. & Campillo, M., 2020a. Observation and explanation of spurious seismic signals emerging in teleseismic noise correlations, *Solid Earth*, **11**(1), 173–184.

Li, L., Boué, P., Retailleau, L. & Campillo, M., 2020b. Spatiotemporal correlation analysis of noise-derived seismic body waves with ocean wave climate and microseism sources, *Geochem. Geophys. Geosyst.*, **21**(9), e2020GC009112. doi:

Lu, Y., Pedersen, H.A. & Stehly, L. & AlpArray Working Group, 2021. Mapping the seismic noise field in Europe: spatio-temporal variations in wavefield composition and noise source contributions, *Geophys. J. Int.*, **228**(1), 171–192.

Lu, Y., Schmid, S.M., Wang, Q.-Y. & Bokelmann, G., 2023. Mapping the mantle transition zone discontinuities across South-Central Europe using body waves from seismic noise correlations, *Earth planet. Sci. Lett.*, **624**, 118457.

Nishida, K. & Takagi, R., 2022. A global centroid single force catalog of P-wave microseisms, *J. geophys. Res.: Solid Earth*, **127**(4), e2021JB023484.

Nissen-Meyer, T., Van Driel, M., Stähler, S.C., Hosseini, K., Hempel, S., Auer, L., Colombi, A. & Fournier, A., 2014. AxiSEM: broadband 3-D seismic wavefields in axisymmetric media, *Solid Earth*, **5**(1), 425–445.

Pedersen, H. & Colombi, A., 2018. Body waves from a single source area observed in noise correlations at arrival times of reflections from the 410 discontinuity, *Geophys. J. Int.*, **214**(2), 1125–1135.

Pedersen, H.A., Mattern, F., Poli, P. & Stehly, L., 2022. Imaging with seismic noise: improving extraction of body wave phases from the deep Earth through selective stacking based on *H/V* ratios, *Geophys. J. Int.*, **232**(2), 1455–1467.

Pham, T.-S., Tkalčić, H., Sambridge, M. & Kennett, B.L.N., 2018. Earth's correlation wavefield: Late coda correlation, *Geophys. Res. Lett.*, **45**(7), 3035–3042.

Poli, P., Pedersen, H.A. & Campillo, M. & the POLENET/LAPNET Working Group, 2012. Emergence of body waves from cross-correlation of short period seismic noise: Body waves from noise correlation, *Geophys. J. Int.*, **188**(2), 549–558.

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

Ruigrok, E., Campman, X. & Wapenaar, K., 2011. Extraction of P-wave reflections from microseisms, *C. R. Géosci.*, **343**(8–9), 512–525.

Sager, K., Tsai, V.C., Sheng, Y., Brenguier, F., Boué, P., Mordret, A. & Igel, H., 2021. Modelling *P* waves in seismic noise correlations: advancing fault monitoring using train traffic sources, *Geophys. J. Int.*, **228**(3), 1556–1567.

Shearer, P.M., 1991. Constraints on upper mantle discontinuities from observations of long-period reflected and converted phases, *J. geophys. Res.: Solid Earth*, **96**(B11), 18 147–18 182.

Snieder, R., Wapenaar, K. & Larner, K., 2006. Spurious multiples in seismic interferometry of primaries, *Geophysics*, **71**(4), SI111–SI124.

Stutzmann, E., Ardhuin, F., Schimmel, M., Mangeney, A. & Patau, G., 2012. Modelling long-term seismic noise in various environments: modelling seismic noise, *Geophys. J. Int.*, **191**(2), 707–722.

Tauzin, B., Debayle, E. & Wittlinger, G., 2008. The mantle transition zone as seen by global Pds phases: no clear evidence for a thin transition zone beneath hotspots, *J geophys. Res.: Solid Earth*, **113**(B8).

Tolman, H.L., 2009. User manual and system documentation of WAVEWATCH III TM version 3.14, *Technical note*, MMAB contribution, **276**(220). Available at: https://polar.ncep.noaa.gov/mmab/papers/tn276/MMAB_276.pdf

Tomasetto, L., Boué, P., Stehly, L., Ardhuin, F. & Nataf, H., 2024. On the stability of mantle-sensitive *P*-wave interference during a secondary microseismic event, *Geophys. Res. Lett.*, **51**(8), e2023GL108018.

Tomasetto, L., Boué, P., Ardhuin, F., Stutzmann, E., Xu, Z., De Plaen, R. & Stehly, L., 2025. Wmsan python package: From oceanic forcing to synthetic cross-correlations of microseis mic noise, *Seismica*, **4**(1).

Wang, S. & Tkalčić, H., 2020. Seismic event coda-correlation's formation: implications for global seismology, *Geophys. J. Int.*, **222**(2), 1283–1294.

Wapenaar, K., 2004. Retrieving the elastodynamic Green's function of an arbitrary inhomogeneous medium by cross-correlation, *Phys. Rev. Lett.*, **93**(25), 254301.

Zhang, R., Boué, P., Campillo, M. & Ma, J., 2022. Catalog of *P*-wave secondary microseisms events.

Zhang, R., Boué, P., Campillo, M. & Ma, J., 2023. Quantifying *P*-wave secondary microseisms events: a comparison of observed and modelled backprojection, *Geophys. J. Int.*, **234**(2), 933–947.

## APPENDIX A:

In this appendix, we explore the potential effect of the PKP phase source bounce (i.e. Pv410pPKP) in the cross-correlation function using ray theory. In the case of perfect symmetry of the velocity model, the Pv410pPKP and PKPPv410p phases have exactly the same slowness and arrival times at a given distance $\Delta$. Consequently, model sensitivity is shared by both phases, that is, 50/50. In the case of a spatially extended source, we present in Fig. A1 a diagram plotting the rays of the two phases Pv410pPKP and PKPPv410p for $\Delta$ varying from 150° to 153° distance from the orange station, with $x = 2^{\circ}$. Pv410pPKP phases are shown in purple, and PKPPv410p phases in blue. We note that the reflection points on the 410 km discontinuity for phases Pv410pPKP move like the source point (i.e. if $\Delta$ varies by 1°, the reflection point will also move by about 1°). We can see that for the PKPPv410p phase, even if the source moves, the reflection point is virtually unaffected. As we saw in Section 3.2, we expect the extended source to act to erase unwanted cross-terms (i.e. $CC_{Pv410pPKP,PKP}$). In real data, we suggest exploring cross-correlations for different ocean microseismic sources, for a given position of station pairs. If the observation of discontinuity depth is independent of source location, this would mean that the effect of the source is minimal.
<figure><img src="image_12.png" /><figcaption>**Figure A1.** Source and receiver side rebound on the 410 km discontinuity. Ray tracing zoom on the source (left) and receiver (right) side for $\Delta$ ranging from 150° to 153° and $x = 2°$. The four sources are represented by stars, emitting PKP (orange), PKPPv410p (blue) and Pv410pPKP (purple) phases. The orange station is the arrival point of the PKP phase, the blue/purple station is the arrival point of PKPPv410p and Pv410pPKP phases.</figcaption></figure>