Ambient seismic noise imaging of the lowermost mantle beneath the North Atlantic Ocean

SUMMARY

Body waves can be extracted from correlation functions computed from seismic records even at teleseismic distances. Here we use P and PcP waves from the secondary microseism frequency band that are propagating between Europe and the Eastern United States to image the core–mantle boundary (CMB) and D” structure beneath the North Atlantic. This study presents the first 3-D image of the lower mantle obtained from ocean-generated microseism data. Robustness of our results is evaluated by comparing images produced by propagation in both directions. Our observations reveal complex patterns of lateral and vertical variations of P-wave reflectivity with a particularly strong anomaly extending upward in the lower mantle up to 2600 km deep. We compare these results with synthetic data and associate this anomaly to a Vp velocity increase above the CMB. Our image aims at promoting the study of the lower mantle with microseism noise excitations.

1 INTRODUCTION

Constraining the chemical, thermal and mechanical processes that occur in the deep mantle is critical to better understand the Earth's past and future dynamics. One of the Earth's major discontinuities, the core–mantle boundary (CMB) is conventionally located at 2981 km depth. It marks the limit between the turbulent flow of liquid molten iron alloy in the outer core and the viscous silicate rocks of the lower mantle. It is also, with the Earth's surface, the region that shows the most contrasts in density, viscosity and elastic wave properties. The lowermost mantle, usually defined as the region between the CMB and 200–300 km above it and designated as the D'' layer (Bullen 1949), is known to be a complex and discontinuous shell that plays a significant role in fundamental geodynamic processes such as mantle convection (Lay 2007) and dynamo effect operating in the outer-core (Glatzmaier et al. 1999). The D'' layer is generally associated with a thermal boundary layer produced by heat fluxing from the much hotter core to the colder lower mantle. A simple thermal discontinuity is not the only candidate that could lead to a seismic wave reflection/diffraction and multiple factors and combinations of factors were rapidly proposed in the literature to explain seismic observations such as a presence of sinking slab debris or thermochemical layering (Boyet & Carlson 2005; Hernlund & McNamara 2015). The presence of subducted material in the deep mantle as a possible origin for plumes and hotspots has been debated for several years (Workman & Hart 2005). In contrast, laboratory and computational mineral physics experiments demonstrated the existence of a phase transition from perovskite to a denser post-perovskite at pressure and temperature close to lower mantle conditions (see Cobden et al. 2015, for a review). But even if this strong mineral physics argument seems compatible with most seismic observations, the consensus is moving towards a multifactor explanation. Consequently, multidisciplinary analyses have been carried out for the past decades to obtain reasonable structural and geodynamic models of the deep mantle and outermost core. To construct a realistic picture, one should use all available inputs from each of these fields to deal with inherent tradeoffs and/or uncertainties of each method.

Seismology plays a major role by bringing direct measurements of elastic properties of the Earth structure with depth. Numerous global tomographic inversions have defined the large-scale structures of the lowermost mantle using various seismic data sets: body wave arrival times, normal modes and/or surface wave dispersion (e.g. Mégnin & Romanowicz 2000; Panning & Romanowicz 2006; Houser et al. 2008; Kustowski et al. 2008; Simmons et al. 2009; Ritsema et al. 2011; Lekic et al. 2012; Moulik & Ekström 2014; French & Romanowicz 2015; Durand et al. 2016; Garnero et al. 2016; Koelemeijer et al. 2016; Durand et al. 2017). These studies have highlighted the presence of two Large Low Shear Velocity Provinces (LLSVPs) beneath Africa and the Pacific. Given the input data and the inversion methods, these tomographic models are usually depicting smooth deep Earth structures. The lack of resolution within the lower mantle (generally > 1000 km) prevents consensual conclusions on the fine structures and more importantly on physicochemical processes. Seismic velocities are obviously important attributes for identifying the physical nature of heterogeneities, but velocity contrasts alone are really informative only if the spatial geometry of these heterogeneities is properly constrained across scales. Tomographic images are too smooth to reveal sharp details and other constraints are needed. For instance, better azimuthal coverage of the illumination of the deep mantle (together with a denser ray coverage) helps to measure the anisotropy and better constrain mineralogical structure at depth (Cobden et al. 2015). Observations and modelling of various specific seismic phases in local regions of the globe have been critical to identifying smaller-scale heterogeneities such as Ultra Low Velocity Zones (ULVZ, Williams & Garnero 1996) and small-scale (<10 km) scatterers (e.g. Rost & Earle 2010). Using high frequency scattered/reflected P and S waves is one of the best ways to detect and probe these heterogeneities (e.g. Earle & Shearer 1997). The top of the D” discontinuity has been extensively studied in various places around the globe (e.g. Lay & Helmberger 1983; Thomas et al. 2004a, b; Gassner et al. 2015). Overviews on the lower mantle and D’’ can be found in Wysession et al. (1998), Garnero (2000), Lay & Garnero (2011) and Lay (2015).

Deep-Earth seismology surveys rely on moderate to large earthquakes to generate impulsive and sufficiently energetic seismic waves that can interact with the fine structure of the Earth and finally reach a sensor at the Earth's surface. Seismologists use various types of wave travel-paths to optimize the Earth coverage, but the regions probed by these waves are ultimately controlled and limited by the earthquake distribution. However, we note the significant effort and constant improvement made to exploit the signals from the earthquakes as fully as possible. For instance, Lai et al. (2019) developed a method to build high-quality traveltime data set including multiply reflected waves such as SSS and ScSScS waves. If we look carefully at previous studies made in our region of interest, which is beneath the North Atlantic, the lower mantle was probed by few studies using earthquakes mostly distributed in South and Central America. Closer to our study, Yao et al. (2015) observed clear S wave reflected waves (SdS) related to D’’ under the North Atlantic Ocean from a single earthquake in the South of Spain and measured on an array of sensors in the USA (transportable USArray). They concluded the existence of a ∼300 km thick D” with a positive velocity jump ranging from 2 to 4 per cent, matching previous studies in neighboring regions. Durand et al. (2019) recently extended that study by detecting both SdS and PdP reflections. This latest article also concluded that the D” layer region beneath the North Atlantic Ocean shows a strong positive velocity jump (ranging from 2.7 to 3.8 per cent) at 300 km above the CMB with significant lateral variations. A recent review of previous studies in this region can be found in Li et al. (2019). This kind of studies would benefit significantly from an increase in the number of reflection points at the CMB and above, with a broader offset and azimuthal coverage, which would help constrain the 3-D geometry of the observed discontinuity.

In this study, we focus on specific types of seismological data that are retrieved from ambient seismic noise. At period band longer than 1 s the ambient seismic noise is linked mostly to the interactions between the ocean and the solid Earth (Longuet-Higgins 1950; Hasselmann 1963). The strongest energy, in the 3 to 10 s period band, is called secondary microseism (Kedar et al. 2008; Ardhuin et al. 2011; Hillers et al. 2012; Stutzmann et al. 2012). It is generated by the non-linear interaction between oceanic waves travelling at the same period with opposite directions. This interaction occurs either within a storm, between two storms, or when ocean waves generated by a storm propagate to a coast and are reflected back in the direction of the storm. The resulting standing wave generates a pressure field at the bottom of the ocean that is converted to seismic waves.

The permanent weak and complex wavefield resulting from these oceanic sources can be turned to useful deterministic signals using the correlation operator (seismic interferometry). Noise-based seismic methods constitute a valuable supplement to earthquake data since they are mostly constrained by the sensor locations (Shapiro & Campillo 2004; Sabra et al. 2005; Shapiro et al. 2005). Their application was first limited to the uppermost layers of the Earth, as surface waves dominate the extracted signals. However, it was demonstrated that crustal body-wave reflections (Zhan et al. 2010; Ruigrok et al. 2011; Poli et al. 2012a) can be extracted from ambient noise correlations. Poli et al. (2012b) then made a first observation of the upper mantle discontinuities (at 410 and 660 km of depth) using a similar method across a cratonic region. Following these observations, different studies generalized the methodology to teleseismic phase detections (Lin et al. 2013; Nishida 2013b; Boué et al., 2013, 2014a). In these studies, it was shown that the retrieved seismic phases carried structural information of the deep Earth. Poli et al. (2015) successfully identified a D’’ reflector at 2530 km depth by extracting a PdP phase on correlation functions computed between Finland and Japan. While the signal-to-noise ratio (SNR) recovered from correlation of ambient noise is lower, it benefits from array processing on both the source and receiver side (as will be developed here). This new class of signals has been proven to be a good supplement to earthquake data for deeper and deeper targets. In parallel, studies have extracted body wave sources of seismic signals linked to oceanic activity in the ambient noise (Gerstoft et al. 2008; Landès et al. 2010; Zhang et al. 2010; Retailleau & Gualtieri 2019; Li et al. 2020). These studies lead the way to a better understanding of the body waves observed by correlation methods and thus to their better use for imaging. Other studies are developing processes to invert ambient noise correlations for the structure and the sources simultaneously (Fichtner et al. 2016; Sager et al. 2018). Alongside the development of body wave retrievals from ocean and solid Earth interactions (Boué et al. 2013; Nishida 2013), different studies point out the possibility to use longer period body waves retrieved from large earthquake late coda arrivals (e.g. Boué et al. 2014a; Huang et al. 2015; Phạm et al. 2018; Wu et al. 2018). This particular method usually shows a richer wavefield with very good SNR as compared to a noise-based approach, but one should be careful because the correlation does not necessarily converge toward the natural Earth impulse response (Poli et al. 2017; Kennett & Pham 2018). As compared with noise-based correlations derived from a ballistic wavefield emanating from Ocean waves, signals retrieved from coda waves between station pairs are less straightforward to exploit to image 3-D heterogeneities of the Earth. Nonetheless, coda wave correlations have proven to be very efficient to extract the weak signature of velocity deviation from 1-D reference models for the deepest part of the Earth (Tkalčić and Phạm 2018). Wang & Tkalčić (2020) tend to demonstrate, with very promising synthetic examples, that this type of signals could also be used to resolve 3-D structures of the inner core.

In this paper, we produce an original 1250 km × 950 km aperture (at the CMB 680 km × 520 km) image of the deep mantle beneath the North Atlantic Ocean using P and PcP body waves extracted from correlations of the secondary microseim computed from about 900 stations in the United States (US) and Europe (Fig. 1). This area is located North of the region studied by Weber & Körnig (1992) and East of the previously mentioned studies by Yao et al. (2015) and Durand et al. (2019). These two study areas are represented in Fig. 1. We first show that P and PcP phases are observed on both the causal and anticausal part of the correlation functions in the secondary microseism period band (3–8 s) after stacking a 1-yr data set. Following this, we use this very large correlation data set (about 100 000 paths, i.e. turning points at depth) to constrain regional reflectivity variations in the few hundreds of kilometres around the CMB using a double-array analysis.

2 RETRIEVING P-WAVE ARRIVALS FROM NOISE CORRELATIONS

Our data set consists of 1 yr (2014) of continuous seismic signals recorded by the vertical component of roughly 900 broadband stations in Europe and the United States (Fig. 1). The daily time series were downloaded from the Incorporated Research Institutions for Seismology (IRIS, http://www.iris.edu/mda) and European Integrated Data Archive (EIDA, http://www.orfeus-eu.org/eida/eida.html) data services. Details about the downloading and networks can be found in the Supplementary Materials S1 and Retailleau et al. (2017).

The correlation functions are computed between each station in the United States and each station in Europe, provided the stations were operating simultaneously for at least 100 d. The correlations are evaluated on 4-hr segments, frequency normalized, for all segments that do not contain strong earthquake signals. All correlations for a given station pair are then stacked over the entire year to optimize the convergence of the correlation function toward the Earth impulse response. The resulting correlation functions are filtered between 3 and 8 s of period, that is in the secondary microseism period band. Figs 2(a) and (b) show all these correlations as a function of interstation distance and lag-time (the number of correlations stacked at all distances can be found in Fig. S2). The two panels correspond to waves travelling from the United States to Europe (Fig. 2a, anticausal part of the correlations) and the opposite direction (Fig. 2b, causal part). Teleseismic P and PcP arrivals are clearly visible for both directions and agree very well with the arrivals predicted using the TauP toolkit (Crotwell et al. 1999) using the IASP91 model (Kennett & Engdahl 1991). No significant other arrivals are visible on Figs 2(a) and (b) leading to two options for a potential D” discontinuity: (1) No such discontinuity exists in this area or (2) the 3-D structure is sufficiently complex to stack destructively in these averaged sections. Previous studies support the second option.

A significant point is the apparent symmetry of the correlations, at least on average, observed for this P-wave propagation. We represent both the causal and the anticausal parts of the correlations in Fig. 2(c) for a direct observation. While a good symmetry of the correlations is not necessary, it is a meaningful argument to support the sufficient convergence of the correlation towards the Green function in this frequency band. In other words, a similar bias that would be observed in both directions is very unlikely. Fig. 3 shows the PcP source regions (Fresnel or stationary zones, e.g. Snieder 2004) that match the stationary phase condition for the interference between PcPPcP and PcP at 7 s of period. The two Fresnel zones correspond to the two propagation directions between a station in the United States and a station in Europe. These two source areas are clearly distinct and it can therefore be said that anticausal and causal signals do not share the same microseism origin. We can also note that both areas are large compared to our study area. This implies that different locations in our study area are illuminated by very similar source areas (for causal and anticausal signals, respectively), making it possible to confidently interpret relative variations of amplitude in our results. This phase combination (PcPPcP and PcP) is one of the possibilities to extract the PcP phase. For instance, a combination between a PKPPcP and a PKP phase is another possibility. Isolated energetic noise sources do not seem to play an important role over the entire year; no spurious phase is visible on this average section near the P-wave arrival (Li et al. 2020). We speculate that noise sources are sufficiently well distributed during the year and that scattering is strong enough at this period band to randomize the wavefield leading to relatively unbiased correlation functions. While this result is very promising we still choose to be careful with our interpretations and prioritize robustness.

3 MAPPING THE REFLECTIVITY OF THE LOWERMOST MANTLE

Figs 2(a) and (b) show spatially averaged signals meaning that sharp lateral variations naturally cancel out over the whole region. However, our data set is sufficiently large to extract regional variations of reflectivity in the area highlighted by yellow dots in Fig. 1. Our objective is to map the contribution of arrivals precursory to the PcP and thus potential reflection above the CMB. We describe here the main steps of the method; more details can be found in Supplementary Material S1. To verify that our workflow does not include any bias and correct from the amplitude variations with geometrical and propagation effects, we apply the same process to synthetics computed in the 1-D ak135 model with attenuation (ak135 of Kennett et al. 1995; synthetics from IRIS DMC 2015; van Driel et al. 2015; Krischer et al. 2017). We generate the signals from a vertical force source located at each station location in Europe and recorded at all the U.S. stations. We thus reproduce the propagation estimated by the correlation process applied on the real data. Because our correlations signals are fairly symmetric, we stack them to produce a single image, but the following workflow is also applied on both the causal and anticausal parts of the correlations to get independent images. As compared to earthquake-based signals, a noise correlation technique that turns any seismic station into a virtual source allows the use of array processing on both the (virtual-) source and receiver side. This is particularly important to overcome low SNR of single station-to-station detections. At that scale, and with similar targets, this so-called double-beam analysis was first applied using nuclear explosion test sites as the source array (Krüger et al. 1993; Scherbaum et al. 1997), and then it proved its usefulness in the noise-imaging framework (Boué et al. 2014b).

Our workflow, applied on the data and on the synthetics, can be decomposed as follows (a detailed diagram is available in Fig. S1):

1. Constitution of subarrays in both the United States and Europe by selecting, on each point of a grid of spacing 110 km × 140 km, all stations within a 150 km radius. Subarrays are then paired with all subarrays from the other continent. The corresponding reflection points below the Atlantic Ocean are represented with yellow dots in Fig. 1.

2. Computation of double-beam transforms (time-slowness analysis, Rost & Thomas 2002; Boué et al. 2014b) using each combination of sub-arrays provided they contain at least 200 pairs of stations (i.e. correlation functions). Fig. 4(a) shows a vespagram (time-slowness representation) obtained for a single combination of subarrays. For sake of simplicity, we here consider that the slowness is the same at both subarrays for a given phase.

3. Construction of time and slowness values pairs of potential reflection points below the lowest point of the P-direct wave are computed using the TauPtoolkit (Crotwell et al. 1999) for a propagation in a 1-D model ak135 (Kennett & Engdahl 1991). These hypothetical reflection points, each corresponding to a given depth, are represented on top of the vespagram in Fig. 4(a) as green dots. The amplitude of the beam is extracted at each depth as shown in Fig. 4(b). We name R_z this depth dependent reflectivity function derived from our beamforming analysis. This figure shows the PcP phase at the correct depth (∼2890 km, corresponding to the CMB) along with a precursory arrival labelled as P*P. We name I_z the function of depth z obtained through the same process using synthetic waveforms computed in ak135. The generalization of this analysis to all subarray combinations (each of which consists of at least 200 stations) is presented in Fig. 4(c). Each R_z reflectivity function, extracted from each array pair, is represented as a function of the interarray distance along the great circle. The first signal, at the beginning of each curve, corresponds to the P arrival. The associated depth varies with the turning point of the phase, thus with the interarray distance.

The second signal, at about 2890 km depth, corresponds to the detection of the CMB. This signal is obtained with most subarray combinations. This first-order observation confirms that we can extract large-scale structural information even using a subset of correlation signals. Fig. 4(c) also reveals that this signal is not depicting a laterally homogeneous interface; indeed very little signal is observable at the CMB at some distances, and some precursory reflected phases appear at about 2600 km. Absolute amplitudes are difficult to assess using correlation functions and would require a better understanding of the noise sources to be unambiguously interpreted. Nevertheless, we note that for a given pair of stations the same noise sources are illuminating both the CMB and potential shallower reflectors. It is therefore justified to discuss the relative amplitude of the reflected phases with respect to the reflection at the CMB. As mentioned earlier in the paper (Fig. 3), one can also note that the size of our target is smaller than the noise source Fresnel zone of the PcP signals by a factor 4. This suggests that the local comparison between the subarrays is also justified.

To build the reflection maps, we first normalize each reflectivity function R_z by the average amplitude <R₂₈₉₀> at the CMB (2890 km of depth) with <.> denoting the average over all paths, and we proceed identically with the synthetics analysis by dividing each function I_z by <I₂₈₉₀>. We then divide each function R_z and I_z by the amplitude I₂₈₉₀ of the corresponding path at the CMB to correct the amplitude variations due to geometrical spreading and attenuation. The final R_z and I_z functions can be seen as a proxy of the P-wave reflectivity relative to the CMB discontinuity. After this normalization, the theoretical CMB within I_z functions should show up with a value of 1.

Most studies that look for lower mantle heterogeneities use shear waves propagating between about 70° and 85° of epicentral distances. In this distance range, precursory arrivals are usually well separated from S and ScS branches (about 10 s separation, e.g. Kito et al. 2007; Lay & Garnero 2011). In our case, using the P waves emerging alone from correlations (no S waves) we work between 51° and 73°. For the largest distance, we expect to have less than a dominant period lag between PcP and potential precursory arrival (Fig. S3). To verify the geographical distribution of the heterogeneities as seen in Fig. 4c (as a function of distance) functions R_z (and similarly I_z) are projected on a latitude-longitude 950 km × 1250 km (680 km × 520 km at the CMB) grid to form a 3-D final image. We sample the grid by a longitude and latitude spacing of 0.5 degrees. At each point of the grid we average the functions R_z reflected at the location provided there are at least four R_z functions (i.e. array combinations). Fig. 5(b) shows the interpolated values R_z at six different depths from 2470 km in the mantle to 3170 km in the outer core. The same process is performed from the synthetic functions I_z whose results are shown in Fig. 5(a) as a reference. This reference figure shows that our processing does not introduce lateral bias and is calibrated to resolve the CMB. We perform a bootstrap analysis to assess the quality of the 3-D image we reconstruct (Fig. S4). We perform 1000 random samplings of 75 per cent of the subarray combinations. Figs S4(a) and (b) represent the average value and standard deviation obtained between all tests at the different locations can. The averaged 3-D image (Fig. S4a) is very consistent with the overall result in Fig. 5(b). The low value of the standard deviation, that is to say various combinations of arrays give similar estimations for the same reflector point, is a good indication that our results are robust; small scale heterogeneities of the mantle are not dominating the image.

4 DISCUSSION

As expected, the result obtained from synthetics generated using a 1-D model, in Fig. 5(a), shows a laterally homogeneous reflection at the CMB (with a value of 1 due to our normalization) and no reflection above or beneath. This shows that our array processing and overall imaging workflow does not introduce artefacts. The result obtained with noise-based data, in Fig. 5(b), also shows a major reflection at the CMB but with strong lateral variations and reflected energy at smaller depths. This coincides with the observations made in Fig. 4(c). At depth of the CMB, R₂₈₉₀ shows a clear high reflectivity anomaly (blue area) of about 200 km × 100 km in the northwest region. Since the amplitudes are normalized so that the average amplitude at the CMB is set to one, this observed high amplitude anomaly implies a strong reflectivity increase relative to theoretical CMB. This structure extends upward in the lower mantle to 2600 km depth. We can also note that, as was seen in Fig. 4(c), there is no energy reflected in the outer core, except from the amplitude introduced by finite frequency effect of the PcP pulse width that is centred at the CMB depth (this effect already vanished on the first slice below the CMB at 3030 km depth). One could question if the precursory signals could be an effect of broadening of the PcP pulse. We expect a broadening effect of stacking to affect depths above and below the CMB after migration.We note in Fig. 4(c) that the wiggles detected prior to the PcP have a much larger amplitude than the small amplitude oscillations visible below the CMB. The comparison with ak135 synthetics confirms that clear reflected waves originate from above the CMB.

The process we use to detect arrivals around the CMB combines the envelope vespagrams generated at the multiple reflection points. While this process makes the observations more robust, the stack smooths lateral variations. As mentioned earlier, we do not expect to recover absolute amplitudes with ambient noise correlations and we aim at detecting the precursory arrivals and to compare their relative amplitude with PcP. The CMB amplitude is also normalized by its average value, meaning that a value of 1 corresponds to the expected CMB contrast (defined by ak135). The observed structure, which extends upward to 2600 km, is too thick for a ULVZ, however, it is consistent in terms of size with other D’’ observations in other regions (Lay 2015; Li et al. 2019). In contrast, the southeast region of the image shows a lower reflectivity at the CMB. The transition between these two areas reveals complex patterns with a typical size of about 50 km. Fluctuations of the amplitude at the CMB range from 0.5 to 1.5. We will interpret these variations of about 50 per cent of the reflectivity cautiously due to the difficulty to properly calibrate correlation amplitudes (as explained above) and also due to the possible impact of a critical reflection as discussed later. However, we show a raw image without further processing to highlight the spatial variability of the structures. Indeed, the characteristic size of the heterogeneities we see on Fig. 5(b) are in good agreement with previous observations in neighboring areas of roughly 300 km thick and laterally discontinuous D” reflectors (Yao et al. 2015; Durand et al. 2019). A full comparison would probably require combining earthquake and microseism data as well as the S waves extracted from microseism excitation correlations. Such microseism based S wave signals have been proved to exist for very energetic oceanic events (Nishida & Takagi 2016), but appear to be more difficult to extract than P waves

To verify that our results are not significantly impacted by noise source distribution or transient signals, we compare the images produced by the anticausal and causal part of the correlation separately (i.e. waves travelling from the United States to Europe and in the reverse direction). The results are presented in Fig. S5. The large-scale variations as discussed above are visible with a positive anomaly in the northwest area in both cases. These images also help to estimate a lateral resolution at the CMB for Fig. 5 of roughly 100 km. Indeed, at lower scales there are some differences between the results obtained with the causal and the anticausal part of the correlations. We represent the difference between the results obtained with the causal and anticausal parts of the correlations in Fig. S6. Other parameters such as the intersubarray distance or the number of stacked signals per bin to construct the image are also possible sources of artefacts and biases in our final result. To check that there is no such bias from these parameters we compared them to the image at the CMB (see Fig. S7). No significant bias emerges from these comparisons. For instance, isovalues of the mean inter-subarray distance are crossing the main features of the image. Finally, the bootstrap analysis (Fig. S4) shows that no subarray pair dominates the results and creates unreliable observations.

Until now, we only used the envelope of the vespagrams for our analyses. We also perform the same process while keeping track of the phase information, meaning not computing the envelope, to analyse the relative polarity of the precursory arrivals. Fig. 6(b) shows these observations of R_z/R₂₉₈₀ averaged over two different areas (zone A and zone B) and for both the data and the synthetic waveforms. As already discussed, the limited bandwidth of the secondary microseism limits the temporal resolution of the pulses we extract resulting in a possible overlap of the CMB arrival with the precursory signal. It is thus not straightforward to conclude on the phase of precursors. There does not seem to be a polarity flip, but we cannot conclude about a potential phase shift to better constrain the real reflection coefficient. Looking at the averaged pulse observed for the positive anomaly in the northwest area (zone A, red rectangle in Fig. 6) it seems that precursory arrival amplitude is high when the CMB reflection is also high. In terms of reflection coefficients, this coincides with a high-velocity D” layer for both P and S waves (as it will be developed later). Looking now at the pulse observed for the other region, which shows a lower reflectivity at the CMB (zone B, black rectangle, Fig. 6a), the existence of a precursory signal is less evident. It is interesting to note that the synthetic signals, for which the same processing workflow was applied, show in addition to a clear CMB reflection a small precursory signal at 2750 km depth, where ak135 experiences a velocity gradient decrease. Our simple migration scheme allows us to properly capture this anomaly. Finally, both data-based averaged pulses (Fig. 6b) show an asymmetric pattern with apparent lower frequencies above the CMB compared to below it. This is mostly due to a stretching effect introduced by the depth migration of the correlation functions (see green dots on Fig. 4a). This effect appears on the synthetic pulses shown on Fig. 6(b), but it is particularly clear if we take a single signal (Fig. 7g).

To go deeper in our interpretation of the precursory signals we are observing, we now consider very simple 1-D models of positive P-wave velocity increase in the last kilometres of the lower mantle. Starting from ak135 as a reference (in black, Fig. 7), we consider four different perturbed models to describe a hypothetical D” layer (Fig. 7a): two models with a simple positive step discontinuity 250 km above the CMB, and two models with a smoother gradient starting around 300 km above the CMB. In both cases a P-wave velocity perturbation (δVp), is fixed to 2 and 4 per cent, respectively (δVs = 0), relative to ak135 within the entire D” layer. We then recreate the configuration of one set of subarrays, with a source on one side of the ocean and 249 stations on the other side (assimilating the virtual station on one side of the ocean and all the others on the other side). The average distance is 62 degrees. Synthetic signals are generated using axisem (Nissen-Meyer et al. 2014) for a vertical point force source at the surface. The signals are filtered in the 3–8 s period band, similarly to the data. We compute the vespagram for each model following the process used to generate Fig. 4(a), without taking the envelope (Fig. 7b–f). We extract each function I_z from the vespagrams, using the migration described in Fig. 4. We represent the different I_z functions in Fig. 7(g) (the colours follow Fig. 7a). Fig. 7(g) is comparable to Fig. 6(b) but for a single reflection point (no spatial averaging). Indeed, in this new and simpler case, the stretching effect due to the migration is clearly visible.

Both step models generate a clear precursory signal centered around 2650 km and visible on Fig. 7(g). We can observe a tenuous PdP arrival in Fig. 7(d), but similar observation is difficult for a 2 per cent step increase (Fig. 7c). The signal of the 4 per cent step increase visible on Fig. 7(g) is the closest of our observation for zone A on Fig. 6(b) (red rectangle). Both models with a smoother D” velocity increase (grad 2 per cent and grad 4 per cent) show a slight stretching towards lower depth (expected for a velocity increase) with no clear precursors. These simple examples show that, using our workflow, we can detect a 250-km-thick D” and distinguish between sharp (step) and smooth (gradient) velocity increase. By comparing Figs 7(g) with 6(b), the strong reflectivity anomaly on the northwest edge of our study area (visible on Fig. 5 and discussed as zone A on Fig. 6) would probably give a best fit with an intermediate model between step 4 per cent and grad 4 per cent. The absence of a clear precursory arrival on many of the signals from other regions of our study zone (Fig. 4c) could indicate smoother velocity variation.

Based on theoretical reflection coefficients on both the CMB and on a hypothetical D” discontinuity, we can estimate the minimal requirements in terms of structural variation to explain our observations. First, theoretical reflection coefficients are computed for the CMB and a 250 km thick (step) D” for all combinations of parameters that define the anomaly, including δVp, δVs and δρ (density). These parameters are ranging from −20 to +20 per cent. In order to match our signal processing workflow, all reflection coefficients are then normalized by a not-perturbed reference CMB reflectivity at each epicentral distance. Fig. S8 gives an example of representation for these 3-D parametric normalized reflection coefficients. Fig. 8 shows reflection coefficients as a function of incidence angles and epicentral distances (without normalization) for 5 models including ak135. By comparing normalized reflection coefficients to what we observe Fig. 5(b) and for the high reflectivity anomaly (zone A, Fig. 6a) with a high reflectivity above (R₂₆₁₀ = 1) and at the CMB (R₂₈₉₀ = 1.2) the smallest perturbation needed is about δVs = +8 per cent and δVp = +4 per cent if we consider a 250-km-thick D” step perturbation of ak135 (similar to the Fig. 7, and here named model step1). Using the same approach we can evaluate that a decrease of Vs of about δVs = −4 per cent within the last 250 km of the lower mantle (model step2) could explain zone B (Fig. 6a) where we observe a low D” reflectivity (R₂₆₁₀ = 0.2) and a low CMB reflectivity (R₂₈₉₀ = 0.8). Fig. 8 shows theoretical reflection coefficients for both step 1 and step 2. These two simple models correspond to a first-order evaluation that tends to show that Fig. 5(b) contains structural information. However, we note that density is poorly constrained by this approach and can stay almost the same with less than a per cent of variation. Vs perturbation is also not very well constrained especially for a positive Vp anomaly, which leads to a critical reflection at short epicentral distances (see Figs 8 and S8). A higher sensitivity on Vp is not surprising since we are only dealing with P waves.

The velocity perturbations estimated here are a little high and difficult to explain in the framework of a perovskite/post-perovskite transition where we expect weak velocity perturbation for D” (Cobden et al. 2015). But two remarks can be made to support our observations. First, a P-wave velocity increase of about δVp = +4 per cent is in good agreement, although in the upper range, with previous seismological studies in neighbouring regions (Weber & Körnig 1992; Durand et al. 2019). A greater Vs discontinuity (δVs = +8 per cent), on the other hand, farther away from previous observations (Yao et al. 2015; Durand et al. 2019), but again not well resolved in our case. For a fixed value of δVp = +4 per cent, a positive jump of Vs has a relatively small influence on the observed P-wave reflectivity of the D” layer; the observed δVs is in fact mostly controlled by the slight increase of the CMB reflectivity (see the difference between model step1 and step4 in Fig. 8), which is less resolved due to PdP pulse broadening and possible overlapping with PcP. Secondly, the fact that we are observing what appears to be relatively strong velocity perturbations could be related to the presence of critical reflection at relatively short epicentral distances. The different models with a positive δVp presented in Fig. 8 (step1, step3 and step4) show sharp increases of reflectivity between 50° and 70° of epicentral distances (similar range that we are using in this study). By possibly mixing pre- and post-critical reflections during the required spatial averaging of the beamformed correlation functions (Fig. 4c), we are aware that we may introduce non-physical reflectivity values in our result (Fig. 5b). Even if the spatial averaging also helps us to smooth these possible biases, the quantification of the influence of such critical reflections would require an a priori knowledge of the D” structure and could be a subject of a separate study (with synthetic data). Finally, we speculate that the amount of data involved in our method might also be a very efficient way to precisely measure the critical angle for precursory reflections and thus give a very strong constraint on δVp and thickness of the D”, but this is beyond the scope of the present article, and will be addressed in future work.

5 CONCLUSIONS

We have built a 3-D image of the CMB region under the Atlantic Ocean from the correlations computed between 1 yr records of seismic stations in Europe and in the United States. Fig. 5(b) constitutes a first tentative, 3-D image of the CMB region derived from microseism excitation. Even if critical reflection may introduce a bias of amplitude in the image, we are confident in the shape and spatial extension of our detection. We compared our observations to different possible models containing a hypothetical D” layer. A model containing a Vp velocity increase of 4 per cent is consistent with our observations and neighboring studies.

Earthquake data and careful synthetic tests would help to better constrain the structural properties. Our methodology would also greatly benefit from other observations such as secondary microseism S-wave reconstruction (Nishida & Takagi 2016) and/or high-frequency signals from specific storms (for example typhoon Ioke, Retailleau & Gualtieri 2019).

The lowermost mantle is known to incorporate multiscale heterogeneities with direct consequences for both mantle and core dynamics (Tkalčić et al. 2015). Detailed studies of high frequency scattered waves and lower frequency waves used for global tomography already give a good picture of these heterogeneities where earthquake data are available. We believe that the present approach can help link up these different scales to form a more comprehensive image of the lowermost mantle by providing complementary information where earthquake data are scarce.

SUPPORTING INFORMATION

supplementary_review.pdf

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

ACKNOWLEDGEMENTS

The data set was downloaded and processed using python and the seismological community scientific library obspy (Krischer et al. 2015). The figures were produced using Python and MATLAB and the Python module basemap was used to produce the map. The facilities of IRIS Data Services, and specifically the IRIS Data Management Center, were used for access to waveforms, related metadata, and/or derived products used in this study. IRIS Data Services are funded through the Seismological Facilities for the Advancement of Geoscience (SAGE) Award of the National Science Foundation under Cooperative Support Agreement EAR-1851048.

The computation of the correlation functions presented in this paper was performed performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is partly supported by the Equip@Meso project (reference ANR-10-EQPX-29–01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche. This work has been supported by a grant from Labex OSUG@2020 (Investissements d'avenir—ANR10 LABX56) and Fondation Simone et Cino Del Duca, Institut de France (Prix scientifique 2013). We acknowledge the support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 742335, F-IMAGE). The research of LR was supported by Pacific Gas and Electric. We thank K. Schaukowitch and the Hume centre for writing and speaking at Stanford for useful comments that helped us to improve the manuscript. We thank the Editor, Michael Ritzwoller, Stéphanie Durand and two anonymous reviewers for their comments and suggestions which helped to improve the quality of this paper.

REFERENCES

APPENDIX

Tables A1 and A2.