# Imaging with seismic noise: improving extraction of body wave phases from the deep Earth through selective stacking based on H/V ratios

Helle A. Pedersen, Flavien Mattern, Piero Poli, Laurent Stehly

► To cite this version:

Helle A. Pedersen, Flavien Mattern, Piero Poli, Laurent Stehly. Imaging with seismic noise: improving extraction of body wave phases from the deep Earth through selective stacking based on H/V ratios. Geophysical Journal International, 2023, 232 (2), pp.1455-1467. 10.1093/gji/ggac388. hal-04968408

HAL Id: hal-04968408

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

Submitted on 26 Feb 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
# Imaging with seismic noise: improving extraction of body wave phases from the deep Earth through selective stacking based on *H/V* ratios

Helle A. Pedersen,*, Flavien Mattern, Piero Poli and Laurent Stehly

Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTERre, 38000 Grenoble, France.

E-mail: helle.pedersen@univ-grenoble-alpes.fr

Accepted 2022 September 30. Received 2022 September 30; in original form 2022 May 3

## SUMMARY

Generating high-resolution images of the deep Earth remains a challenge. Body waves extracted from noise correlations hold high promise to complement earthquake-based studies, but data processing and interpretation are still under development. We develop a methodology to improve signal-to-noise ratio (SNR) of *P410P* and *P660P*, waves reflected at the top and bottom of the mantle transition zone, using data from the greater Alpine area and focussing on the second microseismic peak (2.5–10 s period). Rather than stacking all available data, we only stack correlations for days with a low ratio of amplitudes between the horizontal plane and vertical direction (*H/V*). Due to an improved SNR we can stack over fewer correlation pairs, with the result that horizontal resolution is significantly improved. We propose a systematic approach to determine at each study point the optimal combination of station pairs and the *H/V* threshold. We observe that the optimal choice of parameters is location dependent and that it is generally different for *P410P* and *P660P*. Additionally, we show that in our study area the maximum interstation distance needs to be reduced to ~150 km for *P410P* to avoid that this arrival is contaminated by surface waves. Applied to the greater Alpine area we demonstrate a significant improvement of signal extraction: while *P410P* and *P660P* were only sporadically identified in standard stacks, with the new processing scheme these arrivals are clearly identified with coherent phases across large distances. We also show that amplitudes of *P660P* decrease drastically around longitude ~11°E to ~12°E, indicating that the lower discontinuity of the transition zone in that area is too broad to have a significant reflexion coefficient for *P* waves in the second microseismic peak.

**Key words:** Composition and structure of the mantle; Europe; Tomography; Body waves; Seismic interferometry; Seismic noise.

## 1 INTRODUCTION

Surface wave imaging based on noise correlations (Shapiro & Campillo 2004; Sabra et al. 2005a) has become standard, also in monitoring temporal evolutions (Brenguier et al. 2008a,b). The surface waves extracted from noise correlations can also be used, for example, to characterize anisotropy (e.g. Lin & Ritzwoller 2011) and the coda quality factor (Soergel et al. 2020). These methods are overall robust, with efforts ongoing to clarify choices of data processing and possible biases due to uneven source distributions (e.g. Pedersen et al. 2007; Fan & Snieder 2009; Froment et al. 2010; Sager et al. 2018).

The use of noise correlations to detect body waves for imaging the deep Earth has proven much more challenging. It is possible to extract body waves at a global scale using cross-correlation, both at

long periods (Nishida 2013) and shorter periods, in particular those corresponding to the second microseismic peak (Boué et al. 2013). The interpretation of global correlations is particularly difficult. First, major earthquakes strongly contribute to the emergence of some of the seismic phases (e.g. Lin & Tsai 2013; Boué et al. 2014; Xia et al. 2016; Poli et al. 2017), with resulting spurious phases and timing bias. The strong influence of earthquake coda on noise correlations has encouraged the development of new seismic imaging methods based on the cross correlation of the coda waves of major earthquakes (for a recent review on these methods, see Tkalčić et al. 2020). Additionally, even if major earthquakes are excluded, the interpretation of the deep phases implies disentangling cross terms from properly reconstructed parts of the Green's function (Pedersen & Colombi 2018; Li et al. 2020a). Specific processing might alleviate such problems (e.g. Spica et al. 2017). Another alternative to
standard cross-correlation is to exploit specific network and noise
source geometries to illuminate target areas (e. g. Poli et al. 2015;
Retailleau et al. 2020; Zhang et al. 2020). It may even be possible
to use seismic waves from strong storms much in the way that seis-
mologists use earthquake data, as suggested by for example Zhang
et al. (2010) and Fan et al. (2019), even though it is challenging to
predict which source areas are efficient to extract specific seismic
phases (Li et al. 2020b).

The discontinuities associated with the mantle transition zone, located at approximately 410 and 660 km depth (depths from the Iasp91 model, Kennett 1991) are key targets for understanding mantle dynamics. The vertical extension and depth of each of the phase transitions are indicators of local temperature and mineralogical composition. Pioneering work (Shearer 1991) demonstrated that weak phases related to the mantle discontinuities can be extracted by stacking many seismograms. For example, topside reflexions of earthquake phases and SS precursors from underside reflections on the mantle discontinuities have small amplitudes, but they can be aligned and stacked so as to constrain lateral variations in arrival times of the two discontinuities, see for example Shearer & Buehler (2019) and references within. Surface wave overtones can also be used to constrain lateral variations of the mantle transition zone (e.g. Meier et al. 2009). In this study, we target the European Alps, beneath which there are strong lateral 3-D variations in the mantle structure over small lateral distances. We therefore need to apply or develop methods that can provide information on lateral variations across distances of tens of kilometres. Receiver functions from teleseismic events may provide at least part of the answer, for example, Liu et al. (2022) used this technique to image the two mantle discontinuities in the Western Alps, based on permanent networks and two dense temporary deployments.

Here we explore the opportunities offered by seismic noise. Sub-
vertical P-wave reflections in seismic noise correlations can offer
an attractive alternative when a dense network is available in the
study area. This technique has the advantage that, at least in theory,
the reflection points are located beneath the central point between
each station pair. Poli et al. (2012) showed that these reflections,
*P410P* and *P660P*, can be extracted from noise correlations. So
far, only a handful of studies have successfully managed to use this
technique for imaging purposes (Poli et al. 2012; Feng et al. 2017,
2018, 2021; Li et al. 2019) in areas with lateral variations over
much larger scale than those in the Alps. Our ambition is to use this
technique to image the two mantle discontinuities over the whole
Alpine area, benefiting from the AlpArray temporary deployment
(AlpArray Seismic Network 2015) which, together with the perma-
nent networks in the area, provides unprecedented data coverage of
the Alps.

Initial tests on our data set showed that the technique proposed by Poli et al. (2012) did not produce reliable *P410P* or *P660P* signals in most locations across our study area, unless the area over which the signals were stacked had a lateral extension which would render geodynamic interpretations meaningless. Therefore, improving the technique proposed by Poli et al. (2012) is a necessity in our study area. Such improvements could potentially make the technique more generally applicable. In parallel with these initial tests, Lu et al. (2022) showed that seasonal variations in the ratio of horizontal to vertical component amplitude (*H*/*V*) in the second microseismic peak of the noise in Europe can be explained by the relative influence of surface waves generated by noise sources in the north Atlantic Ocean and body waves from noise sources in the Southern Hemisphere.

The aim of this work is to develop an approach for body wave
extraction from seismic noise based on *H*/*V* ratios, which can be
easily adapted to a variety of applications. We apply the method to
the European Alps and demonstrate that the extraction of *P410P*
and *P660P* can be significantly improved by stacking correlations
only from those days where the *H*/*V* ratio is low compared to the
local average. A better signal extraction might solve problems of
spurious phases close to the arrival time of *P410P* or *P660P* from
regional distance areas of high noise generation (e.g. Pedersen &
Colombi 2018), but will not solve the issue of spurious phases in
general. Signal extraction is however a prerequisite for imaging,
independently of the imaging technique used. We believe that the
approach that we develop here is sufficiently generic to be adapted
also to the extraction of other body wave phases from the deep
Earth.

## 2 DATA

### 2.1 Waveform data set, initial processing and correlation data set

We base our analysis on the preprocessed data set downsampled to
1 s that Lu et al. (2022) used for their study of the characteristics
of the seismic wavefield in Europe. This data set covers the period
2011 (January 1st)–2019 (October 31st). This data set is based on
publicly available data (through EIDA and through the individual
EIDA nodes, see http://www.orfeus-eu.org/data/eida/). We use data
from permanent stations within the greater Alpine region as our
analysis is dependent on the availability of continuous daily records
from each station pair over several years. We include only one
temporary network: AlpArray (network code Z3; AlpArray Seismic
Network 2015; Hetényi et al. 2018) as it was in operation for several
years with each station operating for at least 3 yr.

To decrease the contributions of earthquakes and oceanic storms, each daily record (vertical component only) was first normalized in the time domain in several period bands using the same comb filter as Soergel et al. (2020) and Nouibat et al. (2022). Taking into account the sampling frequency of the signals and the focus on the second microseismic peak, the effective period range of this study is 2.5–10 s; we apply a fourth order zero phase bandpass (2.5–10 s) Butterworth filter to all correlations before further processing. We correlated all available 4 hr time windows (Z–Z correlations) and normalized and summed the up to six correlations for each day for each station pair. Each such correlation was folded and stored for further analysis in the time window 0–300 s. We used only station pairs with an interstation distance of up to 200 km and the midpoint between the stations located between latitudes 44°N–50°N and longitudes 2.5°E–17.5°E. With these criteria, the correlation data set has a volume of only 41 Gbytes, making it possible to perform many tests to achieve an optimized recovery of *P410P* and *P660P*. The networks used are listed in Section Data Availability.

Fig. 1 shows the station distribution as well as the station mid-
points that are used in this study. We use only station pairs for
which the interstation distance is less than 200 km and for which
we have more than 300 daily correlations. During this study we
perform tests on correlation stacks using station pairs with mid-
points located within circular areas of varying radius around equally
spaced points which we will name *SCP*, for Stack Centre Point.
Three such examples are shown in Fig. 1. Around each of the
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Stations used in this study (white circles). (b) Station pair midpoints (black dots) and Stack Centre Points (SCP, yellow crosses) for correlation stacks. Only midpoints used in this study are shown, that is, midpoints with a distance between stations of up to 200 km and for which there were at least 300 daily correlations. Blue crosses indicate SCP for a section of the Alps which will be referred to subsequently. The three areas (Zones 1–3) used as examples for *P410P* and *P660P* identification are shown. The zones are defined as circles with radius of up to 100 km and containing up to 600 midpoints. Zone 1 (green) has 600 midpoints and a radius of 48.8 km, Zone 2 (blue) has 600 midpoints and a radius of 46.3 km and Zone 3 (red) has a radius of 100 km and 445 midpoints.</figcaption></figure>

three *SCPs* we defined the smallest circle, up to 100 km radius, that has 600 midpoints within the circle. In one of the examples (Zone 3), only 445 midpoints were available within a radius of 100 km. Note that the centre of mass of the midpoints within a circle may not be located at the *SCP*, that is at the centre of the circle (see Zone 3).

The amount of available data (see Fig. S1) increased rapidly during the installation phase of the AlpArray temporary deployment (AlpArray Seismic Network 2015), but even before then, the number of permanent stations increased significantly in Europe and both data quality and access increased for the existing ones. The spatial variation in station density influences the density of midpoints between station pairs used in the study. We illustrate this in Fig. S2 which shows the variations in circle radius if we wish to include 600 midpoints within a circle around each *SCP*, using station pairs for which we have more than 300 daily correlations.

## 3 P410P AND P660P EXTRACTION USING CONTINUOUS NOISE RECORDS: STANDARD METHOD AND DEFINITION OF SIGNAL-TO-NOISE RATIO

Previous work (Poli et al. 2012) demonstrated that *P410P* and *P660P* can be successfully extracted by stacking the traces with a constant slowness, across all times, and identifying high-amplitude arrivals at the times of *P410P* (e.g. 105.85 s for iasp91, Kennett 1991) and *P660P* (e.g. 152.01 s for iasp91). The first outcome of this processing is a vespagram, that is a time-slowness diagram. This method is standard in seismology (see for example Gu & Sacchi 2009) and has fast numerical implementations making it possible to carry out a wide range of tests.

In this work, we focus on *P410P* and *P660P* extraction. In the method set out in the following section, we will carry out many stacks of daily records, followed by many spatial stacks, on a large data set. Some optimization and simplification of the analysis is therefore required. A first question is whether the *P410P* and *P660P* reconstitution, through stacks of correlations within a small area, are significantly reduced in quality if we stack records from stations at different altitudes, or above different crustal structures. Altitude corrections remain minor: for example a station pair where both stations are at 1 km altitude will have vertical P-wave reflexions delayed by less than 5 per cent of the signal period as compared to a station pair where both stations are at sea level. Using a difference of up to ±10 km Moho depth within an area across which correlations are stacked, *P410P* and *P660P* will be less than 1 s in advance or delayed as compared the average arrival times for the area, to be compared with a signal period of 7 s. For these reasons, a first simplification of the processing is to not apply such corrections. A second question is whether correcting for hyperbolical (normal) move out prior to linear stack improves the quality of signal extraction. Initial tests showed no improvement for the distance range and period interval we consider.

Using the Poli et al. (2012) technique, we obtain the E–W seismic section in Fig. 2. We observe that the quality of *P410P* and *P660P* extraction varies strongly within the study area. We also observe that the *P410P* and *P660P* cannot be identified reliably in many locations, making it unrealistic to interpret these seismic sections in terms of Earth structure, either in arrival time or in relative amplitude of *P410P* and *P660P* arrivals. The potential for improvement is illustrated, for example, by the coherent arrivals in some areas (e.g. for *P410P* between 12°E and 15°E, and for *P660P* between 6.5°E and 10.5°E). It is not easy to conclude on the absence in certain areas of *P410P* and *P660P*, their variability is largely independent of the number of station pairs or the radius used for each *SCP*. To quantify the impact of our improved stacking method and the different stacking choices, we first define an estimate of the signal-to-noise ratio (SNR).

Fig. 3 shows the method of the calculation of SNR using vespagrams, with examples from Zones 1, 2 and 3 (for locations see Fig. 1): the time-slowness space is divided into cells of a width (slowness range) of 0.05 s km⁻¹ and height (time range) of 13 s, that is approximately twice the dominant period of the second microseismic peak. We calculate the SNR as the ratio of the sum of the total energy (squared amplitudes) in the target cell, that is the cell containing the theoretical arrival time (*P410P* or *P660P*, using the global model iasp91, Kennett 1991), and the mean energy in the eight cells surrounding the target cell. This definition of the SNR gives more stable values than those obtained for example on the
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Example of a seismic section of zero-slowness stacks using SCPs on latitude 46.5°, obtained using the method of Poli *et al.* (2012). (a) SCPs used, evenly spaced at latitude 46.5°N. (b) Number of midpoints (in blue) and radius *R* (in red) used for the stack. (c) Seismic section corresponding to zero-slowness stack for each SCP. Only midpoints located less than 100 km from the SCP and for which there was at least 300 daily correlations were taken into account. The arrows at times ~100 and ~152 s correspond to the expected arrival time of *P410P* and *P660*, using Iasp91 (Kennett 1991).</figcaption></figure>

zero-slowness stack, due to the averaging across the surrounding cells. In the three examples of Fig. 3, the SNR of *P410P* is 1.66 (Zone 1), 1.07 (Zone 2) and 4.70 (Zone 3) and the SNR of *P660P* is 4.70 (Zone 1), 1.15 (Zone 2) and 0.92 (Zone 3). The SNR captures well when *P410P* and *P660P* are poorly defined within the target cell and/or have small amplitudes as compared to neighbouring amplitude peaks. Note that minor shifts in the location of the target cells did not significantly change SNR. Fig. 3 demonstrates that within two very similar network geometries (Zones 1 and 2), the quality of the *P410P* and *P660P* identification can be very different, without it being clear what causes those differences in quality. Similarly, Zone 3, which has fewer midpoints and a larger radius, has the clearest instance of *P410P* arrival within the three examples, but no clear *P660P*.

## 4 IMPROVING *P410P* AND *P660P* EXTRACTION THROUGH SELECTIVE STACKS BASED ON *HIV*

The key target in this work is to improve extraction of *P410P* and *P660P* from noise correlations. From a theoretical point of view, assuming an isotropic distribution of random noise sources, the SNR of the reconstructed Green's function increases as the square root of the duration of the seismic noise record used to compute the correlation (Weaver & Lobkis 2005; Sabra *et al.* 2005b). With a given SNR of the input signals (correlations), a similar behaviour is expected when the number of seismic stations is increased for the *P410P* and *P660P* signals.

The workflow to improve the extraction *P410P* and *P660P* for a given SCP is the following:

(i) Step 1: Calculate one correlation for each station pair by stacking correlations over an optimal subset of days with available data, for the period 2011–2019.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Vespagrams using standard processing for three examples (geographic locations shown in Fig. 1). Top (a)–(c) Zone 1; Middle (d)–(f) Zone 2; Bottom (g)–(i) Zone 3. Left (a), (d) and (g) show the input correlations, stacked in 3 km bins and normalized separately for times before and after *t* = 70 s. Centre (b), (e) and (h) vespagrams (amplitude of radon transform) between times 70 and 200 s, and slowness between –0.1 and 0.1 s km⁻¹. The target cells for *P410P* and *P660P* are indicated by the white rectangles. The cells used for the SNR calculation are the 8 cells around the target cell. Right (c), (f) and (i) zero slowness stack of seismograms between times 70 and 200 s, with black arrows indicating theoretical arrival times for *P410P* and *P660P* (Iasp91, Kennett 1991). Only midpoints located less than 100 km from the SCP with at least 300 daily correlations were taken into account.</figcaption></figure>
(ii) Step 2: Calculate the vespagram and SNR for each SCP using radon transform (Radon 1917). For this stack we use station pairs for which the interstation distance is in a given range [0-Dmax], and the midpoint is within a given radius R from the SCP.

The choices in Steps 1 and 2 will jointly determine how much data enters the *P410P* and *P660P* extraction, so the optimal stack will result from a combination of these choices, on which there are some additional constraints. First, increasing *R* deteriorates lateral resolution. We therefore cannot increase *R* to very large values in Step 2, even though we observe (not shown here) that increasing *R* leads to improved detection of *P410P* and *P660P*. Using the method of Poli *et al.* (2012) leads, for most SCPs, to such large values of *R* that meaningful geodynamic interpretations beneath the Alps are impossible, given the 3-D structural variations that occur in the area.

The scope for this study is to improve *P410P* and *P660P* extraction so that the lateral resolution is useful in terms of interpretation. The aim of this section is therefore to find a procedure to optimize the three following choices for our study area by considering three parameters:

(i) **Dmax** (maximum interstation distance). Increasing *Dmax* will increase the number of potential correlations to stack, but carries the risk of interference with surface waves or surface wave coda within the time window of *P410P* and *P660P*.

(ii) **Radius R** around a given SCP defining the correlations with midpoints within the circle over which the radon transform is calculated. Increasing *R* will increase the number of correlations to stack and deteriorate lateral resolution.

(iii) **Number of days to stack** based on *H/V*. Increasing the acceptable *H/V* level increases the number of days to stack, see discussion associated with Fig. 5. Note that the chosen days can have a minor impact on *R*, because of data coverage variations over time.

The non-uniform station coverage of the study area significantly complicates this optimization task. For a start, the overall station density varies over time, with a strong increase in some areas after 2015. Secondly, the coverage of permanent stations for which data are available has significantly increased over the last decade. Thirdly, the station density (hence midpoint density) is spatially highly variable between the external and internal parts of the Alps (see Fig. 1). Finally, the number of days of data available for permanent stations is very variable across the networks, because of network specific instrumentation and maintenance issues over time. We therefore accumulate several difficulties with spatial and temporal variations of station distribution. While the parameter choices will depend on the station distribution, the amount of available data and the tectonic setting, we believe that the method we set out in this work can be transposed to other study areas as well as to other body wave arrivals than those considered here.

As a prerequisite, we will accept to optimize the parameters for *P410P* and *P660P* independently, rather than applying the same set of choices for both. This means that we favour interpretations based on arrival times, while the relative amplitudes of *P410P* as compared to *P660P* are meaningless because the stacks are not based on exactly the same data. It is however possible to adapt our strategy to find identical parameters for *P410P* and *P660P*, using a compromise in terms of SNR of the two arrivals, this is shown at the end of this study.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Method to choose *Dmax* for *P410P* and *P660P*: Density plot of SNR for *P410P* (a) and *P660P* (b) as a function of *Dmax*, using circular areas with a radius ensuring 600 midpoints within the circle. We cap *R* at 100 km at a few SCPs, where 600 midpoints are not reached at this value of *R*. We include only midpoints for which there are at least 300 daily correlations. To focus on the evolution of SNR as a function of *Dmax* over all SCPs, with the objective of defining the optimal *Dmax*, the SNR as a function of *Dmax* for each SCP was normalized prior to calculating the density plot. The colour scale indicates the percentage of normalized SNR curves that go through the pixel.</figcaption></figure>

## 4.1 Choosing the optimal maximum interstation distance *Dmax*

A maximum interstation distance *Dmax* of 200 km has routinely been chosen in previous studies (Poli *et al.* 2012; Feng *et al.* 2017; Li *et al.* 2019). This value of *Dmax* means that the surface waves and *P410P* should be separated in time. In the case of the greater Alpine region, the vespagrams (Fig. 3) hint at strong arrivals immediately prior to *P410P*. To check whether these arrivals are linked to surface wave coda we extended the vespagram to slowness values relevant for surface waves, revealing significant remaining energy at late times. This observation is coherent with the study of Soergel *et al.* (2020) who demonstrated that late surface wave coda in noise correlations can be used for estimating coda-Q. The late surface wave arrivals have two (linked) causes: (i) the presence of major sedimentary basins (e.g. Po Plain), characterized by low velocities, and (ii) strong lateral heterogeneities in the study area, which create scattering, hence late arrivals out of which a small proportion will be body waves. We therefore first investigate the choice of *Dmax* on SNR. The method we use to find the optimal value of *Dmax* can easily be applied to other studies.

Fig. 4 shows how SNR evolves as a function of *Dmax*. We observe that SNR for *P410P* reaches a plateau for *Dmax* between 150 and 160 km, with a decrease of SNR for a significant number of correlations beyond that distance. The decrease of SNR beyond this distance is linked to the surface wave coda, as discussed above. On the contrary, SNR for *P660P* increases up to 190–200 km distance, indicating that there is less or no influence of surface waves at these later times, 50 s after the *P410P* arrival. At 200 km distance, surface waves have an apparent velocity of 2 km s⁻¹ at the arrival time of *P410P*, and of 1.3 km s⁻¹ at the arrival time of *P660P*. This difference is apparently enough to significantly reduce the influence of the surface wave coda on *P660P*. Using this approach we can objectivize the choice of *Dmax* = 155 km for *P410P* and *Dmax* = 200 km for *P660P*. Our results indicate that in cratons, where *Dmax*
<figure><img src="image_6.png" /><figcaption>**Figure 5.** (a) Average *H/V* ratio from Lu *et al.* (2022) for network CH (Swiss Seismological Service (SED) at ETH Zurich 1983). (b) Cumulative number of days and percentage of days with *H/V* below a given value. The cut-off values depending on the number of chosen days are indicated as dotted lines.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Main plot: SNR of *P410P* as a function of the number of selected *H/V* days (based on *H/V* as shown in Fig. 5). SNR is calculated using four circles of decreasing size, with 600, 400, 200 and 100 midpoints for Zone 1 (geographic location shown in Fig. 1). For this zone, the radius of the circle decreases by approximately 10 km between successive numbers of midpoints. The vespagrams above the main plot correspond to 600 midpoints and four different choices of *H/V* days, shown as white circles in the main plot. The vespagrams to the left of the main plot are calculated for the four different number of midpoints and for a selected number of *H/V* days equal to 200, represented by white squares in the main plot.</figcaption></figure>

= 200 km has been shown to be a useful choice for *P410P*, it might even be possible to increase *Dmax* beyond 200 km for *P660P*.

## 4.2 Using *H/V* to choose the optimal combination of days and radius around each SCP: principles

Most ambient seismic noise techniques are based on stacking as much data as possible for a given target. Beam forming (Rost & Thomas 2002) or double beam forming (e.g. Nakata *et al.* 2016) have proven successful for extracting specific waves, therefore improving SNR. In the case of *P410P* and *P660P*, directional beams

are not possible because the waves have close to zero slowness, and the direction is therefore not determined. The equivalent to beam forming for *P410P* and *P660P* is therefore non-directional slowness stacks (Fig. 3), which result in the seismic section of Fig. 2.

We here take an approach which, rather than increasing the amount of data to stack, selects a subset of data by selecting time windows (daily correlations) that contribute the most to the *P410P* and *P660P* reconstitution. This selective stacking is based on *H/V* ratios. The absolute value of *H/V* is dependent on the geological structure (Nakamura 1989). Lu *et al.* (2022) showed that variations of *H/V* over time in Europe can be explained by the relative strength of surface waves generated in the Northern Hemisphere
<figure><img src="image_8.png" /></figure>

**Figure 7.** SNR in Zones 2 (top panel) and 3 (bottom panel) as a function of number of stacking days and number of station pairs (capped at 600, and R capped at 100 km). To the left, (a) and (d), is shown how the radius varies with the number of pairs, in the case the number of days being 200 days and 500 days. The centre and right columns show SNR. The white circle indicates the combination with the highest SNR, with its value given below. The red dotted lines delimit areas with SNR > 3. (b) SNR for *P410P*, Zone 2; (c) SNR for *P660P*, Zone 2; (e) SNR for *P410P*, Zone 3; (f) SNR for *P660P*, Zone 3. Note that for areas with low mid point density (fewer seismic stations), it is not always possible to reach 600 midpoints with value of R up to 100 km. This is indicated by the red line in plots (a) and (d), and lead to identical values of SNR in some cells (top left-hand corner of b).

(mainly North Atlantic) and subvertical *P* waves from the Southern Hemisphere. Therefore, *H/V* can potentially be used as a proxy to identify the ‘most’ adequate days for stack, that is those with a significant amount of subvertical *P* waves in the noise as compared to surface waves.

Data for the CH network in Switzerland [Swiss Seismological Service (SED) at ETH Zurich 1983] are available over the whole of 2011–2019, and this network is located centrally in the greater Alpine region. We therefore use the mean *H/V* across the CH network (broad-band stations), calculated by Lu *et al.* (2022), to choose which days to stack. These daily *H/V* ratios at each station are obtained by first filtering between 2.5 and 10 s, then calculating *H/V* (square root of energy in the horizontal plane over the energy in the vertical plane) in 5 min nonoverlapping windows, calculating the median value for the day, and finally averaging over all of the CH network. Further detail can be found in Lu *et al.* (2022). Fig. 5 shows the mean variations of *H/V* for CH between 2011 and 2019. As indicated by Lu *et al.* (2022), *H/V* overall decreases in summer, where storm (wave) activity in the Northern Hemisphere decreases (decreasing the amount of surface waves in the signal) and storm activity in the Southern Hemisphere increases (increasing the amount of subvertically propagating *P* waves).

Initial tests showed that *P410P* and *P660P* identification in the Alps improved if we stacked correlation time windows for June-August rather than for all months; this preliminary result was the original motivation for the present study. Fig. 5 shows that small values of *H/V* only occur during a subset of days during the summer. The strategy of this study is therefore to choose those days where *H/V* is small, but with the *H/V* threshold set to have enough time windows to stack for the *P410P* and *P660P* signals to emerge. *P410P* and *P660P* only emerge when correlations are also stacked over space. There are combinations of stacking days and number of station pairs for which *P410P* and *P660P* emerge clearly, while other combinations will not give satisfactory signal emergence. The aim is here to optimize the number of days to stack so that the *P410P* and *P660P* emerge with a good SNR, whilst keeping a small radius *R* around the SCP.

## 4.3 Using *H/V* to choose the optimal combination of days and radius around each SCP: implementation

Fig. 6 shows an example (*P410P* for Zone 1) of SNR as a function of chosen *H/V* days and for four different circle radii, with 100, 200, 400 and 600 midpoints within. We use constant steps in number
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Seismic sections for optimal stacks across the Alps, using the layout of Fig. 2. Left (d–f) show results for *P410P* and right (g–i) for *P660P*. For easy comparison, (a)–(c) show the outcome of standard stacks (Fig. 2). The number of station pairs (blue line) used in the stacks, and the associated value of *R* (red line), are shown above each seismic section. The arrows at times ~100 and ~152 s correspond to the expected arrival time of *P410P* and *P660P*, using Iasp91 (Kennett 1991).</figcaption></figure>

of midpoints rather than constant steps in *R*, to have comparable amounts of data within the stacks for different SCPs. We observe that SNR in general decreases (for 600 midpoints) with an increasing number of days, but the curves have no simple behaviour. In this precise case, the reduction of midpoints from 600 to 400 (reducing the circle radius *R* from approximately 50 km to approximately 40 km) does not significantly decrease SNR when the number of stacked days based on *H/V* is well chosen. Similar plots for SCPs with fewer midpoints show that with a careful choice of stacking, the improvement of resolution can be significant, even for areas with sparse station density. For a same location but for *P660P* (see

Supplementary Material S3), the SNR curves are quite different, illustrating that it can be useful to have different stacks for two reflectors, even at the same location.

It is possible to generalize this approach to a systematic optimization of the parameter choice, by calculating SNR as a function of number of stacking days and number of midpoints (Fig. 7). The combination of *R* and number of days to stack that gives the highest SNR are indicated by a circle in each plot, we will refer to this as the optimal stack. We observe that in a given area, the optimal stack for *P410P* and *P660P* correspond to different values of *R* and number of days. Similarly, for a given reflector, the optimal stack parameters
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Seismic section (f) using for the stack, for each SCP, a combination of number of stacking days and number of midpoints with SNR > 3 for both *P410P* and *P660P*, and using *Dmax* = 155 km. The stack for *P410P* and *P660P* is therefore identical (Choice 3). The number of station pairs (blue line in e) used in the stacks, and the associated value of R (red line), are shown above the seismic section. The arrows at times ~100 and ~152 s correspond to the expected arrival time of *P410P* and *P660P*, using IASP91 (Kennett 1991). For easy comparison, (a)–(c) show the outcome of standard stacks (Fig. 2).</figcaption></figure>

are different in different geographical areas. The figure also shows
that there are different possible parameter choices depending on
scientific objective, for example:

(i) **Choice 1:** *R* and number of days chosen separately for *P410P* and *P660P* so that SNR is maximum for each reflector. This corresponds to the optimal stacks (white circles) of Fig. 7.

(ii) **Choice 2:** *R* and number of days chosen separately for *P410P* and *P660P* so as to use the smallest possible value of *R* for which a given SNR threshold is reached for each reflector, for example SNR > 3. If SNR > 3 is not reached, we use the combination of *R*

<figure><img src="image_11.png" /><figcaption>**Figure 10.** SNR for *P410P* and *P660P* using an optimal stack designed to maximize the SNR (Choice 1). The colour code indicates SNR, and the circles (SCPs) are relocated to the centre of mass of the station geometry for the optimal stack in each SCP. (a) *P410P* and (b) *P660P*.</figcaption></figure>

and number of days as defined in Choice 1, that is use the optimal
stack.

(iii) **Choice 3:** *R* and number of days chosen jointly for *P410P* and *P660P* so as to use the smallest possible value of *R* for which a given SNR threshold is reached for both reflectors. With this choice it is possible to compare relative amplitudes between *P410P* and *P660P*, assuming that an identical value of *Dmax* has been chosen in the previous step. If SNR > 3 is not reached for both *P410P* and *P660P* we choose to normalize the SNR analysis for each of them and calculate the parameter combination which gives the highest mean SNR for *P410P* and *P660P*. In this way both of them have an equal weight in the choice of stacking parameters.

We will show example of cross-sections resulting from each of
these three choices.

4.4 Seismic sections using different stacking choices

The same seismic section as the one in Fig. 2, but this time using
optimal stack (Choice 1) is shown in Fig. 8. In this case the stacks for
*P410P* and *P660P* are optimized separately, leading to two seismic
sections (Figs 8f and i). By comparison between these optimized
sections and standard stacks (for ease of comparison, included as
Fig. 8c) the improvement of the quality of the signal is spectacular.
At the same time, resolution is improved, that is with values of *R*
(Figs 8e and h) which are significantly smaller than those (Fig. 8b)
used in the standard stacks. We also observe that the optimization
parameters are not the same for the two discontinuities.
<table><caption>**Table 1.** Network codes and references for seismic networks used in this study.</caption><thead><tr><th>Network Code</th><th>Reference</th></tr></thead><tbody><tr><td>BW</td><td>Department of Earth and Environmental Sciences, Geophysical Observatory, University of Munchen (2001)</td></tr><tr><td>CH</td><td>Swiss Seismological Service (SED) At ETH Zurich (1983)</td></tr><tr><td>CR</td><td>University of Zagreb (2001)</td></tr><tr><td>CZ</td><td>Charles University in Prague (Czech), Institute of Geonics, Institute of Geophysics, Academy of Sciences of the Czech Republic, Institute of Physics of the Earth Masaryk University (Czech) and Institute of Rock Structure & Mechanics (1973)</td></tr><tr><td>FR</td><td>RESIF (1995)</td></tr><tr><td>G</td><td>Institut De Physique Du Globe De Paris (IPGP), and Ecole Et Observatoire Des Sciences De La Terre De Strasbourg (EOST) (1982)</td></tr><tr><td>GE</td><td>GEOFON Data Centre (1993)</td></tr><tr><td>GR</td><td>Federal Institute for Geosciences and Natural Resources (1976)</td></tr><tr><td>GU</td><td>University of Genoa (1967)</td></tr><tr><td>HU</td><td>Kövesligethy Radó Seismological Observatory [Geodetic And Geophysical Institute, Research Centre for Astronomy and Earth Sciences, Hungarian Academy of Sciences (MTA CSFK GGI KRSZO)] (1992)</td></tr><tr><td>IV</td><td>INGV Seismological Data Centre (2006)</td></tr><tr><td>MN</td><td>MedNet Project Partner Institutions (1990)</td></tr><tr><td>MT</td><td>French Landslide Observatory—Seismological Datacenter/RESIF (2006)</td></tr><tr><td>NI</td><td>OGS (Istituto Nazionale di Oceanografia e di Geofisica Sperimentale) and University of Trieste (2002)</td></tr><tr><td>RD</td><td>RESIF (2018)</td></tr><tr><td>RF</td><td>University of Trieste (1993)</td></tr><tr><td>SI</td><td>Province Südtirol network (ZAMG)</td></tr><tr><td>SK</td><td>ESI SAS (Earth Science Institute of the Slovak Academy of Sciences) (2004)</td></tr><tr><td>SL</td><td>Slovenian Environment Agency (1990)</td></tr><tr><td>ST</td><td>Geological Survey-Provincia Autonoma di Trento (1981)</td></tr><tr><td>Z3 (2015)</td><td>AlpArray Seismic Network (2015)</td></tr></tbody></table>

Another stacking option would be to favour the smallest possible value of *R*, to optimize resolution. In that case one can use the combination of parameters that gives SNR larger than a given threshold with the smallest possible value of *R* (Choice 2), for each of *P410P* and *P660P* separately. With this choice, and using the threshold SNR > 3 (see Fig. S4), *R* is less than 50 km in most locations. Again, the stacks are significantly better than the standard stacks of Fig. 2.

Finally, it is possible to use a strategy (Choice 3) where the stacks are identical for *P410P* and *P660P* (same traces and number of days used) and still significantly improved. With this choice of processing, it is possible to compare the amplitudes of the two phases. Fig. 9 shows an example of such processing, where we searched for the minimum value of *R* for which SNR > 3 for both *P410P* and *P660P*. In this example, we used a maximum interstation distance *Dmax* of 155 km.

Care must be taken as to the interpretation of the different seismic sections, as we have not applied any corrections for altitude and crustal thickness. Scaling the study up from seismic sections to delay measurements for each reflector across the study area is the next natural step of this study, including also data from 2020 to

2022. One observation stands out which is independent of such pre-stack corrections: SNR for *P660P* drops east of approximately 11°E, independently of the choice of stack, while SNR of *P410P* is approximately stable along the seismic sections of Figs 8 and S4. This difference of behaviour is confirmed in Fig. 9, where the stack is the same for both reflectors.

Additional information is available by analysing SNR for each discontinuity, as illustrated by Fig. 10. Optimal SNR for *P410P* is particularly poor around the Po Plain of northern Italy, and in random locations in peripheral areas. Overall, SNR of *P410P* is higher in the eastern part of the study area than in the western part. For *P660P* there are also low SNR is some peripheral areas, but the Po Plain has good SNR (possibly due to use of larger *Dmax*, so integrating more Alpine stations). The SNR drops abruptly across longitude ~11°E. This observation is stable independently of the parameter combination we use for the stacks (see also Fig. 9). We therefore suggest that there is an abrupt and significant E–W change of the characteristics of the 660 km discontinuity beneath the Alps.

CONCLUSIONS

In this work, we propose a strategy to optimize the SNR of *P410P* and *P660P* extracted from seismic noise correlations. We observe that an adequate choice of three parameters for each reflector and location can significantly improve the SNR of these two phases as compared to the technique of Poli *et al.* (2012). The key element in our approach is to only stack a subset of days with a low *H/V* ratio, that is days where the amount of subvertically propagating body waves is high in the noise as compared to surface waves from the surrounding seas (Lu *et al.* 2022). With this approach, we can reduce the effective area over which the stacks are performed and improve lateral resolution. We also adapt the maximum distance between station pairs independently for each of *P410P* and *P660P*. The methodology that we use for choosing the three parameters can easily be transposed to other locations, seismic noise characteristics and network characteristics (operating time and network geometry). It may also be applicable to other seismic phases. Our method is flexible and can be optimized for different scientific goals: optimal observation of phases (Choice 1), best lateral resolution (Choice 2) or identical stacks for the two phases (Choice 3). Consequently, the compromises made between signal quality and resolution are explicit and scientifically driven.

In this work, we focused on extracting the body wave reflections from the noise. For geodynamic interpretations, further steps need to be taken. First the correlations should at least be corrected for altitude and crustal structure prior to stack. To transform the arrival times to depth, it would additionally be necessary to correct for lateral variations in mantle velocities. Once these corrections are carried out, the combination of noise based studies and receiver functions from earthquake waves open up for new constraints on the mantle discontinuities. These discontinuities are in reality strong velocity gradients and the thickness of the ‘discontinuities’ significantly influences the P-P reflection coefficient as well as the coefficient of wave conversions. The joint interpretation of noise based (~7 s period) and earthquake based (typically 12–15 s period) data could provide valuable constraints on the lateral variations of temperature and composition in the mantle transition zone. One stable observation in the study area is that for the period range considered here (2.5–10 s), the SNR of the reflections on the 660 km discontinuity strongly decrease in amplitude east of approximately 11–12°E.
Further work needs to be carried out to understand if the reflection points are located underneath the study area. With a reduced number of sources, it may well be that the retrieved reflections in the correlations are not located in the central point between each station pair and that the stacks should be attributed to another geographical location other than the centre of mass of the theoretical reflection points. Independently of these challenges, we argue that selective stacking of correlations can be a fundamental asset in the use of seismic noise to image the deep Earth with body waves.

## ACKNOWLEDGMENTS

We first of all thank all the institutions and personnel involved in operating the seismic networks in Europe, and ORFEUS/EIDA and the EIDA nodes for the data distribution. Special thanks to the many teams involved in the AlpArray Project. We thank Yang Lu for providing the H/V data. Piero Poli and Flavien Mattern were supported from the European Research Council (ERC) under the European Union Horizon 2020 Research and Innovation Program (grant agreement 802777- MONIFaults). All authors are part of Labex OSUG@2020 (ANR10 LABX56). We thank Christian Stanciu and an anonymous reviewer for constructive comments to the initial version of the manuscript.

## DATA AVAILABILITY

This study is based on publicly available data (through EIDA and through the individual EIDA nodes, see http://www.orfeus-eu.org/data/eida/). Table 1 shows the arrays that contribute data to the correlation data set used in this study.).

## REFERENCES

AlpArray Seismic Network, 2015. AlpArray Seismic Network (AASN) temporary component [Data set]. AlpArray Working Group. Other/Seismic Network, https://doi.org/10.12686/alarray/z3_2015.

Boué, P., Poli, P., Campillo, M., Pedersen, H., Briand, X. & Roux, P., 2013. Teleseismic correlations of ambient seismic noise for deep global imaging of the Earth, *Geophys. J. Int.*, **194**, 844–848.

Boué, P., Poli, P., Campillo, M. & Roux, P., 2014. Reverberations, coda waves and ambient noise: correlations at the global scale and retrieval of the deep phases, *Earth planet. Sci. Lett.*, **391**, 137–145.

Brenguier, F., Campillo, M., Hadzioioannou, C., Shapiro, N. M., Nadeau, R. M. & Larose, E., 2008a. Postseismic relaxation along the San Andreas fault at Parkfield from continuous seismological observations, *Science*, **321**, 1478–1481.

Brenguier, F., Shapiro, N. M., Campillo, M., Ferrazzini, V., Duputel, Z., Coutant, O. & Nercession, A., 2008b. Towards forecasting volcanic eruptions using seismic noise, *Nat. Geosci.*, **1**(2), 126–130.

Charles University in Prague (Czech), Institute of Geonics, Institute of Geophysics, Academy of Sciences of the Czech Republic, Institute of Physics of the Earth Masaryk University (Czech), & Institute of Rock Structure and Mechanics. 1973. Czech Regional Seismic Network [Data set]. International Federation of Digital Seismograph Networks, https://doi.org/10.7914/SN/CZ.

Department of Earth and Environmental Sciences, Geophysical Observatory, University of Munchen. 2001. BayernNetz [Data set]. International Federation of Digital Seismograph Networks, https://doi.org/10.7914/SN/BW.

ESI SAS (Earth Science Institute Of The Slovak Academy Of Sciences). 2004. National Network of Seismic Stations of Slovakia [Data set]. Deutsches GeoForschungsZentrum GFZ, https://doi.org/10.14470/FX099882.

Fan, W., McGuire, J. J., de Groot-hedlin, C. D., Hedlin, M. A. H., Coats, S. & Fiedler, J. W., 2019. Stormquakes, *Geophys. Res. Lett.*, **46**, 12 909–12 918.

Fan, Y. & Snieder, R., 2009. Required source distribution for interferometry of waves and diffusive fields, *Geophys. J. Int.*, **179**(2), 1232–1244.

Federal Institute for Geosciences and Natural Resources. 1976. German Regional Seismic Network (GRSN) [Data set]. Bundesanstalt für Geowissenschaften und Rohstoffe, https://doi.org/10.25928/mbx6-hr74.

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**, 8328–8335.

Feng, J., Yao, H. & Wang, W., 2018. Imaging mantle transition zone discontinuities in southwest China from dense array ambient noise interferometry, *Earthq. Sci.*, **31**, 301–310.

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

French Landslide Observatory – Seismological Datacenter /RESIF. 2006. Observatoire Multi-disciplinaire des Instabilités de Versants (OMIV) [Data set]. RESIF - Réseau Sismologique et géodésique Français, https://doi.org/10.15778/RESIF.MT.

Froment, B., Campillo, M., Roux, P., Gouedard, P., Verdel, A. & Weaver, R., 2010. Estimation of the effect of nonisotropically distributed energy on the apparent arrival time in correlations, *Geophysics*, **75**(5), SA85-SA93.

GEOFON Data Centre. 1993. GEOFON Seismic Network [Data set]. Deutsches GeoForschungsZentrum GFZ, https://doi.org/10.14470/TR560404.

Geological Survey-Provincia Autonoma di Trento. 1981. Trentino Seismic Network [Data set]. International Federation of Digital Seismograph Networks, https://doi.org/10.7914/SN/ST.

Gu, Y.J. & Sacchi, M., 2009. Radon transform methods and their applications in mapping mantle reflectivity structure. *Surv. Geophys.*, **30**, 327–354.

Hetényi, G. et al., 2018. The AlpArray Seismic Network: a large-scale european experiment to image the alpine orogen, *Surv. Geophys.*, **39**, 1009–1033.

INGV Seismological Data Centre. 2006. Rete Sismica Nazionale (RSN) [Data set]. Istituto Nazionale di Geofisica e Vulcanologia (INGV), Italy. https://doi.org/10.13127/SD/X0FXnH7QfY.

Institut De Physique Du Globe De Paris (IPGP), & Ecole Et Observatoire Des Sciences De La Terre De Strasbourg (EOST). 1982. GEOSCOPE, French Global Network of broad band seismic stations [Data set]. Institut de physique du globe de Paris (IPGP), Université de Paris, https://doi.org/10.18715/GEOSCOPE.G.

Kennett, B.L.N. (Compiler and Editor). 1991. IASPEI 1991 Seismological Tables. Bibliotech, Canberra, Australia, pp. 167, https://doi.org/10.1111/j.1365-3121.1991.tb00863.x.

Kövesligethy Radó Seismological Observatory (Geodetic And Geophysical Institute, Research Centre For Astronomy And Earth Sciences, Hungarian Academy Of Sciences (MTA CSFK GGI KRSZO)). 1992. Hungarian National Seismological Network [Data set]. Deutsches GeoForschungsZentrum GFZ, https://doi.org/10.14470/UH028726.

Li, L., Boué, P. & Campillo, M., 2020a. Observation and explanation of spurious seismic signals emerging in teleseismic noise correlations, *Solid Earth*, **11**, 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), doi:10.1029/2020GC009112.

Li, S., Huang, J., Liu, Z. & Zhang, Y., 2019. Study on the mantle discontinuity structures beneath Northeast China with time-frequency phase-weighted stacks of ambient noise correlations, *Geophys. J. Int.*, **218**, 1490–1501.

Lin, F.-C. & Ritzwoller, M. H., 2011. Helmholtz surface wave tomography for isotropic and azimuthally anisotropic structure, *Geophys. J. Int.*, **186**(3), 1104–1120.

Lin, F.-C. & Tsai, V.C., 2013. Seismic interferometry with antipodal station pairs, *Geophys. Res. Lett.*, **20**, 4609–4613.
Liu, D. et al., 2022. Receiver function mapping of the mantle transition zone beneath the Western Alps: new constraints on slab subduction and mantle upwelling, *Earth planet. Sci. Lett.*, **577**, doi:10.1016/j.epsl.2021.117267.

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

MedNet Project Partner Institutions. 1990. Mediterranean Very Broadband Seismographic Network (MedNet) [Data set]. Istituto Nazionale di Geofisica e Vulcanologia (INGV), https://doi.org/10.13127/SD/fBBBtDtd6q.

Meier, U., Trampert, J. & Curtis, A., 2009. Global variations of temperature and water content in the mantle transition zone from higher mode surface waves, *Earth Planet. Sci. Lett.*, **282**, 91–101.

Nakamura, Y., 1989. A method for dynamic characteristics estimation of subsurface using microtremor 1151 on the ground surface, *QR Railway Tech. Res. Inst.*, **20**, 1.

Nakata, N., Boué, P., Brenguier, F., Roux, Ph., Ferrazzini, V. & Campillo, M., 2016. Body and surface wave reconstruction from seismic noise correlations between arrays at Piton de la Fournaise volcano, *J. geophys Res.*, **123**, 6911–6921.

Nishida, K., 2013. Global propagation of body waves revealed by crosscorrelation analysis of seismic hum, *Geophys. Res. Lett.*, **40**, 1691–1696.

Nouibat, A., et al., Cifalps Team & Alparray Working Group 2022. Lithospheric transdimensional ambient-noise tomography of W-Europe: implications for crustal-scale geometry of the W-Alps, *Geophys. J. Int.*, **229**(2), 862–879.

OGS (Istituto Nazionale di Oceanografia e di Geofisica Sperimentale) and University of Trieste. 2002. North-East Italy Broadband Network [Data set]. International Federation of Digital Seismograph Networks, https://doi.org/10.7914/SN/NI.

Pedersen, H. A. & 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**, 759–772.

Pedersen, H. A. & Krüger, F., SVEKALAPKO Seismic Tomography Working Group, 2007. Influence of the seismic noise characteristics on noise correlations in the Baltic Shield, *Geophys. J. Int.*, **168**, 197–210.

Poli, P., Campillo, M. & Pedersen, H. A., the POLENET/LAPNET Working Group, 2012. Body wave imaging of the Earth's mantle discontinuities from ambient seismic noise, *Science*, **338**(6110), 1063–1065.

Poli, P., Campillo, M. & van de Hoop, 2017. Analysis of intermediate period correlations of coda from deep earthquakes, *Earth planet. Sci. Lett.*, **477**(2), 147–155.

Poli, P., Thomas, C., Campillo, M. & Pedersen, H.A., 2015. Imaging the D' reflector with noise correlations, *Geophys. Res. Lett.*, **42**(1), 60–65.

Radon, J., 1917. Über die Bestimmung von Funktionen durch ihre Intergralwerte längs gewisser Mannigfaltkeiten, *J. Radon. Akad. Wiss.*, **69**, 262–277.

Research Centre For Astronomy And Earth Sciences, Hungarian Academy Of Sciences (MTA CSFK GGI KRSZO)). 1992. Rete Sismica Nazionale (RSN) [Data set]. Istituto Nazionale di Geofisica e Vulcanologia (INGV), Italy, https://doi.org/10.13127/SD/X0FXnH7QfY.

RESIF 1995. RESIF-RLBP French Broad-band network, RESIF-RAP strong motion network and other seismic stations in metropolitan France [Data set]. RESIF - Réseau Sismologique et géodésique Français, https://doi.org/10.15778/RESIF.FR.

RESIF 2018. CEA /DASE broad-band permanent network in metropolitan France [Data set]. RESIF - Réseau Sismologique et géodésique Français, https://doi.org/10.15778/RESIF.RD.

Retailleau, L., Boué, P., Li, L. & Campillo, M., 2020. Ambient seismic noise imaging of the lowermost mantle beneath the North Atlantic Ocean. *Geophys. J. Int.*, **222**(2), 1339–1351.

Rost, S. & Thomas, C., 2002. Array seismology: methods and applications, *Rev. Geophys.*, **40**(3), 2–12–27.

Sabra, K. G., Gerstoft, P., Roux, Ph., Kuperman, W. A. & Fehler, M. C., 2005a. Surface wave tomography from microseisms in Southern California, *Geophys. Res. Lett.*, **32**(14), doi:10.1029/2005GL023155.

Sabra, K. G., Roux, P. & Kuperman, W. A., 2005b. Emergence rate of the time-domain Green's function from the ambient noise cross-correlation function, *J. acoust. Soc. Am.*, **118**, 3524–3531.

Sager, K., Boehm, C., Ermert, L., Krischer, L. & Fichtner, A., 2018. Sensitivity of seismic noise correlation functions to global noise sources, *J. geophys. Res.*, **123**, 6911–6921.

Shapiro, N.M. & Campillo, M., 2004. Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise, *Geophys. Res. Lett.*, **31**, L07614.

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

Shearer, P. M. & Buehler, J., 2019. Imaging upper-mantle structure under USArray using long-period reflection seismology. *J. geophys. Res.*, **124**, 9638–9652.

Slovenian Environment Agency. 1990. Seismic Network of the Republic of Slovenia [Data set]. International Federation of Digital Seismograph Networks, https://doi.org/10.7914/SN/SL.

Soergel, D. et al., 2020. Coda-Q in the 2.5–20 s period band from seismic noise: application to the greater Alpine area, *Geophys. J. Int.*, **220**(1), 202–217.

Spica, Z., Perton, M. & Beroza, G. C., 2017. Lateral heterogeneity imaged by small-aperture ScS retrieval from the ambient seismic field, *Geophys. Res. Lett.*, **44**, 8276–8284.

Swiss Seismological Service (SED) At ETH Zurich. 1983. National Seismic Networks of Switzerland [Data set]. ETH Zürich, https://doi.org/10.12686/sed/networks/ch.

Tkalčić, H., Pham, T.-S. & Wang, S., 2020. The Earth's coda correlation wavefield: rise of the new paradigm and recent advances, *Earth Sci. Rev.*, **208**, doi:10.1016/j.earscirev.2020.103285.

University of Genoa. 1967. Regional Seismic Network of North Western Italy [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/GU.

University of Trieste. 1993. Friuli Venezia Giulia Accelerometric Network [Data set]. International Federation of Digital Seismograph Networks, https://doi.org/10.7914/SN/RF.

University of Zagreb. 2001. Croatian Seismograph Network [Data set]. International Federation of Digital Seismograph Networks, https://doi.org/10.7914/SN/CR.

Weaver, R. L. & Lobkis, O. I., 2005. Fluctuations in diffuse field-field correlations and the emergence of the Green's function in open systems, *J. acoust. Soc. Am.*, **117**, 3432–3439.

Xia, H. H., Song, X. & Wang, T., 2016. Extraction of triplicated PKP phases from noise correlations, *Geophys. J. Int.*, **205**, 499–508.

Zhang, J., Gerstoft, P. & Shearer, P. M., 2010. Resolving P-wave travel-time anomalies using seismic array observations of oceanic storms, *Earth planet. Sci. Lett.*, **292**(3-4), 419–427.

Zhang, L., Li, J., Wang, T., Yang, F. & Chen, Q.-F., 2020. Body waves retrieved from noise cross-correlation reveal lower mantle scatterers beneath the Northwest Pacific subduction zone, *Geophys. Res. Lett.*, **47**, e2020GL088846, doi:10.1029/2020GL088846.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

**Figure S1.** Evolution of the amount of data for the data set, (a) in number of station pair midpoints during the year for which at least 1 d of correlation is available and (b) in total number of daily correlations available for the year.

**Figure S2.** Illustration (circles not to scale) of the geographical distribution of the maximum radius (SCP—midpoint distance) reached by considering 600 midpoints for each SCP in the study area. Only station pairs for which there were at least 300 daily correlations were taken into account. Note that the the radius of the circles should be multiplied by 6 to be to scale.
**Figure S3.** Main plot: SNR of *P660P* as a function of the number of selected *H/V* days (based on *H/V* as shown in Fig. 5). SNR is calculated using four circles of decreasing size, with 600, 400, 200 and 100 midpoints for Zone 1 (geographic location shown in Fig. 1). For this zone, the radius of the circle decreases by approximately 10 km between successive number of midpoints. The vespagrams above the main plot correspond to 600 midpoints and four different choices of *H/V* days, shown as white circles in the main plot. The vespagrams to the left of the main plot are calculated for the four different number of midpoints and for a selected number of *H/V* days equal to 250, represented by white squares in the main plot.

**Figure S4.** Sections obtained using, for each SCP and independently for *P410P* and *P660P*, a combination of number of stacking days

and number of midpoints with SNR > 3 and for which R is as small as possible (Choice 2). Left (d–f) show results for *P410P* and right (g–i) for *P660P*. The number of station pairs (blue line) used in the stacks, and the associated value of R (red line), are shown above each seismic section (e, h). The arrows at times ~100 and ~152 s correspond to the expected arrival time of *P410P* and *P660P*, using IASP91 (Kennett 1991).

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.