Adjoint inversion of antipodal PKPab waveforms for P wave velocity anomaly at the base of the lower mantle

We perform an adjoint inversion by using antipodal PKPab phases to estimate the heterogeneous Vp structure at the base of the lowermost mantle. We have carefully examined antipodal stations with high S/N ratios during the past 30 years and selected 20 source-receiver pairs with the epicentral distances > 179.0 degree and the Mw < 7.0. We have used the spectral-element method on the Earth Simulator of JAMSTEC to calculate the synthetic seismograms for a heterogeneous mantle Vp model with the accuracy of period about 8 s (110 km wavelength at the CMB). We have set up the time window to retrieve the PKPab phase from the vertical component of both observed and computed seismograms to calculate the adjoint source to obtain the sensitivity kernels of Vp in the mantle. The computed Vp sensitivity kernel for each event shows the characteristic annulus pattern in the lowermost mantle, which covers a large area of the CMB. The twenty PKPab earthquake-station pairs in this antipodal study contribute the equivalent of about 3140 measurements at the CMB — compared with 1871 previously studied — and provide new data. Therefore, although the number of individual source-receiver pairs is not large, the summed sensitivity kernels of the PKPab phase for Vp structure at the base of the mantle may be sufficient to model heterogeneity of the Vp structure at the CMB. We summed each event kernel to set up a sensitivity kernel of Vp in the lowermost mantle and iterated the inversion to estimate a heterogeneous structure in D ″ . Although we have iterations that dominantly affect both South America and the South Pacific, the sum-mary final model shows features within South of Africa, South Pacific, SE Australia, and Central & South America. Keeping the Vs model fixed, we map the CMB Vp heterogeneity measured by the parameter R s,p = dlnVs/dlnV p and find qualitative, proximal correspondence with the degree-2 pattern of Large Low Velocity Shear Provinces observed in shear tomographic models: in the Pacific and Africa R s,p > 2.0, whereas the surrounding edges of the edges of the Pacific show R s,p < 2.0.


Introduction
The base of the Earth's mantle has long been recognized as nonhomogeneous from the early days of seismology (e.g., Bullen, 1963Bullen, , 1965)).With the deployment of the World-Wide Standardized Seismographic Network (Peterson and Hutt, 2014) in the early 1960's, global travel time data sampling the Core Mantle Boundary (CMB) were shared, and photographic records were hand-digitized for computer analysis.With the development in the late 1980's of the Global Seismographic Network (GSN, Butler et al., 2004) and the international Federation of Digital Seismograph Networks (FDSN, Suárez et al., 2008), a superabundance of CMB data sources were then available from body-wave travel times and waveforms, and Earth's normal modes (spheroidal and toroidal) for discovery of CMB features and application of inversion methods to derive 3D structure.Today, two features of the CMB standout out defining the heterogeneity of the CMB: Ultra-Low Velocity Zones (ULVZ) and Large Low Shear Velocity Provinces (LLSVP).
ULVZ's were discovered in the late mid-1990's (e.g., Garnero and Helmberger, 1995;Mori and Helmberger, 1995) and identified as such by Wen and Helmberger (1998).The discovery of a degree-2 pattern at the base of the mantle beneath the South Pacific and Africa was also recognized in the 1990's in modal inversions (S12, Su et al., 1994;and SAW12D, Li and Romanowicz, 1996), and later termed LLVSP (Deschamps et al., 2010).Subsequently, both of these global CMB features have been modeled in greater detail and considered in how the CMB region's chemistry, phase and lateral variation relate to mantle convection and plate tectonics(e.g., reviews in Yu and Garnero, 2018;McNamara, 2019).Additionally, McNamara (2019) notes, "To answer these questions will require continued additional observations from seismology.Each new seismic observation provides a small piece to that puzzle." Many data sets and inversion methodologies have been globally applied toward resolving the structure.These include-though by no means comprehensive-the following references which have contributed to our knowledge of the CMB: Kustowski et al. (2008), Li et al. (2008), Ritsema et al. (2011), Young et al. (2013), Chang et al. (2015), Koelemeijer et al. (2016), Bozdag et al. (2016), Hosseini et al. (2020), Muir and Tkalčić (2020), and Lei et al. (2020).To our knowledge, no antipodal data have been used in prior inversions, and thus new information is available from the CMB.
PKPab has been employed by Young et al. (2013) and Muir and Tkalčić (2020), comprising 1871 measurements.Each of these sampled two wavelengths at the core entry and exit path which constitute small circles at the CMB.Antipodal PKPab is not a ray, but rather comprises a ray sheet sampling the whole path, integrating all azimuths between source and receiver.Since the shortest wavelength is about 110 km, each antipodal PKPab samples the entire circumference of the small circle (~17,230 km), comparable to about 157 ray paths.The twenty PKPab earthquake-station pairs in this study contribute the equivalent of about 3140 measurements at the CMB-compared with 1871 previously studied-and provide new data.
For synthesis of antipodal data, we have utilized the Spectral-Element Method in prior studies of the structures at the Inner-Outer Core Boundary (Butler and Tsuboi, 2010, 2020, 2021).In this study we follow the adjoint inversion methodology of Bozdag et al. (2016) and Lei et al. (2020).

Antipodal seismic observation
The seismic wavefield recorded at the antipodal point of an earthquake is characterized by the annulus region surrounding the diameter which connects the epicenter and the antipode.These seismic waves, such as PKIIKP or PKPab, are focused near the antipodal point where the seismic energy is amplified-excepting only paths traversing a diameter through the earth, i.e., the PKIKP wave.Various studies have been conducted on the Earth's core structure using these characteristic properties of PKIIKP and related phases at the antipode (Rial and Cormier, 1980;Butler, 1986;Niu, 2008;Butler and Tsuboi, 2010, 2020, 2021;Cormier, 2015;Attanayake et al., 2018;Tsuboi and Butler, 2020).
In this study, we use the PKPab wave recorded at the antipode to investigate the three-dimensional structure of the CMB.Basically, we use data from earthquake-receiver pairs studied in Butler and Tsuboi (2021) -Tonga to a station in Algeria (TAM), Sulawesi to Amazon (PTGA), northern Chile to Hainan Island (QIZ), and central Chile to mainland China (ENH and XAN).In addition to these earthquakereceiver pairs, we reviewed antipodal stations with high S/N ratio during the past 30 years and selected additional source-receiver pairs with the epicentral distances >179.0 degree and the Mw <7.0.We included Spanish data (CART, ECAL EMUR, and E098) from earthquakes in New Zealand, Venezuela (SDV) from an earthquake in Java, United States (AGMN and RTBA) from earthquakes in the Indian Ocean, French Guiana (MPG) from an earthquake in Banda Sea, Indonesia (PSI) from an earthquake in Ecuador and Russia (ARU) from an earthquake near the Pacific-Antarctic Ridge.Fig. 1 shows raypaths of PKIKP, PKIIKP and PKPab.Fig. 2 shows each epicenter, and the antipodal observation point, and Table S1 shows the hypocenter and source-mechanism information of the earthquake used.We have 20 earthquake-receiver pairs in total used in the present study.

Sensitivity kernels by adjoint method
Seismic tomography provides an efficient tool to explore the Earth's internal structure by using travel times of body waves (e.g., Dziewonski, 1984;Dziewoński et al., 1977;Fukao and Obayashi, 2013), surface wave dispersion (e.g., Ritzwoller and Levshin, 1998;Ferreira et al., 2010;Ekström, 2011;Ritsema et al., 2011;Chang et al., 2015), and normal modes (e.g., Kustowski et al., 2007).Recent advances in numerical techniques to compute seismic waveforms in complex 3D structures and the development of massively parallel super computers have made it feasible to apply full waveform inversion to estimate 3D structure of the Earth.Adjoint waveform tomography (Tromp et al., 2005;Fichtner et al., Liu and Tromp, 2008;Tape et al., 2007) is one such technique.
In adjoint waveform tomography, synthetic waveforms for realistic 3D Earth structure computed by numerical technique, such as spectralelement method, are combined with the computation of sensitivity kernels of waveform misfit to the 3D seismic velocity structure of the Earth.These sensitivity kernels will be used in an iterative inversion to update a model, 3D velocity structure of the Earth without using the raytheoretical approximation in the tomography study.Adjoint waveform tomography has been applied to regional, continental, and global scales to improve the 3D images of the Earth and demonstrated its potential and future directions (e.g., Bozdag et al., 2016;Fichtner et al., 2008Fichtner et al., , 2009;;Lei et al., 2020;Tape et al., 2009Tape et al., , 2010;;Tromp, 2020;Rodgers et al., 2022).
Here, we calculated the sensitivity kernel of P-wave velocity to be used in the waveform inversion for heterogeneous P-wave velocity structure at the base of the lower mantle.To calculate the sensitivity kernels, we use the Spectral Element Method (Komatitsch and Vilotte, 1998;Komatitsch et al., 2002;Tsuboi et al., 2003;Komatitsch et al., 2005), as in our previous studies (Butler and Tsuboi, 2010;Tsuboi and Butler, 2020;Butler and Tsuboi, 2020).The Earth's internal structure model used incorporates a 3D tomography model SGLOBE-rani (Chang et al., 2015) for the Earth's mantle, a crustal model CRUST2.0(Bassin  (Tsuboi and Butler, 2020).PKIKP (black) travels a ray path along the diametral axis between earthquake and seismic station, but is not antipodally focused.PKP-AB (blue) traverses the D″ region at the base of the Mantle and the upper Outer Core.PKP-C diff (red) diffracts around the Inner Core at the base of the Outer Core.PKIIKP (green) enters the Inner Core and reflects once beneath the top of the Inner Core.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)et al., 2000), and ellipticity.Topography at the surface of the Earth, rotation and the effect of the ocean are not included.SGLOBE-rani is a global radially anisotropic shear wave velocity model with radial anisotropy allowed in the whole mantle.It is based on a seismic data set of over 43 M seismic surface wave (fundamental and overtones) and body wave measurements.It models simultaneously crustal thickness and mantle structure, and its reference model is PREM (Dziewonski and Anderson, 1981).We use the anisotropic version of SGLOBE-rani and the scaling relation δVp/Vp = 0.5 δVs/Vs (Chang et al., 2015) to model Vp structure for the mantle and used as a starting model of the inversion.In the present study, we do not use transverse anisotropy to model the lowermost mantle structure but we may include this option in the future inversion.
The spectral element method program, SPECFEM3D_GLOBE, is used, and the number of elements at the surface along the two sides of block, NEX_XI, is set to 640.We use 9600 cores of the Earth simulator to calculate the theoretical seismic waveform which is accurate to a shortest period of about 6.9 s (about 8 s when bandpass filtered).We apply Chebyshev type-2 filter with the number of pole = 4, corner freq.0.02 Hz and 0.125 Hz to both data and theoretical seismograms.We have applied two-pass filtering for a zero-phase filter.The comparison of the theoretical seismogram computed by the spectral-element method and the observation is shown in Fig. 3, which is for June 16, 2000 Chile-Argentina border region earthquake (Mw 6.4) recorded at Xi'an, China (A) and January 25, 2014 earthquake in Java (Mw 6.1) and the receiver is the station SDV, Santo Domingo, in Venezuela (B).Because the epicentral distance for this station is about 179.16 degree and 179.09 degree, respectively, we see the large amplitude of PKPab phase due to the focusing at the antipode.The theoretical seismogram computed for this station by the spectral element method models this focusing generally well and the agreement between the observed seismogram and theoretical seismogram is satisfactory.However, for SDV, there exists a small difference in arrival time of the PKPab phase.This small difference in travel time of the PKPab phase may originate from the Vp model in the lower mantle used to compute theoretical seismograms.Therefore, we may retrieve information on hetereogeneous Vp structure in the lowermost mantle by using the observation of the PKPab phase at the antipode of the earthquake.
We compute the finite frequency sensitivity kernel of the P-wave velocity structure by using the adjoint method (Tromp et al., 2005(Tromp et al., , 2008;;Liu andTromp, 2006, 2008;Tape et al., 2007Tape et al., , 2009;;Peter et al., 2011).Computation of the finite frequency sensitivity kernel is done in two steps (see method for details).First, we use the global CMT solution (Ekström et al., 2012) of the earthquakes shown in Table S1, to calculate by using the spectral element method the theoretical seismic waveform for the seismic station at the antipodal point corresponding to each earthquake.Then, at the end of the time step, the global wave fields in which the theoretical seismic waveform is calculated are saved in the disk.We apply the same band pass filter to both theoretical seismograms and observed seismograms.From the calculated seismic waveforms and observed waveforms, we set the time window of arrival of PKPab waves and cut out these waveforms.We calculate the adjoint source for the travel time of the PKPab phase from the difference between the observed waveform and the theoretical waveform by using MEASURE_ADJ program (Tape et al., 2009) provided with the SPECFEM3D_GLOBE program package.Selection of an adjoint source to be used in the iterative inversion depends on what type of measurement is done for the computation of the misfit function.We calculate a cross-correlation of a travel time difference for an event kernel.The antipodal stations we have selected in the present study locate very close to the antipode (epicentral distance is 180 degree) but there are slight deviations from the antipode.Therefore, when cutting out the waveform, we need to include both major arc PKPab and minor arc PKPab arrivals and should set the window to include both of these phases.Then, using the adjoint source obtained in this way, the theoretical seismic waveform that propagates backward from the observation station to the seismic source is calculated by reversing the time.We calculate the finite-frequency kernel of the P-wave velocity in the lower mantle for the travel time differences of the PKPab wave for the path combining the propagation of the seismic waves from the source to the observation station used from the theoretical seismic waveform calculated beforehand.Tromp et al. (2005) has claimed that the sensitivity of the event kernel becomes high at regions near the source and receiver.To accommodate this characteristic, we normalize the event kernel by the value at deeper structure such as the lower mantle, keeping the sensitivity for the target region at the base of the lower mantle.
We plot an example of the cross-sectional view of the sensitivity kernel of the P-wave velocity structure in the lower mantle for the travel Fig. 2. The locations of seismic stations (blue circles) used for the inversion and their corresponding antipodal earthquakes (red stars, see Table S1).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)time differences of the PKPab wave in Fig. 4A.The event is the January 25, 2014 earthquake in Java and the receiver is the station SDV, Santo Domingo, in Venezuela.We show the sensitivity kernel for the minor arc of the propagation path.In Fig. 4B we show both the theoretical seismic waveform and the observed waveform for this earthquake-station pair.The adjoint source of the travel time differences of the PKPab wave is shown in Fig. 4C.The Fig. 4A demonstrates that the kernels are sensitive along the raypaths of PKPab phases especially at the core-mantle boundary.The adjoint kernel calculated for each earthquake and its antipodal observation point is called an event kernel.Each of the twenty earthquake-receiver pairs used in this study covers a ray-surface encompassing nearly 10% of the core-mantle boundary.Accounting for redundant coverage, most of the core-mantle boundary is substantially sampled, even though the number of the earthquake-receiver pairs is not large.By adding the event kernels calculated for each earthquakereceiver pair, the gradient of the misfit function to improve the initial model was calculated (see Method in Supplementary Information for details).

Inversion of PKPab for the CMB Vp structure
We use the gradient of the misfit function to find the change δm that updates the model (see Method for more information).There are several options for the schemes to calculate updates of the model, such as conjugate gradient method and steepest decent method.Here we use a steepest descent method with a step length determined by the given maximum update percentage.We found that the amplitude of sensitivity kernels for Vp structure in the lower mantle is generally small, which means that the model update is small in each iteration, probably because the heterogeneity in the lower mantle is not so large.Then we assume that the standard technique, such as steepest decent method is sufficient to find the updates of the model.Here, assuming that the initial model is updated with a step length of 3% for each iteration, the amount of modification of the model was calculated.When we update the model and calculate the synthetic by forward simulation, we calculate crosscorrelation travel time difference for PKPab phase between the data and synthetic for each station.Then we compute the root mean square (RMS) of the cross-correlation travel time differences for all stations to check the performance of the updated model.
We have iterated seven times to update the Vp model at the base of the lower mantle.Fig. 5 shows the variation of RMS computed for all stations after each iteration.We have aligned the data and synthetics with the arrival of PKIKP and measure the cross-correlation travel time differences for PKPab.After the fifth iteration, we have realigned the PKIKP because the arrival time of PKIKP phase in the synthetic seismograms slightly vary due to the modification of Vp structure in the lower mantle.Fig. 5 shows that RMS of the travel time differences for PKPab decreases reasonably well after the iteration.We also have shown in Fig. S1 the amount of modification of the P-wave velocity structure to be made at the base of the lower mantle after each iteration.These figures show that the distribution of sums of the event kernels for each iteration and may be used to examine which portion of the lower mantle the modification is made with the inversion.The shape of the kernel shows that the sensitivity is generally high in the southern hemisphere, which may be due the distribution of earthquake-receiver pairs used in this study.As we did not apply any weights in the summation of event kernel, the sum of event kernels may reflect the earthquake-receiver distribution.We may discuss the location of those features at CMB, such as ULVZ, and Fig. S1 demonstrates that the resolution is sufficient to discuss the reliability of the model for southern hemisphere.Fig. S2 shows comparisons of the data and synthetics computed for the Vp model after the inversion.Fig. S2 demonstrates that the updated model generally improves the arrival of PKPab in most of the stations.There are some stations which do not show improvements of waveform match after the iteration.These cases may be characterized by the inaccurate modeling of the earthquake source mechanism.Although we have selected events so that we do not include large earthquakes to avoid complex source mechanisms, we may need to consider reducing the magnitude threshold to avoid complex earthquake sources or determination of moment tensor solution for these events in the future study.

Results of the inversion
We show our results of inversion for Vp structure at the base of the lower mantle in Fig. 6.We plot (A) (Vp-Vp0)/Vp0, where Vp0 is the spherically symmetric model.In addition, we plot (B) Vp(SGLOBE-rani)-Vp0)/Vp0 and (C) (Vp-Vp(SGLOBE-rani))/Vp0 for comparison.We use PREM as Vp0 and the model plot is made for CMB.The grid data used to make this plot are also provided as supplementary data with this paper.The heterogeneity at the base of the mantle is characterized by ULVZs (e. g., Garnero and Helmberger, 1995;Mori and Helmberger, 1995;Wen and Helmberger, 1998) and LLVSP (Su et al., 1994;Li and Romanowicz, 1996;Deschamps et al., 2010) beneath the South Pacific and Africa.LLSVP is seen in the SGLOBE-rani Vs model and also seen in the Vp model, since Vp model is derived by the scaling relation.We use this scaled Vp model in the lower mantle as the initial model of our inversion.Our final model also shows these large scale features at CMB.Our model demonstrates that these large scale features of Vs at CMB also exist in Vp heterogeneity at CMB.Because the sensitivity kernel of our inversion has a good resolution in the southern hemisphere at CMB, the similarity of Vp heterogeneity with Vs at CMB may be reliable.Also, in the present study for Vp model at CMB, we use PKPab phase, which is not used in the SGLOBE-rani Vs model inversion.Use of a different dataset and the good resolution of the sensitivity kernel at CMB indicate that the large scale heterogeneity at CMB also exists for Vp heterogeneity.

Discussion
We use PKPab phase observed at the antipode of the earthquake to study the heterogeneous Vp structure at CMB.As we mentioned in the introduction, antipodal PKPab is not a ray, but rather comprises a ray sheet sampling the whole path, integrating all azimuths between source and receiver.It is necessary to use full waveform theory to model propagation of PKPab phase to get the theoretical waveform at the antipode.We use the spectral element method (e.g., Komatitsch et al., 2005) to compute theoretical waveform of the PKPab phase for 3D Earth model.The spectral element method is a widely used numerical scheme to compute theoretical seismograms for a realistic 3D Earth model.It combines the flexibility of the finite-element method with the accuracy of the pseudospectral method and is suitable for large scale computation on parallel supercomputers.We use 9600 cores of the Earth Simulator at JAMSTEC to compute theoretical waveform of PKPab with the accuracy of 6.9 s.This accuracy may not be sufficient to model fine structure of CMB with the accuracy of body waves, such as 1 s, but may be used to estimate the smooth structure at CMB.We try to estimate the limit of the accuracy which we will be able to achieve on this machine.The current Earth Simulator is equipped with AMD EPYC 7742 (64 cores) and NEC SX-Aura vector engine.Here we use AMD EPYC for our study and 9600 core amounts to 10.4% of the total system.Because we use 20 earthquake-receiver pairs, we need to compute the sensitivity kernel for Vp for each earthquake-receiver pair during the iteration.It took about 6 h to compute the forward and backward simulation to calculate one event kernel.We could use up to about 40% of the total system of the Earth Simulator simultaneously, which means that we may compute 4 event kernels at one time.Therefore, we may obtain 20 event kernels for each earthquake-receiver pair within at most 2 days, considering the I/O time and waiting time for computational jobs to execute.If we perform iterative inversion to update the initial model, we may obtain the final model within a half month, if we iterate 10 times, for example.We may conclude that this is a feasible amount of computation for our supercomputer system.
Then, we examined the minimum period of the accuracy achievable in our adjoint inversion by using our supercomputer system.We may use about 40% of the total system at one time for each user.Thus, we may use four times larger cores for our simulation at one time, which is 38,400 cores and the number of divisions used in SPECFEM3D_GLOBE, NPROC, is 80 in this case.If we may use this number of cores at one time, NEX of SPECFEM3D_GLOBE can be 1280 and the minimum period of the accuracy can be 3.4 s.For this number of cores, it may not be possible to compute event kernels for multiple earthquake-receiver pairs, and we need to wait about 6 h to obtain one event kernel.It will take one to two weeks to compute event kernels for 20 earthquake-receiver pairs, and we may be able to obtain the updated model within four to five months at most for 10 iterations.This is a feasible number and it should be worth trying for our next work.To get two times more accurate results, which is 1.7 s accuracy, we may need four times larger cores, 153,600.We do not have these resources at this moment but we should mention that this is not an unrealistic number.We may conclude that the adjoint inversion for CMB Vp structure with 1 s accuracy is within our reach.We also should mention the use of General Purpose Graphic Processing Units (GPGPUs) for scientific computation.Since the Earth Simulator is basically equipped with CPUs and vector engine, we have not tried the use of GPGPUs for our inversion.However, it is well known that the GPU version of SPECFEM3D_GLOBE is thoroughly optimized for the use of GPGPU environment.We know that the use of GPGPUs enables the computation by almost two orders of magnitude faster than the CPUs.If we use 9600 GPGPUs at one time for our computation, we will be able to  compute one event kernel with 6.9 s accuracy within several minutes.Therefore, we estimate that the adjoint inversion of the current scale can be done within a few days.We do not have the computer resources with this number of GPGPUs but we suppose that we should consider this possibility seriously in our future works.
We consider only isotropic Vp structure at the base of the mantle for simplicity.The starting model which we used for our adjoint inversion is SGLOBE-rani (Chang et al., 2015), an anisotropic Vs model throughout the mantle.Since transverse isotropy may be used to distinguish between horizontal and vertical mantle flow, it may be helpful if we can provide a transverse anisotropy model at the base of the lower mantle.Although there is no scaling relation suggested for Vp anisotropy for SGLOBE-rani, we may calculate sensitivity kernels for Vph and Vpv without additional computational cost.We may use a scaled isotropic Vp structure as a starting model and may compute sensitivity kernels for both Vpv and Vph.We have not checked if the horizontal components of observed PKPab phase at antipodal stations match with the theoretical seismograms, and should examine the potential use of three component adjoint sources for adjoint inversion.Dziewonski and Woodhouse (1987) introduced dlnVs/dlnVp "as a new in situ datum for mineral physicists," that in the mantle may be between 2.0 and 2.5, and may arise from lateral temperature differences associated with convective plumes.Since this introduction, dlnVs/dlnVp ≡ R s,p has been linked with many geophysical and mineral physics parameters and processes.These broadly include: partial melt fraction (Agnon and Bukowinski, 1990); a chemical component to large scale heterogeneity (Romanowicz, 2001); heterogeneity due to subducted slabs in the deep mantle (Saltzer et al., 2001); combined harmonic, anelastic and thermal effects on velocity and density (Komatitsch et al., 2002); the pattern of mantle convection in the deep mantle where compositional heterogeneity is away from major downwellings (van der Hilst, 2002); constraining spatial variations in temperature and chemical composition not attributable to temperature effects alone (van der Hilst et al., 2003); an indicator of textural equilibrium in pore aspect ratios (Takei and Nakajima, 2003); Large Low Shear Velocity Province (LLSVP) and Ultra Low Velocity Zone (ULVZ) models (Song et al., 2010;Cottaar and Romanowicz, 2011;Ritsema et al., 2014;Muir and Brodholt, 2014;Tanaka et al., 2015); the presence of the post-perovskite (pPv) phase in the lower mantle (Koelemeijer et al., 2016); thermal and compositional evolution of the earth when combined with insights from mineral physics (Su and Romanowicz, 2021).
In Fig. 6c we presented the difference between our final model and the initial SGLOBE-rani model.2016) R s,p ≈ 1.8 km at CMB for SP12RTS degree-12 model.More quantitative measures of R s,p are considered for future inversions in the context of radius above the CMB and the scaling for R s,p .

Summary
In this study, we applied adjoint inversion method for the antipodal PKPab phase to estimate heterogeneous structure of Vp at core mantle boundary.We used spectral element method to compute theoretical seismograms with the accuracy of 6.9 s and compute sensitivity kernels for Vp to achieve iterative inversion.The result of our inversion shows that there exists LLPVP in the Vp heterogeneity model (i.e., Large Low P Velocity Provinces) at approximately the same LLSVP location as suggested by Vs model, but with a Vp amplitude slightly larger than expected.This result shows that the adjoint inversion method using spectral element method with the accuracy of 6.9 s is feasible, although the required computer resources are moderately large.Because antipodal PKPab is not a ray, but rather comprises a ray sheet sampling the whole path, integrating all azimuths between source and receiver, we need to use full wave propagation scheme to model waveforms and compute sensitivity kernels.By using those phases, such as PKPab at the antipode, which require full wave propagation procedure for the study of Earth's deep structure, we may expand the resolution of heterogeneity models.In this regard, the use of the adjoint inversion method becomes promising when coupled with the numerical technique of the spectral element method to model seismograms.We also have shown that the projection of available computer resources in near future will enable us to perform the adjoint inversion with the accuracy of around 1 s for our cases.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data were obtained from GEOSCOPE and the Earthscope Data Management System.We used the computer program (SPECFEM3D) for Spectral-Element Method.Centroid moment tensor solutions (GCMT) are used for synthetic models.All data were downloaded from the Earthscope Data Management System and ORFEUS.Earthquake parametric data were downloaded from the USGS earthquake catalog and tabulated in the Supplementary Information file.Earthquake source mechanisms were downloaded from the Global Centroid Moment Tensor database.The datasets used and/or analyzed during the current study is available from the corresponding author on reasonable request.

Fig. 1 .
Fig. 1.Principal antipodal ray surfaces are shown in cross section(Tsuboi and Butler, 2020).PKIKP (black) travels a ray path along the diametral axis between earthquake and seismic station, but is not antipodally focused.PKP-AB (blue) traverses the D″ region at the base of the Mantle and the upper Outer Core.PKP-C diff (red) diffracts around the Inner Core at the base of the Outer Core.PKIIKP (green) enters the Inner Core and reflects once beneath the top of the Inner Core.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .
Fig. 3. Comparison of observed (black) and theoretical (red) seismograms for (A) June 16, 2000 Chile-Argentina border region earthquake recorded at Xi'an, China and (B) January 25, 2014 earthquake in Java and the receiver is the station SDV, Santo Domingo, in Venezuela.Broadband vertical component displacement seismograms are bandpass filtered by Chebyshev type-2 filter with the number of pole = 4, corner freq.0.02 Hz and 0.125 Hz. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig. 4. (A)A cross-sectional view of the event sensitivity kernel of the P-wave velocity structure in the lower mantle for the travel time differences of the PKPab phase.The event is the January 25, 2014 earthquake in Java and the receiver is the station SDV, Santo Domingo, in Venezuela.We show the sensitivity kernel for the minor arc of the propagation path.(B) The theoretical seismic waveform (red) and the observed waveform (black) for this earthquake-receiver pair.To calculate the adjoint source, we set the window between 1300 s and 1335 s, which are shown as the vertical bar.(C) The adjoint source is shown for the PKPab phases.These traces are ground displacement and bandpass filtered between periods 8 and 50 s.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .
Fig. 6.Map of the Vp velocity modification.(A) Our model, (B) initial SGLOBE-rani, (C) differences.Note that our model (A) is enhanced with respect the SGLOBErani (B).For example, for SE Australia shows dVp/Vp ~ 0.005 (A), ~ − 0.002 (B), and ~ 0.007 (C), where C = A-B.From the velocity variations in (C), we derive the R s,p = dlnVs/dlnVp = 2.0 for no change (white); otherwise blue implies R s,p < 2.0, and red implies R s,p > 2.0.Note in C also the apparent degree-2 pattern for LLVSPs in the Pacific and South of Africa.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
In this map, the zero values (white) indicate no Vp change, and hence R s,p = 2.0.Blue values indicated increasing Vp which corresponds for decreasing R s,p < 2.0; Red values indicate decreasing Vp corresponding with increasing R s,p > 2.0.The LLVSPs in the Pacific and Africa show R s,p > 2.0, whereas the edges of the edges of the Pacific show R s,p < 2.0.Qualitatively, some of the published R s,p values have proximal fits with Fig. 6c.These include Li et al. (2004) R s,p ≈ 2.1 for lower part of the lower mantle; (Davaille et al., 2005) R s,p ≈ 2.5 at the bottom of the South African mantle; Cottaar and Romanowicz (2011) R s,p ≈ 2.1 from P diff for the northern boundary of Pacific Superplume; and Koelemeijer et al. (