Broadband waveform modeling with the spectral element method on Earth Simulator for the study of the structure at the top of the Earth's outer core

Abstract

The effects of the complex seismic structure in the lowermost mantle on the seismic SmKS phases that propagate beneath the core-mantle boundary are important, but as yet unclear. Thus, in this study, broadband waveform modeling with the spectral element method is conducted using the Earth Simulator. One-hour length seismograms are first synthesized with one-dimesional velocity structure of PREM, and the portions of the SmKS phases are retrieved. The shortest period that the Earth Simulator can achieve is up to 3.5 s, which is too long to reproduce S5KS and higher SmKS phases. To read the differential travel times of SmKS phases accurately, the phase-weighted stack is adopted and the uncertainty is inferred with the bootstrap method. Next, wave fields are simulated with three-dimensional velocity structures of S20RTS with emphasized velocity perturbation at the base of the mantle and SB4L18 expanded by spherical harmonics. The Earth Simulator enables the generation of a three-dimensional (3D) structure using spherical harmonics coefficients of up to 40 degrees. The different models result in different residuals for differential travel times of S4KS-S3KS and S3KS-S2KS and change in the incident azimuths of S3KS with respect to S2KS, even if global tomography models with long-wavelength heterogeneity of several thousand kilometers are used. These results clearly suggest that there are strong effects of heterogeneity in the lowermost mantle on the differential travel times of S4KS-S3KS and S3KS-S2KS. The characteristcs of the uncertainty depend on the 3D-mantle models, which may provide clues to the separation of the effects of heterogeneity at the base of the mantle for SmKS anomalies.

1. Introduction

The seismic structure at the top of the Earth's core is quite important to elucidate the stable stratification in its outermost part, which is due to compositional and thermal origins (Lister and Buffett, 1998; Buffett and Seagle 2010; Buffett and Seagle, 2011). SmKS waves are most suitable for observational constraints, as these waves are propagated in the outermost core as shown in Fig.1. They traverse the mantle as S-waves, and are converted to P-waves in the outer core, reflected m-1 times at the underside of the core-mantle boundary (CMB), and finally return to the mantle as S-waves. Although many recent studies suggest that P-wave velocity in the outermost core is slightly slower than that in PREM (Dziewonski and Anderson, 1981), the thickness of the low velocity region and the amount of the velocity reduction are still controversial (Tanaka, 2007; Alexandrakis and Eaton, 2010; Helffrich and Kaneshima, 2010; Kaneshima and Helffrich, 2013). In these studies, the differential travel times between SmKS waves are generally analyzed to reduce the effects of mislocation and the structural uncertainty in the crust and mantle. However, it is thought that the effects from the strong and unknown structure in the ${\rm D}''$ region, the base of the mantle, certainly remain as a bias in the outer core structure estimation (Garnero and Helmberger, 1995; Sun et al., 2013). For further studies, relevant waveform modeling and comparisons with observed waveforms are required to determine how to distinguish the effects of the structures above and below the CMB. This study is a first step to addressing this issue using a waveform simulation and aims to suggest the requirements for the next update of the Earth Simulator.

Fig.1.

Seismic ray paths of SKS and SmKS waves.

2. Method

To date, the reflectivity method (RM) (Müller, 1985; Kennett, 1988) has been frequently used for the study of SmKS waves (Tanaka, 2007; Kaneshima and Helffrich, 2013). Here calculations have been conducted for a slowness range of 0.02-0.08 s km^－1, a sampling interval of 0.05 s, and an isosceles triangle source time function with a half duration of 1 s. The Earth flattening approximation (Müller, 1985) is applied to the PREM structure that is discretized typically at a thickness of 1 km. Physical dispersion due to anelasticity is considered. RM synthesizes SmKS waves down to a period of 2 s at distances from 150°to 160°in approximately 30 minutes, and can easily reproduce S5KS on velocity seismograms, as previously noted by Helffrich and Kaneshima (2010). To fully incorporate three-dimensional (3D) mantle heterogeneity, the spectral element method (SEM) (Komatitsch and Tromp, 2002) can be used on the Earth Simulator (Tsuboi et al., 2003), which is updated to the Earth Simulator 2 in 2009 and enables us to achieve the shortest period of approximately 3.5 s with 127 nodes (1014 CPUs). Here the setting by Tsuboi et al. (2003) is adopted, in which the sampling interval is 0.05 s and the impulse source time function is considered. Furthermore the waveform modeling for PREM is conducted by incorporating the topographic model ETOPO5 (NOAA, 1988) and crust model CRUST2.0 (Bassin et al., 2000) to simulate seismograms from a South Fiji earthquake (hypocenter information of the event occurred on October 16, 2007) observed by European seismic networks as a single large array comprising 114 stations (Fig.2). The latitude, longitude, and focal depth of the event are 25.755°S, 179.530°E, and 509 km, respectively, and its CMT solution is obtained by the Global CMT project (http://www.globalcmt.org). To synthesize one-hour length seismograms, the SEM is executed three times because of the elapsed time limitaions, which costs approximately 18 hours in total. Subsequently, only the portions of SmKS waves are retrieved.

Fig.2.

Geographical distribution of epicenter (star) and 114 stations (triangles) used for the calculation of synthetic waveforms. The great circle path between the epicenter and each station is represented.

After adjusting the arrival times of S2KS phases using the cross-correlation method and normalizing the amplitude with its maximum for each trace, the phase-weighted stack (Schimmel and Paulssen, 1997) with a power index of 2 is applied. Differential travel times are determied by reading the peaks of concerned phases on the vespagrams obtained from the original and Hilbert-transformed waveforms to consider phase shifts between SmKS waves. To evaluate the uncertainty of the measurement, the Bootstrap method (Efron and Tibshirani, 1993) as per Tanaka (2004) is adopted. One-hundred combinations of the normal and Hilbert-transformed vespagrams are prepared with the record sections through bootstrap sampling, and the peak locations are read as on the vespagrams by the original record section. Finally, averages and standard deviations of all peak locations are used for the results.

Furthermore, additional runs are conducted by incorporating 3D mantle models. The first model is S20RTS (Ritsema et al., 1999), which is originally involved in the mesh generator of the SEM. However, the velocity perturbations in S20RTS are smaller than those in other shear-wave tomography models. In particular, the perturbations at the base of the mantle are quite small: thus, the spherical harmonics coefficients at the base of the mantle are replaced with those multiplied by a factor of three. To consider another model, the spherical harmonics expansion (Becker and Boschi, 2002) up to 40 degrees is applied to the velocity perturbations in the SB4L18 model (Masters et al., 2000), which is originally a block model. To utilize the vertical spline functions used in the mesh generator of the SEM, whose values are normalized to unity at each depth node, spherical harmonics coefficients are obtained for the velocity perturbations at the depths concerned. Although it is possible to prepare any maximum degree of spherical harmonics, the degree of 40 is determined to be the appropriate maximum number because of the elapsed time limitation in the execution of the mesh generator on the Earth Simulator.

3. Result

First, the results of the simulation with a one-dimensional (1D) structure are summarized. The differential travel times of S3KS-S2KS, S4KS-S3KS, and S5KS-S4KS for PREM calculated by the ray theory at the array center (epicentral distance = 155.195°), for a focal depth of 509 km, are 44.40 s, 14.62 s, and 6.39 s, respectively (Table1). These values are compatible with those measured on the vespagram obtained from the phase-weighted stack of the velocity seismograms synthesized with RM, whose residuals with respect to the ray theoretical values are +0.24±0.04 s, +0.25±0.04, and +0.07±0.02 s for S3KS-S2KS, S4KS-S3KS, and S5KS-S4KS, respectively (Fig.3, Table 2). However, the residuals obtained by SEM velocity seismograms are +0.24±0.20 s and +0.77±0.23 s for S3KS-S2KS and S4KS-S3KS, respectively. The residual for S4KS-S3KS is considerably larger than that obtained by RM (Fig.4, Table 3). The large discrepancy observed in the S4KS-S3KS times is possibly due to an invisible S5KS phase that contaminates the peak location of the S4KS phase. The uncertainties become larger than those in the case of RM potentially because of the noise caused by scattering from the crustal structure and topography (Fig.4a).

Table 1. Differential travel times for SmKS phases calculated by the ray theory with PREM

Phase pair	Differential travel time (s)
S3KS-S2KS	44.40
S4KS-S3KS	14.62
S5KS-S4KS	6.39

Fig.3.

(a) Velocity seismograms of S2KS and later phases synthesized by the reflectivity method. PREM is used. (b) The vespagram by the phaseweighted stack with a power index of 2 using the original waveforms in which the peaks of S3KS and S5KS are marked by white crosses. (c) The vespagram by the phase-weighted stack using Hilbert-transformed waveforms in which the peaks of S2KS and S4KS are marked by white crosses.

Table 2. Differential travel times for SmKS phases measured on the vespagram from velocity synthetic seismograms by the reflectivity method with PREM

Phase pair	Differential travel time (s)	Residual wrt ray theory (s)
S3KS-S2KS	44.64±0.04	＋0.24±0.04
S4KS-S3KS	14.87±0.04	＋0.25±0.04
S5KS-S4KS	6.46±0.02	＋0.07±0.02

Fig.4.

Same as Fig.3, except the spectral element method is used for the waveform calculation. PREM incorporating three-dimensional crust and topography models is used.

Table 3. Differential travel times for SmKS phases measured on the vespagram from velocity synthetic seismograms by the spectral element method with PREM including a three-dimensional crust model and topography

Phase pair	Differential travel time (s)	Residual wrt ray theory (s)
S3KS-S2KS	44.64±0.20	＋0.24±0.20
S4KS-S3KS	15.39±0.23	＋0.77±0.23

Next, the displacement seismograms of RM (Fig.5) and SEM (Fig.6) are examined. The differential travel times of S4KS-S3KS and S3KS-S2KS are approximately 44.7 s and 15.2 s, respectively (Tables 4 and 5). Those for RM and SEM coincide within approximately 0.10 s; however, these values are larger than the ray theoretical values by approximately 0.3 s and 0.5 s for S3KS-S2KS and S4KS-S3KS, respectively, which is probably due to the frequency components of the displacement seismograms. As the displacement seismograms contain longer period components than the velocity seismograms, the peak locations may be delayed. This suggests that the observation of the differential travel times on displacement seismograms should be compared with those measured on synthetic displacement seismograms, not with the ray theoretical values. Furthermore, the S5KS phase is not identified on both displacement seismograms by RM and SEM. Thus, the displacement seismogram is not suitable to search for S5KS and higher SmKS phases, and the current power of the Earth Simulator is insufficient for outermost core study using velocity seismograms with SEM. We require a finer mesh spacing in the waveform simulation to obtain higher frequency contents of greater than 0.5 Hz, as opposed to that of approximately 0.3 Hz corresponding to the currently available SEM setting on the Earth Simulator.

Fig.5.

Same as Fig.3, except for the displacement waveforms. The peak of S5KS is skeptical.

Table 4. Differential travel times for SmKS phases measured on the vespagram from displacement synthetic seismograms by the reflectivity method with PREM

Phase pair	Differential travel time (s)	Residual wrt ray theory (s)
S3KS-S2KS	44.69±0.02	＋0.29±0.02
S4KS-S3KS	15.09±0.06	＋0.47±0.06

Fig.6.

Same as Fig.4, except for the displacement waveforms.

Table 5. Differential travel times for SmKS phases measured on the vespagram from displacement synthetic seismograms by the spectral element method with PREM, including a three-dimensional crust model and topography

Phase pair	Differential travel time (s)	Residual wrt Ray theory (s)
S3KS-S2KS	44.67±0.07	＋0.27±0.07
S4KS-S3KS	15.20±0.08	＋0.58±0.08

Next, the simulation results with 3D structures are overviewed. Fig.7a presents the ray paths of SmKS and the background represents S-wave velocity perturbations from the modified S20RTS projected on a vertical cross section from Fiji to Europe, which shows the S-legs of SmKS phases traverse the Pacific Large Low Shear Velocity Province in the mantle. Fig.7b indicates the ray piercing points of S2KS, S3KS, and S4KS at the source-side of the CMB with the horizontal map of shear wave velocity perturbation at the base of the mantle. Fig.7c shows the synthetic waveforms of S2KS and later phases. Fig.7d shows the vespagrams obtained from the original and Hilbert-transformed seismograms. The differential travel times of S3KS-S2KS and S4KS-S3KS are 45.31±0.47 s and 15.16±0.52 s, respectively (Table 6). Interestingly, the uncertainties become considerably larger than the previous cases, which is clearly due to the effect of mantle heterogeneity. The residuals for S3KS-S2KS and S4KS-S3KS times with respect to the times measured on the vespagrams from the displacement seismograms by SEM with PREM are ＋0.64 s and －0.04 s. On the other hand, Fig.8 shows the case for the model SB4L18. The uncertainties of the differential travel times are also large (Table 7). The residuals for S3KS-S2KS and S4KS-S3KS times with respect to PREM are ＋0.25 s and ＋0.33 s, respectively. These values are completely different from those for S20RTS probably because the pattern and strength of the shear velocity anomalies are different.

Fig.7.

(a) Vertical cross section from the hypocenter to a station in Europe with seismic ray paths of SmKS. The background color indicates shear-wave velocity perturbations in the mantle calculated with S20RTS, but the coefficients at the base of the mantle are multiplied by a factor of 3. (b) Geographical distribution of SmKS piercing points at the core-mantle boundary (CMB) beneath the hypocentral area. Triangles, crosses, and circles represent the piercing points of S2KS, S3KS, and S4KS, respectively, at the CMB. Contours indicate shear-wave velocity perturbations with a 0.2％interval at the base of the mantle. (c) Displacement waveforms of SmKS synthesized by the spectral element method with S20RTS, but the coefficients at the base of the mantle are multiplied by a factor of 3. (d) The vespagrams of S2KS and later phases.

Table 6. Differential travel times for SmKS phases measured on the vespagram from displacement synthetic seismograms by the spectral element method with the modified S20RTS including a three-dimensional crust model and topography

Phase pair	Differential travel time (s)	Residual wrt SEM with PREM (s)
S3KS-S2KS	45.31±0.47	＋0.64±0.48
S4KS-S3KS	15.16±0.52	－0.04±0.53

Fig.8.

Same as Fig.7, except for the three-dimensional model of SB4L18 expanded by spherical harmonic functions with coefficients up to 40 degrees.

Table 7. Differential travel times for SmKS phases measured on the vespagram from displacement synthetic seismograms by the spectral element method for SB4L18 including a three-dimensional crust model and topography

Phase pair	Differential travel time (s)	Residual wrt SEM with PREM (s)
S3KS-S2KS	44.92±0.54	＋0.25±0.54
S4KS-S3KS	15.53±0.61	＋0.33±0.62

Additionally, the apparent incident azimuths of S3KS phases are measured by the phase-weighted stacking of the SEM displacement seismograms for three cases: PREM, the modified S20RTS, and the approximated SB4L18, as shown in Fig.9. The peaks of S3KS phases are observed at a lapse time of approximately 45 s, and the uncertainty is also inferred by the Bootstrap method. The obtained incident azimuth for PREM is 36°±5.5°, which unsurprisingly coincides with the back azimuth from the array center to the epicenter. However, those for S20RTS and SB4L18 are 29°±7.5°and 30°±6.4°, respectively, suggesting that both wave fronts rotate counterclockwise (Table 8). Considering Figs.7b and 8b, the dense ray cluster of S3KS traverses a low velocity region in which the western part is slower than the eastern part. This can make the wave front rotate counterclockwise. Additionally, the mode of the incident azimuths for the modified S20RTS, located at approximately 22°, does not coincide with the average value (Fig.9b, right panel), and the scattering in the peak azimuth is asymmetric and the largest among the three models. These effects are probably related to the stronger lateral variation in shear velocity at the base of the mantle in the modified S20RTS than that in the SB4L18 (Figs. 7b and 8b). This suggests that the characteristics of the observational uncertainties may be clues for the considering the effects of the mantle structure on the SmKS observations.

Fig.9.

(Left) Stacked power in dB as functions of incident azimuth and the lapse time from S2KS arrival obtained by SEM synthetic seismograms with (a) PREM, (b) the modified S20RTS, and (c) the SB4L18. The azimuth of the S3KS peak measured on the stacked power diagram constructed with the original record section is marked by a white cross. (Right) Histograms for the azimuths of S3KS peaks measured on the stacked power diagrams constructed with 100 record sections by the bootstrap sampling for (a) PREM, (b) the modified S20RTS, and (c) the SB4L18.

Table 8. Apparent incident azimuths of S3KS peaks

Model	Azimuth (°)
PREM	36.0±5.5
Modified S20RTS	29.0±7.5
SB4L18	30.0±6.4

4. Conclusions

Broadband waveform modeling with the SEM is conducted to investigate the structure at the top of the Earth's outer core using the Earth Simulator. First, attention was focued on the syntheized SmKS phases with the 1D velocity structure PREM. The shortest period of the SmKS phases is currently 3.5 s, which is too long to reproduce S5KS and higher SmKS phases on velocity seismograms. Next, we simulate wave fields using 3D velocity structures of S20RTS with emphasized velocity perturbation at the base of the mantle and SB4L18 expanded by spherical harmonics. The different models result in different residuals for the differential travel times of S4KS-S3KS and S3KS-S2KS, change in incident azimuths, and various characteristics of observational uncertainties. These results clearly confirm that there are strong effects of heterogeneity in the lowermost mantle on the differential travel times of S4KS-S3KS and S3KS-S2KS, even though global tomography models with wavelengths of several thousand kilometers are considered. Furthermore, the characteristics of the observational uncertainties vary from model to model, providing potential clues to the separation of the mantle effects from SmKS anomalies.

Acknowledgments

This study is supported by the project of the whole Earth elastic response simulation (Group ID G5009). The code of the spectral element method v.3.5 is executed on the Earth Simulator. I thank Seiji Tsuboi for coordinating this project and providing me with effective scripts. Many thanks to Masayuki Obayashi and the anonymous reviewer for valuable comments and suggestions, and Daisuke Suetsugu for his editorial efforts. This study is partly supported by Grants-in-Aid for Specially Promoted Research (No. 22000003).

References

•     Alexandrakis,  C. and  D. W.  Eaton (2010), Precise seismic-wave velocity atop Earth's core: No evidence for outer-core stratification, Phys. Earth Planet. Inter., 180 (1-2), 59-65, doi:10.1016/j.pepi.2010.02.011.
•     Bassin,  C. ,  G.  Laske , and  G.  Masters (2000), The Current Limits of Resolution for Surface Wave Tomography in North America, EOS, Trans. Am. Geophys. Un., 81, F897.
•     Becker,  T. W. and  L.  Boschi (2002), A comparison of tomographic and geodynamic mantle models, Geochem. Geophys. Geosys., 3, 2001GC000168.
•     Buffett,  B. A. and  C. T.  Seagle (2010), Stratification of the top of the core due to chemical interactions with the mantle, J. Geophys. Res., 115, doi:10.1029/2009JB006751.
•     Buffett,  B. A. and  C. T.  Seagle (2011), Correction to “Stratification of the top of the core due to chemical interactions with the mantle” , J. Geophys. Res., 116, doi:10.1029/2011JB008376.
•     Dziewonski,  A. M. and  D. L.  Anderson (1981), Preliminary reference Earth model, Phys. Earth Planet. Inter., 25 (4), 297-356.
•     Efron,  B. and  R. J.  Tibshirani (1993), An Introduction to the Bootstrap, Chapman and Hall/CRC.
•     Garnero,  E. J. and  D. V.  Helmberger (1995), On seismic resolution of lateral heterogeneity in the Earth's outermost core, Phys. Earth Planet. Inter., 88 (2), 117-130.
•     Helffrich,  G. and  S.  Kaneshima (2010), Outer-core compositional stratification from observed core wave speed profiles, Nature, 468 (7325), 807-810, doi:10.1038/nature09636.
•     Kaneshima,  S. and  G.  Helffrich (2013), Vp structure of the outermost core derived from analysing large-scale array data of SmKS waves, Geophys. J. Int., 193 (3), 1537-1555, doi:10.1093/gji/ggt042.
•     Kennett,  B. L. N. (1988), Systematic approximations to the seismic wavefield, in Seismological algorithms: Computational methods and computer programs, edited by D. J. Doornbos, pp. 237-259, Academic Press.
•     Komatitsch,  D. and  J.  Tromp (2002), Spectral-element simulations of global seismic wave propagation - I. Validation, Geophys. J. Int., 149 (2), 390-412.
•     Lister,  J. R. and  B. A.  Buffett (1998), Stratification of the outer core at the core-mantle boundary, Phys. Earth Planet. Inter., 105 (1-2), 5-19.
•     Masters,  G. ,  G.  Laske ,  H.  Bolton , and  A. M.  Dziewonski (2000), The relative behavior of shear velocity, bulk sound speed, and compressional velocity in the mantle, in Earth's Deep Interior: Mineral Physics and Tomography from the Atomic to the Global Scale, Geophysical Monograph vol. 17, edited by S. Karato, A. Forte, R. Lieberman, G. Masters and L. Stixrude, pp. 63-87, AGU.
•     Müller,  G. (1985), The reflectivity method - a tutorial, J. Geophys., 58 (1-3), 153-174.
•    NOAA (1988), National Oceanic and Atmospheric Administration (NOAA) product information catalog-ETOPO5 Earth Topography 5-minute digital model, Techinical report, U.S Department of Commerce, Washington D.C., 171 pp.
•     Ritsema,  J. ,  H. J.  van Heijst , and  J. H.  Woodhouse (1999), Complex shear wave velocity structure imaged beneath Africa and Iceland, Science, 286 (5446), 1925-1928.
•     Schimmel,  M. and  H.  Paulssen (1997), Noise reduction and detection of weak, coherent signals through phase-weighted stacks, Geophys. J. Int., 130 (2), 497-505.
•     Sun,  D. ,  D. V.  Helmberger ,  J. M.  Jackson ,  R. W.  Clayton , and  D. J.  Bower (2013), Rolling hills on the core-mantle boundary, Earth Planet. Sci. Lett., 361, 333-342, doi:10.1016/j.epsl.2012.10.027.
•     Tanaka,  S. (2004), Seismic detectability of anomalous structure at the top of the Earth's outer core with broadband array analysis of SmKS phases,　Phys. Earth Planet. Inter., 141 (3), 141-152.
•     Tanaka,  S. (2007), Possibility of a low P-wave velocity layer in the outermost core from global SmKS waveforms, Earth Planet. Sci. Lett., 259, 486-499.
•     Tsuboi,  S. ,  D.  Komatitsch ,  C.  Ji , and  J.  Tromp (2003), Broadband modeling of the 2002 Denali fault earthquake on the Earth Simulator, Phys. Earth Planet. Inter., 139 (3-4), 305-313, doi:10.[2]1016/j.pepi.2003.09.012.