Analysis of core-mantle boundary seismic waves using full\-waveform modelling and adjoint methods

Using spectral-element and adjoint methods, we investigate body waves interacting with the Earth's most dramatic interface, the core-mantle boundary (CMB). Intermediate-to-high frequency seismograms are computed incorporating topography models. We analyse the sensitivity of many seismic phases interacting with the interface. The study aims at showing effects of CMB structure on synthetics and highlights difficulties of imaging this region due to strong trade-off between velocity variations and topography. Synthetic waveforms computed at dominant periods of 6-18 seconds are used in order to observe time shifts due to topography and calculate the adjoint sensitivity kernel. We focus on diffracted, core reflected and refracted P and S waves. The sensitivity kernel depicts first Fresnel zones of these phases and others that may contribute to narrow time window, although unpredicted by ray theory. We perform comparisons between time shifts due to topography models made on full-waveform synthetics to ray theoretical predictions to assess methods usually deployed for imaging CMB. This shows that for most phases ray theory performs well enough with some accuracy loss. We propose that using relevant seismic phases simultaneously in full-waveform inversion may improve CMB topography imaging. However, it seems necessary to jointly invert for velocity variations due to D\"layer, which is so far poorly understood and presents a challenge when it comes to identifying its effects on traveltime data. In our further research, a FWI workflow is being developed and aims at addressing these issues.


Introduction
The core-mantle boundary (CMB) is a major and sharp discontinuity inside the Earth, separating the mantle, composed of solid rocks, from the outer core, made of liquid iron alloy. Specifically, it lies between the mantle rocky aggregate and the liquid iron alloy outer core. Key phenomena affecting the Earth evolution and dynamics occur on each side of this boundary. For instance, on the mantle side, the CMB region is believed to be the source of mantle plumes, which An important hindering factor is the complexity of the seismic structure in the bottom 300-400 km of the mantle. This includes the D" discontinuity, which is not present everywhere and is thought to be related to the phase transition from bridgmanite to post-perovskite, [e.g. Hirose et al., 2017, Cobden et al., 2015, ultra-low seismic velocity zones (ULVZ), which are thin pockets where compressional and shear wave velocities (V P and V S ) drop by up to 10 % and 30 %, respectively, [e.g. Yu and Garnero, 2018], strong, small-scale anomalies induced by cold oceanic slabs that have sunk down to the CMB, [e.g. Hung et al., 2005, Whittaker et al., 2015, Borgeaud et al., 2017, and the large low shear-wave velocity provinces (LLSVPs), defined as regions a few thousands of kilometers across where shear velocity drops by a few percent, and which are located beneath Africa and the Pacific ocean, [e.g. Li and Romanowicz, 1996, Ritsema et al., 2011, Lekic et al., 2012, Garnero et al., 2016, McNamara, 2018. Because they may be closely related to mantle dynamics, LLSVPs may strongly influence CMB topography. Numerical simulations showed that depending on their exact nature, clusters of purely thermal plumes or thermo-chemical piles, LLSVPs would trigger either uplifts or depressions in the CMB [Lassak et al., 2010, Deschamps et al., 2018.
Ideally, one would like to resolve CMB topography structure, while simultaneously obtaining an updated velocity model for the mantle. It has been shown by many researchers that CMB topography and mantle velocity variations exhibit strong trade-off, which makes the inference of a topographic map along the core-mantle boundary a hard task, [e.g. Garcia andSouriau, 2000, Sze andvan der Hilst, 2003].
The nature of difficulties in resolving lowermost mantle structure can be analysed and addressed by investigating seismic waves which are sensitive to lowermost mantle structure and have been used to constrain the topography of the CMB. These seismic phases are characterised by complex wave phenomena occurring while travelling from the source to receivers, which makes it cumbersome to determine which type of seismic waves has a maximum CMB sensitivity. Nonetheless, there have been many studies analysing the compressional and shear wavespeed structure of the lowermost mantle using several related body wave phases as well as differential traveltimes between these, [e.g. Li and Romanowicz, 1996, Garnero, 2000, Garnero et al., 2016, Ritsema et al., 2011, Tanaka, 2002, Lekic et al., 2012, Muir and Tkalčić, 2020. A summary of the progress made during the last decades and relevant literature can be found found in [Gurnis et al., 1998, Garnero, 2000 and [Koelemeijer, 2020]. In a different approach, studies [Colombi et al., 2012, based on axisymmetric spectral-element modelling , have implemented boundary and volume sensitivity kernels within a spectral-element code given spherically symmetric background models. They used this implementation to address the sensitivity of PcP and P dif f and showed a suite of sensitivity kernels to CMB structure [Colombi et al., 2012]. Their research gave insights into how to incorporate sensitivity of traveltimes of P dif f , PKP, PcP and ScS for imaging the CMB  in a joint inversion for boundary and lowermost mantle velocity structure in order to tackle the trade-off between parameters. Our study is stepping further by using the full-waveform synthetics calculated with spectral-elements in a cubed sphere, as implemented in SPECFEM3D_GLOBE, where the global wave propagation effects can be more accurately accounted for and implemented into the full-waveform processes in a consistent way.
The goal of this work is to investigate the finite-frequency sensitivity of frequently used seismic phases for lowermost and core-mantle boundary structure and provide a kernel gallery showing their sensitivity to boundary and volumetric isotropic structure. This is following the formulation of finite-frequency traveltime sensitivity given by [Marquering et al., 1999, Dahlen, 2005 and [Tromp et al., 2005] with the use of adjoint methods. The main contribution of this work is to better understand how to optimally use traveltime data of particular Pand S-waves for imaging the CMB and show how their sensitivities can be incorporated in a full-waveform inversion scheme that is targeted for resolving the lowermost mantle and CMB topography. We investigate a suite of Pand S-wave phases and their traveltimes as observed in synthetic seismograms computed for intermediate frequencies (0.01-0.08 Hz), which are at the higher end of frequencies used in global waveform inversion studies. We focus on the seismic phases PcP, PKP, PKKP, PcS, ScS, SKS and SKKS. We deploy a full-waveform, finite-frequency approach for analysing CMB related seismic waves, using well-established techniques based on spectral-element waveform modelling and adjoint methods for calculating the Fréchet sensitivity kernels, [e.g. Komatitsch and Tromp, 1999, 2002a,b, Tromp et al., 2005.
To address the effects of adding topography and/or 3-D velocity variations on the traveltimes of these phases, we compare existing ray-based methods of inferring core-mantle boundary topography to cross-correlation time shift measurements made on the spectral-element synthetics. For this, we use full-waveform synthetic seismograms for events selected especially for CMB topography studies. The analysis is made by comparing traveltime anomalies predicted using ray theory, for the CMB topography models mentioned above, to residual traveltimes measured on full-waveform synthetics computed for the same CMB topography models. The traveltime comparisons are done similar to [Bai et al., 2012] and [Koroni and Trampert, 2016]. This comparison could provide insights into the usability of the traveltimes of our selected phases within the ray theoretical framework. Using full-waveform modelling as realistic recordings, we can address the effects of varying the model on specific observables. This analysis should benefit from investigating the separate effects of 3-D mantle velocity structure and boundary topography on the traveltimes, when only these model parameters are perturbed during the simulations. The separate effects are examined by computing synthetics in a 3-D model and in a 1-D velocity model plus core-mantle topography. Our results show that for 1-D background the ray-theoretical framework works better for S phases than for P phases, while it deviates from a good prediction when the background velocity model is 3-D. We also find that the boundary sensitivity of almost all phases is much less pronounced than their clearer sensitivity to volumetric parameters, both compressional and shear velocity. This confirms the reported trade-off between CMB and 3-D velocity structure, and also indicates that all types of kernels shown here could be used to provide better constraints on relevant structures, as the CMB topography cannot solely be responsible for the traveltime shifts and the velocity structure imprint is significant.
The paper is structured as follows: We first present the methods for carrying out the traveltime analyses and comparisons. Then we explain the specifics of Fréchet sensitivity analysis. We present our results and we discuss them by addressing the effects of different CMB topography models, the influence of the background velocity model (1-D or 3-D). We also consider the different effects on Pand Sseismic phases. Finally, we explain the resulting Fréchet kernels in the context of full-waveform inversion and trade-off between the two parameters.

Computation of synthetic seismograms
In order to enhance our knowledge about core-mantle boundary structure by investigating its effects on Pand Sseismic phases, we perform full-waveform modelling and focus on frequencies ranging from 0.01-0.08 Hz, given the computational limitations, still lying within an acceptable observational resolution. The theoretical raypaths of the chosen seismic phases are shown in Figure 1, which is produced using 1-D ray tracing with tauP [Crotwell et al., 1999] and using model PREM [Dziewoński and Anderson, 1981] as reference. The computation of synthetic waveforms is done using the community code SPECFEM3D_GLOBE [Komatitsch and Tromp, 2002a,b] in 1-D and 3-D mantle velocity models, namely PREM [Dziewoński and Anderson, 1981] and S20RTS [Ritsema and van Heijst, 2002], respectively. The synthetic waveforms are used to imitate realistic waveforms for the comparisons between ray theory and measurements and to assess the direct effects of CMB topography and 3-D velocity variations on the traveltime of the seismic phases.
For the comparisons of traveltimes derived using ray theoretical methods to measurements done on full-waveform synthetics, existing core-mantle boundary topography models are added, resulting in PREM or S20RTS plus CMB topography. A model with very low topographic variations is used, derived by [Tanaka, 2010] referred here as TK. In order to examine the effects of a variety of topographies with larger perturbations, we compute synthetic waveforms by also adding the models by [Xiang-Dong Li et al., 1991], denoted here as LM91, with peak-to-peak variation scaled up to 8 km, and two topography models derived from geodynamic simulations T1-pPv and TC1-pPv from [Deschamps et al., 2018], denoted here as T1 and TC1 for simplicity. In order to make the geodynamics models more applicable to  Table 1: Main features of the earthquake events used in this study. These are recorded at the station locations of the Global Seismographic Network.
our source-receiver geometry, model TC1 is rotated so that its spherical harmonic degree-2 pattern match best that of S20RTS S-velocity variations right above the CMB (this is done by performing a simple grid search). The same is done for model T1, but with the sign of S-velocity variations inverted. This is because S-velocity anomalies at the CMB, and CMB topography variations are expected to be positively correlated for thermo-chemical models (TC1), but negatively correlated for purely thermal models (T1) [Deschamps et al., 2018]. These topography models are derived from purely thermal (T1) and thermochemical (TC1) simulations of mantle convection modelling. They were filtered for spherical harmonic degrees = 0 to = 20, and have absolute peak-to-peak amplitude around 20 km, for both models. Most of this amplitude is however accommodated by deep depressions caused by downwellings reaching the CMB, while positive topography associated with plume clusters in T1 and depressions associated with thermochemical piles in TC1 have height or depth, respectively, of about 2 km. Having various possible topography models allows us to investigate the reliability of ray theoretical methods in a more conclusive manner. All CMB topography models used here along with the dataset constructed by the earthquake events shown in table 1 are displayed in Figure 2.
For the simulations, we set the resolution such that the resulting noise-free waveforms have dominant period of about T = 11 s. Some modelling parameters which are present in real data can be important, such as attenuation, and can have complex effects on traveltimes of seismic phases, thus they are not considered for our simulations, since we want to focus on the effects of CMB topography and 3-D lateral isotropic variations without further complicating the waveforms.
To obtain a realistic dataset, we use seven events with a range of focal depths, source time function duration and moment magnitudes in the range M w = 5.5 − 7.5. The events are summarised in table 1. The receiver locations are set to those of the Global Seismographic Network and we make measurements for each source-receiver pair as well as ray theoretical predictions for the same pairs. This is done in time windows around the seismic phases of interest. The synthetic data are processed by filtering using a band-pass between the frequencies 0.01-0.08 Hz. The synthetic waveforms are used for two different purposes: Firstly, for comparisons of traveltime shifts caused by topographic variations in 1-D or 3-D mantle structure models and, secondly, for computing the sensitivity kernels for selected time windows encapsulating the maximum peak of the energy of our selected seismic phases.

Comparison of ray theoretical time shift prediction and cross-correlation measurements
In this section, we describe the prediction of traveltime anomalies due to perturbations of CMB topography using ray theory (RT), and as measured on full-waveform synthetics computed using the spectral element method (denoted by FW, which stands for full-waveform) using cross-correlation. The motivation of this analysis is to check whether we can reliably predict a topographic perturbation given the ray-theoretical framework. If this is true, the measured by cross-correlation time shift and the predicted time anomaly for a specific topographic variation should be identical. In the opposite case, where the two values diverge, there are two implications: firstly, the ray-theoretical framework is not sufficient for translating the time shift measured to topographic variation; secondly, there are velocity effects which are not adequately taken into account when using a linearised inversion. Traveltime perturbations due to the CMB topography are computed in the framework of ray theory following [Morelli andDziewonski, 1987, Tanaka, 2010], using the expressions given in eq. 1 when the ray reflects off the CMB (e.g., ScS), and eq. 2 when the ray transmits through the CMB (e.g., SKS): Phase Before selection Selected (PREM) Selected (S20RTS)  ScS  1273  795  550  SKS  2262  932  905  SKKS  2296  587  666  PcP  570  137  151  PKP  375  322  264  PKKP  1715  700  775  PcS  889  449  379  Table 2: Number of time windows used in this study, for each phase before selection based on cross-correlation coefficients, and after selection for the PREM and S20RTS models. Note the smaller number of PcP time windows before selection compared to ScS, which is in part due to P merging with PcP at smaller epicentral distances than for S and ScS.
where δr (in km) denotes the variation of the CMB radius at the bouncing points, p (in s/rad) is the ray parameter for the 1-D model without perturbations of CMB topography, and η = r cmb /v(r CM B ) (in s/rad) is the vertical ray parameter, with η + , and η − representing the values of η just above, and below the CMB, respectively. For PREM, which is our selected model, η + = 479.03 for S phases, and η + = 253.71 for P phases, while η − = 431.51 (only P phases are supported in the core).
Traveltime perturbations for the synthetics are computed by cross-correlation between synthetics for the reference model (either PREM, or S20RTS) and for the same model with added perturbations of CMB topography. For each of the seismic phases ScS, SKS, SKKS, PcP, PKP, PKKP and PcS, we first compute preliminary time windows from 5 s before, to 20 s after the arrival time of the phase as predicted using tauP [Crotwell et al., 1999] and PREM. We exclude time windows that overlap with the phases S, SS, SSS, sS, P, pP, PP, PPP, sP and pS, since overlapping phases affect the traveltime measurement by cross-correlation, which makes the comparison with ray theory less meaningful.
For a given reference model, we then compute improved time windows centered on the peak of each phase following the adaptive stacking algorithm of [Rawlinson and Kennett, 2004], in order to: 1) correct for the shift in travel-time for the case of a 3-D reference model, and 2) further eliminate time windows with overlapping phases. The improved time windows have a length that ranges from -15 s to +15 s after the arrival time of the peak for each phase. The output of the adaptive stacking algorithm is the time windows aligned on their respective peak amplitude and an average source wavelet that represents the apparent source time function of the phase of interest. An example of the adaptive stacking procedure is shown in Figure 4. The selection of time windows without overlapping phases is done by imposing a minimum cross-correlation coefficient between a given time window, and the inferred source wavelet. In the remaining of this section, the minimum cross-correlation coefficient is set to 0.8 for the P phases (PcP, PKP, PKKP and PcS), and to 0.95 for the S phases (ScS, SKS, SKKS), due to the higher signal-to-noise ratio for Sthan for P-waves. For the PcP phase, we further exclude all records at epicentral distances smaller than 40 • , since such records show consistent, large disagreement between RT and FW. The total number of time windows before and after selection is shown in Table 2. The amount of time windows that pass our quality criteria becomes significantly smaller, however, we prefer to base our conclusion on a more robust dataset. The smaller numbers of high correlation coefficients time windows in noise-free synthetics also shows that the isolation of these phases can be considerably cumbersome in real datasets.

Computation of Fréchet sensitivity kernels
The sensitivity kernels are computed for time windows of 30 s length, centered on the predicted peak arrival of the Pand S-wave phases under investigation. The sensitivity is computed based on the Born approximation [e.g. Marquering et al., 1999, Dahlen et al., 2000, Dahlen, 2005] and using adjoint methods as implemented in the software package SPECFEM3D_GLOBE for boundary [Dahlen, 2005, Liu and and volumetric kernels [Tromp et al., 2005]. With this analysis, we conduct an investigation of the traveltimes by visualising their sensitivity to volumetric, namely shear and compressional velocity mantle structure, and boundary parameters, this is interface denoting CMB. The main goal is to distinguish these phases which can truly help us to improve core-mantle boundary topography models in a full-waveform inversion scheme, judging by the effect this has on their traveltime.
More specifically, the kernels are computed with respect to the CMB radius and compressional and shear wavespeed in 1-D background model PREM. We compute the so-called "banana-doughnut" traveltime kernels , [e.g. Luo and Schuster, 1991, Marquering et al., 1999, Dahlen et al., 2000, which associate a traveltime measurement on the waveforms to the changes in the finite-frequency sensitivity depending on the background model. The traveltime refers to the time window around the arrival of the phase according to a given model, which we have fixed to PREM, for a single pair of source-receiver. The relationship connecting a traveltime window to volumetric structure defined by seismic properties of compressional (α) and shear (β) wavespeed are shown below. The following expressions are given for a reference seismic receiver (r) at position x r . The reader is referred to [Tromp et al., 2005] for the complete derivations according to the Born approximation and the adjoint method. The volumetric kernels are related to the time window via: The units of volumetric wavespeed kernels are s · km −3 .
The expression relating the traveltime window which encloses the seismic phase to the boundary perturbation on the core-mantle boundary (solid-fluid interface denoted by ΣF ) is based on theory derived by [Dahlen, 2005] and implemented in the spectral-element method: whereby K h is the time integration for variations along the boundary [Dahlen, 2005]. The units of boundary kernels are s · km −2 .
We select the phases for an ideal set-up of source-receiver given their epicentral distance also according to literature and optimal theoretical observational range. Specifically, for phases ScS, PcP, and PcS the epicentral distance is chosen to be 53 • , as at this distance they are less likely to interact with SS, PP, and P, respectively. The phases PKKP, SKS, and SKKS are observed at an epicentral distance equal to 110 • for better isolation of their theoretical arrival time. The PKP phase is selected at a distance of 160 • as, according to 1-D ray theory and tauP analysis, at this distance, it is better separated from PKIKP. The phases that travel through the outer core as well as PcS are time windowed on radial components of the synthetic seismograms, as they are expected to have maximum energy at this component, while core reflected phases are observed and analysed on the vertical (for PcP) and transverse component (for ScS). The traveltime curves are shown in figure 1 with annotated characteristic distances for each seismic phase window.
For the visualisation of the sensitivity kernels, we use a characteristic scale length, which is appropriate for the range of data values of the kernels and allows us to observe the theoretical finite-frequency sensitivity. It is helpful to consider that the volumetric and boundary kernels have different units and are explicitly calculated for different parameters; therefore the model update using these kernels will be different in the context of full-waveform inversion, depending also on the optimisation algorithm chosen for the potential model updates. Considering the expression given above, i.e. eq. 3 and eq. 4, it is not easy to directly compare boundary to volumetric kernels and assess how they will change the initial model. It is however meaningful to consider that the traveltimes of the selected seismic phases will likely show sensitivity to both parameters, as expected by the potential trade-off between velocity and CMB structure, and thus it is worthwhile to judge the qualitative effects of volumetric and boundary structure separately.

TRAVELTIME ANALYSIS USING FULL-WAVEFORM MODELLING
In this section, we show the comparisons of traveltime anomalies measured on FW synthetics to RT prediction for phases ScS, SKS, SKKS, PcP, PKP, PKKP, and PcS. In Figures 5 and 6, we present the results of this comparison for all topography models and for all the selected phases, for synthetic waveforms calculated using PREM, and S20RTS (1-D, and 3-D, respectively) as velocity background model. Figures 5 and 6 show that there is part of the measured traveltime shift which is not fully recovered when translating a topographic variation to time delay using ray theory. This essentially means that the ray theoretical expression cannot fully explain the topographic variation as time anomaly measured in full-waveforms with and without topography perturbation. The discrepancy seems to increase linearly with dt RT , and thus with the amplitude of the CMB topography, and varies significantly for different topographic models and phases. If the recovery using RT could accurately explain the traveltime shifts on these seismic phases, all points in Figures 5 and 6 would lie on the line with slope equal to one. However, the slopes of the linear regression lines have values different from one, as shown in Table 3. For the ScS phase (which gives the most corresponding results between RT predictions and WF measurements) with PREM as a background model, deviations from a perfect agreement are between 2% to 7%, while other well-resolved phases (SKS, SKKS, PcP and PKP) show deviations of up to 26% (for the SKS phase for model T1). This is true when the background velocity model is identical between the compared waveforms, namely with and without the topography model. In this case, we use PREM which is a 1-D model and is also consistent with the ray theoretical predictions made for the given topography model. Figures 5 and 6 show that the agreement between a ray theoretical prediction of time anomaly caused by topography and the equivalent measurement on full-waveform synthetics is better for S-wave phases than for P-wave phases. More specifically, ScS, SKS, SKKS (and PKP from the P-waves) traveltime shifts due to topography seem to be better Phase T1  TC1  TK10  LM91  phase T1  TC1  TK10 Table 3: Slopes of the linear regression lines for the dt RT vs. dt W F diagrams in Figures. 5 and 6. Each cell shows two slopes, which correspond to S20RTS (left side of the forward slash), and PREM (right side) as a background model, respectively.
predictable with ray theory. However, it should be noted that the difficulty in predicting an adequate time anomaly may not be solely due to ray theoretical limitations. Additional difficulties may also arise when trying to isolate the P-waves.
Focusing now on the comparisons of ray theoretical predictions to measurements on full-waveforms using S20RTS as a background model (blue dots on Figures 5 and 6), we see that the agreement between RT and WF is worse than when using PREM as a background model. Table 3 summarises the slope values, for example, for the ScS phase, deviations from a perfect agreement between RT and WF are ∼1.6 to ∼3.5 times larger (i.e., 5% to 18%) for S20RTS as a background model than for PREM, depending on the topographic model. The slopes vary significantly with the topographic model, and seem to be larger for models with bigger absolute topographic variation or smaller-scale amplitudes of CMB topography.
We further investigated this behaviour by computing traveltime residuals for three scaled versions of model LM91, with maximum topographic amplitudes of ±1 km, ±2 km, and ±4 km, in addition to the original LM91 model, whose maximum amplitude is ±8 km. These are shown in Figure 7. The results for the case of a 3-D background model show that the slope values are very close for all the scaled versions and the original model. This implies that the amplitude of the topography has some effect, but that it does not crucially affect the agreement between RT and WF, at least for amplitudes up to ±8 km. The differences in slopes for different topographic models could therefore be due to differences in the scale of the topographic anomalies. We note that S20RTS is a relatively smooth 3-D model compared to more recent global tomographic models. We expect even larger discrepancies for 3-D models with stronger, smaller-scale velocity anomalies (e.g., S40RTS [Ritsema et al., 2011] or SEMUCB-WM1 [French and Romanowicz, 2014]).
Our results illustrate the key role played by 3-D velocity variations in determining traveltimes. While ray theoretical predictions are based on a 1-D ray tracing expression, the result remains that the measurements on full-waveform synthetics show discrepancies with predictions, even though the measurement is done by comparing 3D and 3D+CMBtopo synthetics (keeping a fixed velocity model), indicating a non-linear effect of the two parameters to the traveltime. This also shows the highly non-linear effects of topography and 3-D variations on traveltimes of these phases, when topography is the only parameter added to the simulation. It is expected that some of the mismatch could properly be accounted for when a 3-D velocity correction is made, but it is not expected that it will fully be compensated, due to non-linear trade-off between CMB topography and the 3-D structure of the lowermost mantle, and to finite-frequency effects discussed in the next section. Figure 8 shows the complete sensitivity of the P-wave phases under investigation. The Fréchet derivatives are computed for compressional (V P ) and shear wavespeed (V S ) as well as for boundary of the interface corresponding to core-mantle boundary at approximately 2890 km for the chosen background model, PREM. All P-wave phases show their theoretical, finite-frequency, hence broader Fresnel zone, sensitivity to V P . This indicates that the time window chosen around the predicted arrival according to PREM encapsulates the phases well enough and they present little or no interference with other waves. It is worth noting that, specifically for the PcS and PcP phases, they seem to exhibit considerable sensitivity to shear wavespeed as well, showing the more complex contributions to their traveltimes and perhaps explaining the weakening of their sensitivity to topography, as shown in the analysis done in section 3.

FRÉCHET KERNEL GALLERY FOR CORE-MANTLE BOUNDARY RELATED PHASES
Regarding the S-wave phases, displayed in Figure 9, their traveltime sensitivity to shear and compressional wavespeeds and CMB surface are shown in a similar way as in Figure 8. It is readily observed that the shear wavespeed kernels computed for the time windows around the predicted arrival of these phases also clearly show their theoretical path, in this case broader and more extended due to the finite frequency effects which are in this case properly accounted for. The seismic phase SKS shows a broad and noisy sensitivity. The travelling paths of SKS, SKKS through the mantle are very similar, with their major differences appearing within the outer core. This is why their differential time has frequently been used for shear wave splitting studies in the lowermost mantle, [e.g. Restivo and Helffrich, 2006, Favier and Chevrot, 2003, Sieminski et al., 2008, Souriau and Poupinet, 1991. In the context of our study, it is important to understand their finite-frequency sensitivity and its variability between the two phases when it comes to the CMB surface. The S-wave phases also show clear contributions only from their theoretically predicted path in the time window selected for the calculation. In contrast to the P-wave kernels, there is little to none sensitivity of these phases to compressional wavespeed.
Turning to the boundary kernels, from visual inspection it becomes apparent that contributions to the traveltime sensitivity are low, especially for the traveltime windows of PKP and PKKP. This is in reference to the background model PREM. A low sensitivity to that means that the theoretical traveltime has low values at the interface of the core-mantle boundary, hence the particular seismic phase does not seem to arrive at a considerably different time than that reference traveltime prediction using PREM, due to the interaction with the core-mantle boundary and its topography. The boundary sensitivity of PcS and PcP are more pronounced and seem to sample its structure adequately and with the typical top-side reflection Fresnel zone. Considering the boundary sensitivity of S-waves, from the kernels in Figure 9, we observe that also for these phases, the reflected ScS phase has a clear sensitivity to the CMB. However, for SKS and SKKS, as shown by their corresponding boundary kernels, most of the sensitivity is observed at the piercing points in and out their way to the core, with SKKS having a stronger contribution. We discuss these results in greater detail in the next section.

DISCUSSION & CONCLUDING REMARKS
We have presented a traveltime analysis, where a ray theoretical prediction of traveltime anomaly caused by core-mantle boundary topography variations was compared to a measurement made by cross-correlation on full-waveform synthetics, imitating the common techniques of measuring time shifts on recorded waveforms to infer CMB topography. The retrieval of traveltime anomalies using ray theory shows that we could predict the shift caused by a given topography variation, albeit with some loss of accuracy and most reliably in time windows for the selected S-wave phases. We also showed that the correspondence worsens when the measurements are made in a 3-D background model, implying considerable effects of lateral velocity variations on the traveltimes. This is shown by a deviation from a perfect linear correspondence between the prediction and the measurement on realistic waveforms.
For the P-wave phases, the correspondence between ray theory and measurements on waveforms is much worse than for S-wave phases, except for PKP and PcP. The lack of agreement could be associated with better time window isolation of the given S-wave phases, contrary to the difficulties occurring for P waves. Specifically, multiple branches of PKKP can make its peak separation from its different branches (e.g. PKKPdf ) more difficult than for other phases. It is also generally expected that P-wave propagation is more complex when interaction with strong velocity contrasts occurs [Aki andRichards, 1980, Doornbos andMondt, 1980]. This involves both up-and down-going P and S waves when a P-wave interacts with an interface, and especially the CMB, due to the strong contrasts in density and thermo-elastic properties between the solid rocks defining the mantle and the liquid iron alloy determining outer core composition.
For both the case of PREM ( Figures 5 and 6 red points) and S20RTS ( Figures 5 and 6 blue points) as a background model, we observe deviations from a perfect agreement between traveltime residuals predicted by ray theory, and those measured using cross-correlation on realistic waveforms (slope value equal to unity). The deviations seem to increase linearly with the amplitude of the CMB topography. The rate of increase does not seem to depend on the amplitude of the CMB topography, but varies with topographic models, which could indicate a dependence on the scale of topographic anomalies. In particular, the rate seems larger for CMB topography models T1 and TC1 [Deschamps et al., 2018], which have smaller-scale anomalies, but are not built on seismic data and whose geographical distribution may thus differ significantly from the real CMB topography. The deviations for the case of PREM as a background model are between 2% to 7% for the ScS phase (the phase with the most reliable measurements), and between 5% and 18% for the case of S20RTS as a background model.
The slopes given in Table 3 for each of the phases and each of the background model show that the sign of the topographic variation can be in most cases correctly defined. However, for the PKKP phase and all models tested in our study, the measurement done in both background models seems to fail to infer the correct sign of topography (negative slopes). The linear regression attempted to fit the points of measurement, though, is not efficient as there are many points for which the ray theoretical prediction is close to zero. This perhaps implies that PKKP is indeed difficult to isolate and given its multiple branches, an assignment to its traveltime shift due to CMB topography, neglecting effects from lowermost mantle and outer core velocity variations, is a rather simplified assumption as many effects can lead to almost synchronous arrivals [e.g. Earle andShearer, 1997, Rost andRevenaugh, 2003]. We note that difficulties isolating separate branches of the PKKP phase would be reduced when using higher-frequency waveforms, as usually done when using ray theory, [e.g. Tanaka, 2010]. This is, however, not currently possible when using our spectral-element synthetics, because of the significant increase in computational requirements for higher frequencies.
This implies that phases such as high-frequency PKKP are currently difficult to use in full-waveform inversion for the core-mantle boundary topography.
The smaller deviations when using PREM as a background model suggest that the ray theory expressions in eqs. 1 and 2 are not completely adequate for the time shift measurements on 3-D background. These are derived using the ray parameter for a 1-D model (i.e., without topography), which implies that the observed deviations could be an effect of the small change in ray parameter due to the addition of CMB topography. Most likely though, the deviation could be due to finite-frequency effects when using the full-waveform theory. The larger deviations observed when using S20RTS as a background model indicate trade-offs between the CMB topography, and 3-D Sand Pvelocity anomalies in the mantle. For model T1 with strong, small-scale topographic anomalies beneath cold downwellings at the CMB, the deviations reach 18% for the ScS phase. This could affect the accuracy of inferred global tomographic models, which often do not account for the core-mantle boundary topography, although deviations in our calculations are not so large as to imply significant artifacts in tomographic models.
The finite-frequency sensitivity kernels show that the phases we investigated have strong contribution in their predicted traveltime windows and are not interfering with other phases. This observation shows that measurements on these phases can contribute to a good quality traveltime dataset and misfit measurements on their time windows and could help for a meaningful model update, given that they can be well separated with time windows. We also observe that reflected P-wave phases, i.e. PcP and PcS, are sensitive to shear wavespeed structure, too. In common inversion methods for imaging CMB topography, traveltimes are corrected for mantle velocity structure. The traveltime sensitivity of the S-waves does not seem to present similar patterns to that of P-waves, namely the traveltime is affected almost only from shear wavespeed variations. This may partly explain the poorer agreement between ray theory and measured traveltimes for P-waves than for S-waves, as the ray theoretical prediction does not account for the influence of shear wavespeed structure on the P-wave traveltimes (for example see expression 1). It also implies that linearised corrections for mantle structure may require the correction for shear wavespeed structure as well, similar to a conclusion made by [Koroni et al., 2019] for the, rather difficult to observe, PP precursors.
The sensitivity of the traveltimes to CMB seems to be weak for most of the selected time windows, implying that the contribution of the CMB interface to the traveltime shift is too weak to provide a robust update of the topography. Notable exceptions are the phases which reflect off the CMB, namely the ScS, PcP, PcS and to a lesser extent SKS, SKKS at the piercing points on the way in and out the outer core. This indicates that for purely CMB topography retrieval using traveltimes of these phases, one should consider that their sensitivity differs significantly and some phases offer a more optimal illumination of CMB structure than others (for example compare ScS to SKS).
Trying to explain the discrepancies between predicted and measured traveltime anomaly when working in 3-D background, this time in PREM and S40RTS [Ritsema et al., 2011] (1-D and 3-D mantle models), we computed boundary sensitivity kernels for the same time window isolating the phase SKS, same station-source configuration as in Figure 6. This phase is chosen for exemplification. The resulting kernels show that the 3-D model in the background significantly affects the boundary sensitivity of the traveltime, see Figure 10. This indicates that 3-D effects are imprinted into the traveltime of seismic phases and their sensitivity to the CMB is clearly changing in a substantial way, thus making an interpretation of the traveltime anomaly due to only topographic variations difficult. Additionally, given that our knowledge for the lowermost mantle structure is still lacking sufficient resolution, a joint inversion using Fréchet kernels, where both velocity and CMB topography are updated, may be more appropriate for the nature of this problem. This also confirms the traveltime analysis comparing traveltimes calculated using ray theory to those measured on synthetics, and indicates that not addressing this trade-off could bias both the resulting velocity structure and the inferred CMB topography.
From our sensitivity analysis we can conclude that top reflected seismic phases (e.g. ScS) are adequately informative for properly imaging CMB topography using boundary kernels. We also find that the effect of the velocity variations are significant and should be taken into account using the volumetric finite-frequency kernels as shown in this work. To reliably separate the two effects seems complicated, especially with an one step approach; however, it can be more efficient to utilise an iterative optimisation of a starting model and help make better inferences about CMB topography. High resolution imaging of the CMB and the mantle structure just above this interface requires an interdisciplinary approach. Understanding how subduction processes and hot material accumulated at the bottom of the mantle can affect observations of seismic phases traveltimes on the recordings is fundamental for a clarified selection of seismic traveltime data with stronger sensitivity at these regions.  [Crotwell et al., 1999]. The phases we investigate in this work are drawn with solid lines. More primary phases are plotted in order to justify the chosen epicentral distances for the sensitivity analysis. The characteristic distances chosen for computing the sensitivity kernels are shown by blue circles on the solid lines for each phase. Right: Theoretical raypaths for the same model. Red triangles denote stations at the epicentral distances were we opt observing the phases for the sensitivity analysis. The colours in the two sub-figures are not corresponding. IC denotes inner core, OC denotes outer core.  [Rawlinson and Kennett, 2004] used to define the time windows, for the ScS phase for model S20RTS and event no. 1 in Table 1. (a) Waveforms before alignment, using time windows 5 s before, and 20 s after the ScS arrival; (b) same as (a), but after alignment using the source wavelet (red trace) obtained from the adaptive stacking procedure. P c S P c P P K P P K K P T K L M9 1 T 1 T C 1 Figure 5: Comparison between ray theoretical prediction of time shift given a CMB topography model to the time anomaly measured on each time windowed seismic phase by comparing synthetic waveforms in 1-D and 3-D velocity background, red and blue points on the plots, respectively. Each column corresponds to a different topography model (TK, LM91, T1, TC1). Each row corresponds to a seismic phase (PcS, PcP, PKP and PKKP). Most of the phases show acceptable correspondence between predicted and measured by cross-correlation time shift, implying that the signature of CMB topography can be only partially captured by using RT. This is not the case for phase PKKP, where a prediction fails completely and the time window isolation is insufficient to make a reliable measurement. Each row corresponds to a seismic phase (ScS, SKS and SKKS). Most of the phases show better correspondence between predicted and measured by cross-correlation time shift, implying that S-wave phases have a better sensitivity to CMB topography and this can be accommodated by RT.
P c S P c P P K P P K K P S c S S K S S K K S Figure 7: A test using four version of the same topography model to check whether the amplitude of the variation linearly affects the RT prediction and cross-correlation time shift measurements. The model LM91 is scaled for ±1 km (yellow dots), ±2 km (orange dots), ±4 km (green dots) and ±8 km (blue dots) topographic variations. The slopes of a simple linear regression for the cloud of points remain almost unchanged with the scaling of the model. This indicates that there is no clear dependence on the topographic variation when considering the agreement between prediction of RT and CC measurement.   Figure 10: Comparison between the boundary kernels computed in 3-D background using model S40RTS (red line) and the same seismic phase SKS in 1-D background PREM (black line). the fact that the boundary kernels differ due to the different background model shows that interpretation of the traveltime shift of the phases due to topography should take into account the trade-off between topography and 3-D velocity variation, ideally in a non-linear way (not only by subtracting a time shift correction due to the 3-D mantle model). The units of boundary kernels km −2 · s.