以SmKS波到時探討地球外核頂部的速度結構

(1)

國立臺灣大學理學院地質科學研究所碩士論文

Department of Geosciences College of Science

National Taiwan University Master Thesis

以 SmKS 波到時探討地球外核頂部的速度結構

Using SmKS Traveltimes to Investigate the Velocity Structure near the Top of the Earth’s Outer Core

唐楚欣 Vivian Tang

指導教授：趙里博士洪淑蕙博士

Advisors: Li Zhao, Ph.D. Shu-Huei Hung, Ph.D.

中華民國 103 年 7 月

July 2014

(2)

誌謝

能完成這篇論文，我要特別感謝我的指導教授趙里老師與洪淑蕙老師。從在中研院當暑期生就開始跟著趙里老師學習，趙老師非常的用心指導我，不管是在研究還是課業方面，對於我的問題總是耐心、細心的回答，並且培養我對研究的熱情與信心，對於老師給予我的幫助實在難以用言語表達我內心中的感謝，真的非常謝謝老師!到台大後有了洪淑蕙老師一同指導我，讓我能從不同的方向去思考研究問題，感謝洪老師與我討論並且提供建議，幫助我更透徹的了解問題，也使我的論文更完整。謝謝喬凌雲老師、郁文哲老師和龔源成老師共同擔任我的碩士口試委員，替我批改論文，提出不同觀點的意見，也給予我鼓勵。此外，也非常感謝欣穎、侑頻、柯彥廷和珮如等實驗室的強者學長姐們對我的幫忙與鼓勵，在我遇到問題時總能向他們尋求到答案。羅翊菁、林姿綺和其他的學弟妹，謝謝在我忙碌之時，願意幫我買飯或是拿食物餵食我，也陪我聊天舒壓，讓研究之路多了些趣味性。幸運的是在這裡我還有一群從大學時期一起奮鬥努力的好朋友相伴，李詩婷、吳欣玫、陳力維、秦念祺、葉德揚和蔣盼忻，以及其他在中央大學的好友們，雖然無法一一列舉，但心中時常覺得有這些朋友一起上課、討論問題、

吃飯聊天和互相鼓勵，真的是難得可貴的情誼。尤其是李詩婷和我從大學時期不

僅是摯友也是室友，在台大過的幾乎是24 小時形影不離的生活，特別謝謝她在身

邊，不論何時總是帶給我歡笑，我想我畢業之後會很懷念她在房間唱的那些台語歌! 還有謝謝其他在我周邊給我支持鼓勵的好朋友們，為了配合我的時間，約吃飯地點永遠是公館跟師大。

最後，我要感謝我的父母與妹妹，從小到大在爸媽細心關懷與照料下樂觀開朗的成長，會有今天的我，都要感謝他們的支持與栽培。謝謝爸媽的信任給我很大的發揮空間，也在我趕進度忙碌不能回家時總能體諒。當我心情低落遇到困難時開導我，還有時常北上給我鼓勵，回家也總是有滿滿的一桌好菜替我補身體!爸媽對於我的付出，是我繼續努力的動力。謝謝所有曾經給過我幫助和鼓勵的家人、

老師以及朋友們，沒有大家就不會這般順利且愉快的研究生活!

(3)

中文摘要

地球的固態內核之主要成份為鐵和鎳，另外還有少部分的氧、硫和矽等之比重相對較輕的物質。動力學研究認爲這些較輕的物質會因浮力作用上浮而聚集在地球液態外核的頂部，但目前零星的地震學觀測結果尚未能夠對此說法給出明確的證據。從地震波觀測的角度來看，地球外核頂部若聚集有此較輕的物質，將會使核幔邊界(CMB)以下一部分深度區域内的密度和 P 波（即外核内 K 波）速

度低於它們在該深度上的全球平均模型 PREM 的速度值。因此，通過在全球不

同地區對在外核頂部傳播的 K 波之系統性觀測，並且與不同外核頂部速度模型

所得到的理論模擬結果進行比較，就可以判斷外核頂部是否有低速區域存在。

本研究中，我們首先收集震央距為120°-140°範圍内的測站之寬頻地震波形記

錄，並測量SKKS 波與 SmKS 組合波（包括 S3KS、S4KS 和 S5KS，但主要是 S3KS）

的到時，將SKKS 波與希爾伯特轉換之 S3KS 波透過交互相關函數得到兩者之間觀

測的相對到時差。然後我們使用Direct Solution Method (DSM)來模擬全球平均模型 PREM 中在相同測站的理論波形,並從理論地震圖上得到模型預測的相對到時差。

我們還可以進一步降低CMB 以下一部分深度範圍内的 P 波速度,再計算其 DSM 理

論地震圖，得到外核頂部有低速層存在的模型預測之相對到時差。由於S3KS 波在

外核頂部的K 波部分射綫長度大於 SKKS 波，外核頂部速度對 S3KS 波的影響會

比對 SKKS 影響大。因此外核頂部低速層的存在會使 S3KS-SKKS 相對到時差增

大。我們從全球地震數據中心IRIS 搜集全球分布規模大於 6.0 且深度大於 400 公

里之地震，並系統性地處理時間從1990 年至今之 78 個地震，測量得到 606 筆精

確的S3KS-SKKS 相對到時差。由我們測量之結果顯示 CMB 底下 400 公里的部分

有相對於PREM 模型的低速層存在。最後透過 Bayesian 逆推方法來尋找統計上的

期望值模型。逆推結果也顯示，在CMB 底下 550 公里深度範圍内地球外核的實際

速度低於 PREM 模型將近 0.11%，所以我們的研究結果為外核頂部有輕物質的存

在提供了強有力的證據。

關鍵字: 外核、Direct-Solution Method (DSM)、SmKS 波、Bayesian 逆推。

(4)

ABSTRACT

The solid inner core of the Earth consists of heavy minerals Fe and Ni with a fraction of light elements such as O, S and Si. These lighter elements are expelled from the inner core during its formation, rise up through the outer core as the result of buoyancy, and are trapped below the core-mantle boundary (CMB). Seismological evidence has been presented both for and against the existence of light materials at the top of the outer core. In this study, we use waveforms of recorded and modeled SmKS waves to investigate the effect of velocity perturbation under the CMB on the differential traveltimes between SKKS and S3KS waves. Due to the long propagation distance and interference with neighboring phases, the arrival times of SKKS and S3KS waves are difficult to define accurately in the records. Therefore, in our analysis we measure both the observed and model-predicted traveltimes by cross-correlating the waveform of the Hilbert-transformed S3KS with that of SKKS. We obtained 606 high-quality S3KS-SKKS differential traveltimes from 78 deep earthquakes (depth  400 km). We use synthetic seismograms calculated by the direct-solution method (DSM) in a suite of one-dimensional models with different structural profiles under the CMB to examine the existence of a zone of lowered velocity at the top of the outer core. Then we conduct a Bayesian inversion of the observed differential traveltimes for the velocity structure at the top of the outer core. The Metropolis-Hastings Monte Carlo algorithm is adopted for an efficient sampling of the model space. Inversion result indicates that the seismic velocity in the 550-km layer under the CMB is on average 0.11% lower than that in PREM. The clear depth-dependent velocity profile strongly favors the existence of light elements and chemical stratification at the top of the Earth's outer core.

(5)

LIST OF FIGURES

Fig. 1.1 Cartoon showing outer core with light materials on its top ... 2

Fig. 2.1 SKKS and S3KS ray paths in model PREM ... 6

Fig. 2.2 Distributions of earthquakes and stations ... 7

Fig. 2.3 Great-circle paths and waveforms from Sumatra earthquake to stations in North America ... 9

Fig. 2.4 Model PREM and four PREM-like models ... 12

Fig. 2.5 DSM synthetics for PREM and PREM-OC1 ... 13

Fig. 2.6 Zoom-in of DSM synthetics for PREM and PREM-OC1 ... 14

Fig. 2.7 Measuring S3KS-SKKS differential traveltime from a record ... 16

Fig. 2.8 Measuring S3KS-SKKS differential traveltime from a synthetic ... 18

Fig. 2.9 Recorded and synthetic SKKS/S3KS waveforms with quality score 19 ... 19

Fig. 3.1 Ray parameters and turning depths of SKKS and S3KS in PREM and PREM-OC1... ... 24

Fig. 3.2 S3KS-SKKS TauP differential times in PREM and PREM-OC1 ... 24

Fig. 3.3 Observed S3KS-SKKS differential traveltimes ... 25

Fig. 3.4 Comparison between observed S3KS-SKKS differential traveltimes and predictions by PREM ... 26 Fig. 3.5 Comparison between observed S3KS-SKKS differential traveltimes and

(8)

predictions by PREM-OC1. ... 26

Fig. 3.6 Comparison between observed S3KS-SKKS differential traveltimes and predictions by PREM-OC2. ... 27

Fig. 3.9 Double-differential traveltimes for PREM ... 29

Fig. 3.10 Double-differential traveltimes for PREM-OC2 ... 30

Fig. 3.13 Histograms of the double-differential traveltimes ... 33

Fig. 4.1 Difference in TauP- and DSM-predicted S3KS-SKKS differential times for PREM ... 39

Fig. 4.2 Observed S3KS-SKKS differential traveltimes used in inversion ... 40

Fig. 4.3 Models involved in Bayesian inversion ... 41

Fig. 4.4 Histogram of double-differential traveltimes for PREM-SOC ... 42

Fig. 4.5 Comparison with previous study ... 43

(9)

LIST OF TABLES

Table A.1 List of Deep Earthquakes Used in This Study ... 49 Table B.1 List of Stations Yielding Data ... 52 Table C.1 List of 606 S3KS-SKKS Differential Traveltime Measurements ... 61

(10)

Chapter 1 Introduction

To the lowest order, the seismic structure of the Earth's interior can be simply divided into several layered spherical shells, including an outer solid crust, a highly viscous upper mantle, a solid mid and lower mantle, a liquid outer core that is much less viscous than the upper mantle, and a solid inner core.

It has been widely accepted that the Earth's solid inner core consists mainly of heavy minerals Fe (~85%) and Ni (~6%) with a significant fraction of light elements such as O, S and Si (e.g. Jephcoat & Olson, 1987; Stixrude et al., 1997). As the temperature cools with time, Fe and Ni solidify and the inner core grows in radius. In the process, the light elements O, S and Si are expelled from the inner core into the outer core.

The outer core is composed of the same heavy materials Fe and Ni as the inner core but in a liquid state. There has been much less interest in studying the seismic structure of the outer core in comparison to the other regions. One reason is that the liquid outer core has only compressional wave which puts a stronger limitation on our ability to resolve its structural details. The other reason is that the outer core has always been considered to be more or less chemically homogenized by vigorous internal convection (Masters, 1979). Any variation in the structure of the outer core is thought of purely due to the depth gradient of pressure and temperature. This consideration is also the basis in establishing the outer core structure in the global average one-dimensional (1D) models such as PREM (Dziewonski & Anderson, 1981).

However, the light elements O, S and Si expelled from the inner core may alter this situation. Because of gravitational buoyancy, these lighter materials rise up through the Fe/Ni filled outer core and become trapped beneath the core-mantle boundary (CMB), a

(11)

solid thermal-chemical barrier, as illustrated in Fig.1.1. Such a scenario implies that there are light materials (in comparison to the expected heavy elements Fe and Ni), and therefore a chemical stratification may exist at the top of the outer core (e.g. Lister &

Buffett, 1998; Buffett & Seagle, 2010; Frost et al, 2010), but the validity of this scenario needs to be confirmed by evidence gathered from geophysical (especially seismology) observations, laboratory experiments, and numerical modelings.

Fig. 1.1 Cartoon illustrating the schematic model of an outer core with light materials on its top. The light blue curvy arc indicates the CMB with possible topography. The inner core consists mainly heavy elements Fe and Ni, and a fraction of light materials O, S and Si (red dots), which are expelled from the inner core during its formation, rise up through the outer core, and are eventually trapped beneath the CMB.

The seismic waves propagating through the outer core, such as SKS and its core multiples SmKS, provide us the direct means to estimate the wave speed in the outer core and to infer the structure there. The effect of lighter materials on wave speed seems to be straightforward. In the liquid outer core, the speed of the compressional K wave is

_ ^_,

  (1.1) where  and  are the bulk modulus and density, respectively. Therefore, intuitively a density decrease associated with the existence of lighter elements would result in an increase in wave speed. However, Helffrich (2012) demonstrated that the elastic

(12)

property of a multi-component liquid behaves non-linearly if the mixture has not reached an ideal state, and may offset the effect of density change on the wave speed.

Therefore, chemical stratification at the top of the outer core due to the existence of lighter materials there could lead to the non-intuitive outcome of wave speed decrease.

The rapid increase in global deployments of permanent and temporary networks has provided more and more opportunities of obtaining high-quality records of SmKS waveforms which make the detailed investigation of the outer core feasible. Various seismic observations have been reported. Based on a limited dataset using earthquakes in Indonesia and stations in Africa, Tanaka (2004) reported possible low P-wave speed in a 50-km layer below the CMB. Eaton & Kendall (2006) found that a 12-km thick high-speed/low-density layer below the CMB best explained their SmKS traveltimes obtained from broadband seismograms at stations in Canada from a few earthquakes in Indonesia and Fiji. Tanaka (2007) analyzed a series of global records of SmKS, and proposed a model with a 1.4% lower speed than PREM in the top of 90 km of the outer core. More recently, a pair of studies by two independent groups reported very different findings. Alexandrakis & Eaton (2010) carefully analyzed SmKS traveltimes from 44 teleseismic earthquakes with a variety of source depths and in epicentral distance range 124°-140°, and concluded that there is no evidence for outer core stratification. In contrast, Helffrich & Kaneshima (2010) modeled the SmKS traveltimes along the two paths from South America to Japan and from Fiji to Europe, and suggested that the wave speed in the 300-km layer below the CMB can be up to 0.3% lower than PREM.

The validity of the previous studies on the outer core structure using the SmKS traveltimes may be limited by three aspects. First, all these authors have used relatively few earthquakes and stations, and their data are limited in both quantity and geographical coverage. Secondly, most of them have used the traveltimes of SmKS up

(13)

to S5KS, which are often weak and have low signal-to-noise ratios, which reduced the measurement quality. Finally, the long traveling distances of the SmKS waves lead to attenuation of high-frequency contents and in the recorded signals. Therefore, these waves are usually best observed at frequencies below 0.2 Hz, which makes measurement of the onset times of arrivals very difficult. However, all the previous studies have been using ray theory to model the SmKS traveltimes, but the bias introduced by finite-frequency effect are not addressed.

The purpose of this study is to conduct a careful and thorough investigation of the seismic structure at the top of the outer core by establishing a larger and higher-quality dataset of SmKS traveltimes from globally distributed earthquakes and stations, and a more accurate modeling of these traveltimes. We focus on the differential traveltimes between SKKS and S3KS waves to both ensure the quality of the measurements and eliminate the possible contributions from factors such as source location errors and lateral structural variations along the paths in the mantle. We measure the differential traveltimes by the cross correlation of SKKS and (Hilbert-transformed) S3KS waveforms, and we model the finite-frequency differential traveltimes by the cross correlation of waveforms calculated by the direct-solution method (DSM).

The outline of this thesis is as follows. In Chapter 2, we describe the process for measuring the S3KS-SKKS differential traveltimes, and the DSM method used in modeling these data. Then we present the observed S3KS-SKKS differential traveltimes as well as predictions by PREM and several PREM-like models in Chapter 3. In Chapter 4, we carry out inversion method to search more models to understand the velocity structure below CMB 550 km. Finally we discuss and summarize our results in Chapter 5.

(14)

Chapter 2 Data and Methodology

2.1 Seismic phases for studying the top of the outer core

Similar to the investigation of the other parts of the Earth's interior, the best means to probe the properties of the outer core is by analyzing the seismic waves that passing through it, namely the phases with a K-wave (compressional wave in the liquid outer core) segment, such as the SKS wave. The seismic properties of the outer core revealed by these waves can be combined with thermal-dynamical laws derived from laboratory experiments and/or theoretical studies to infer the chemical compositions and physical conditions of the outer core.

The most direct and robust observation on the seismic property at the top of the outer core can be obtained from the difference between the travel times of SKKS wave and the group arrival of SmKS waves, with m being 3, 4 and 5 (i.e. S3KS, S4KS and S5KS, respectively). The best source-receiver geometry to clearly observe the difference between SKKS and SmKS phases is in the distance range 120º-140º and from a deep earthquake (400 km and deeper) to avoid the interference of near-source surface reflections. The S3KS phase dominates S4KS and S5KS in this distance range, S3KS is weaker but beyond 135º, and the waveform of the SmKS group becomes more complicated and difficult to analyze.

At the distances of 120º, 130º and 140º, SKKS wave turns near the depths of 300 km, 400 km, and 500 km, respectively, below the CMB. The first Fresnel zone width of the SKS wave at this distance range is larger than 500 km at the CMB (Zhao & Chevrot, 2011), where the first Fresnel zones of SKKS and S3KS waves are similar in size. On the other hand, the distance between the CMB-piercing points of SKKS and S3KS is

(15)

about 250 km. Therefore, as shown in Fig.2.1, the ray paths of SKKS and S3KS in the distance range 120º-140º are also almost completely overlap outside the core, and any anomaly in their differential travel time TS3KS-TSKKS is almost entirely due to the perturbations of the K-wave speed in the top ~400 km of the outer core (between dashed and blue semi-circles in Fig.2.1). This makes the differential time TS3KS-TSKKS a very useful observable in studying the top of the outer core.

Fig. 2.1 SKKS and S3KS ray paths in model PREM at an epicentral distance of 120°

from a deep earthquake with a source depth of ~600 km. The dashed semi-circle is ~400 km below the CMB (light blue semi-circle). In the epicentral distance range 120º-140º, ray paths of SKKS and S3KS almost coincide outside the CMB. As a result, anomalies in the differential traveltime TS3KS-TSKKS are mostly due to velocity perturbations in the

~400-km zone below the CMB.

2.2 Deep earthquakes and records used in this study

In this study, we measure the differential travel times of SKKS and S3KS using records from globally distributed deep earthquakes and seismic stations. Our goal is to have as many measurements as possible, in order to increase the redundancy to reduce data uncertainty and expand the geographical coverage. For this purpose, we

(16)

downloaded from IRIS (http://www.iris.edu/wilber3/) the two horizontal-component broadband records at stations in the distance range 120º-140º from earthquakes between 1998 and 2013. The earthquakes have depths greater than 400 km and moment magnitudes 5.7 ≤ MW ≤ 8.3. Fig.2.2 displays the distributions of earthquakes and stations used in this study. Detailed information on the sources and stations are provided in Tables A.1 and B.1 in the Appendix A and B, respectively.

Fig. 2.2 Distributions of the 78 earthquakes (1998-2013) and the stations used in this study. Circles represent earthquakes, with their sizes and colors indicating magnitudes and focal depths, respectively. Red triangles are 744 stations providing waveforms.

Gray circles show 6 events analyzed with no usable data. Thin yellow lines are source-receiver great-circle paths with usable data.

We only use the radial component here because all the relevant waves (SKKS and S3KS) arrive at the station as an S wave, with motion perpendicular to the ray path. At 120º-140º degrees distance, the SKKS and S3KS waves arrive at the station nearly vertically (see Fig.2.1), therefore the amplitudes of SKKS and S3KS will be

(17)

mostly horizontal. Furthermore, during the transmission at the CMB, the compressional K wave in the outer core only converts to P wave and the longitudinally-polarized SV wave in the mantle. Therefore, the SKKS and S3KS waves arriving at the station are both polarized in the radial component.

An example of the radial-component seismograms is presented in Fig.2.3 for the 25 July 2004 deep earthquake in Sumatra recorded at 54 stations in North America. Due to the rather long traveling distance, higher frequency content of the SKKS wave and the SmKS group is strongly attenuated. As a result, at these long distances these phases have relatively low signal-to-noise ratio (SNR), and it is necessary to remove the high-frequency content of the seismograms to ensure sufficient SNR. Through trial-and-error, we found that the best frequency band to have clear waveforms of both SKKS wave and the SmKS group is 0.005-0.2 Hz. The waveforms in Fig.2.3 have all been deconvolved with the instrument responses to velocity, and bandpass filtered.

There are many publicly available programs, such as the TauP Toolkit (Crotwell et al., 1999), which can be used to calculate the travel times of SKKS and SmKS waves based on 1D models (e.g. PREM) easily and quickly. However, these ray-theory travel times are model predictions of the onset times of the corresponding waves at high frequency, which cannot be easily located in the observed seismograms in the distance range 120º-140º used in this study. There are two main reasons for this. One is the fact that we need to filter the recorded seismograms to below 0.2 Hz to achieve sufficient SNR.

Another reason is that the interference of nearby phases immediately before the arrivals of SKKS and SmKS waves often distorts the beginning parts of these waves and makes it impossible to identify their onset times. As a result, the most accurate and reliable approach to measuring the differential travel time TS3KS-TSKKS is through the cross correlation of the recorded SKKS and SmKS waveforms.

(18)

Fig. 2.3 (Top) Focal mechanism (beachball) of the 25 July 2004 Sumatra earthquake (Mw=7.3, depth=600.5 km) and the great-circle paths (red lines) to 54 stations (green triangles) in North America in the distance range 120º-140º. Color dots on the paths are locations where SKKS and S3KS interact with the CMB. (Bottom) Radial-component velocity records from Sumatra earthquake to the 54 stations. Red line is the predicted SKKS traveltime curve based on PREM predicted by the TauP Toolkit (Crotwell et al., 1999), whereas the blue lines are predicted traveltime curves of SmKS based on PREM.

All waveforms have been bandpass filtered between 0.005 Hz and 0.2 Hz.

(19)

Since the SKKS and SmKS differential travel times have to be measured from waveforms, we also need synthetic seismograms of these phases in order to have equivalent model-predicted differential travel times to compare with. In this study, we adopt the numerical algorithm of the direct-solution method (DSM), which provides the exact synthetic seismograms based on 1D models.

2.3 Waveform simulation: the direct-solution method

As mentioned before, the SKKS and S3KS waveforms are best observed at frequencies below 0.2 Hz in records. As a result, their onset times are difficult to pick and have relatively large uncertainties. In order to ensure data quality, we measure the SKKS and S3KS differential traveltimes by cross correlations of their waveforms in records. This requires us to be able to calculate their synthetic seismograms so that we can obtain equivalent measurements to compare. In this study, we use the direct-solution method to compute exact synthetic seismograms based on 1D models (Cummins et al., 1994; Geller & Ohminato, 1994; Kawaii et al., 2006).

The direct solution method (DSM) is closely related to the normal-mode summation method. Both methods solve the seismic wave equation by decomposing the solutions into spherical harmonics horizontally and radial eigenfunctions vertically. However, in normal-mode summation, the synthetics are calculated in two steps: (1) compute the radial eigenfunctions (solving the homogeneous equation with no source); and (2) sum the eigenfunctions with weights determined by the earthquake source to get the synthetics. Eigenfunctions calculated in Step (1) can be archived and used for any earthquakes. This makes the normal-mode summation method more efficient in computing synthetics for many earthquakes.

(20)

The DSM does not solve for source-free radial eigenfunctions, but finds the radial solutions directly for the inhomogeneous wave equation with an earthquake source. It uses triangular splines as radial basis functions, and the solution is specified by the expansion coefficients. The solution is expressed in the frequency domain:

( , , , ) [ ( ) ( ) ( , ) ( ) ( ) ( , )

( ) ( )( ) ( , )],

klm k lm klm k lm

k l m

klm k lm

r A X r Y B X r Y

C X r Y

        

  



 

 



u r

r





(2.1)

where Ylm(θ,ϕ) is the surface spherical harmonic of degree l and order m, and Xk(r) is the k-th radial basis functions of order. The expansion coefficients Aklm, Bklm and Cklm are constrained by the earthquake source and solved by systems of linear equations. Once the frequency-domain solution is obtained, time domain response (waveform) is obtained by an inverse Fourier transform. The steps for calculating DSM synthetics are as follows:

DSM modeling procedure

1. Edit the input file including the Earth model, the earthquake source location and moment tensor; and the station locations (there can be more than one stations at a time). The file also contains parameters for the maximum frequency of the seismograms as well as their length in time.

2. Run the DSM code to compute the frequency-domain responses of the P-SV and SH systems. Since the frequency-domain samples can be divided and calculated independently, the computation can be done in parallel on any number of processors to reduce the CPU time. The P-SV and SH systems calculated on different processors are merged after the parallel computations are completed.

3. Apply an inverse Fourier transform to the merged frequency-domain response to obtain the time-domain seismograms. The DSM results are velocity response.

(21)

4. The results are then convolved with a Gaussian-shaped source-time function given by the time shift and half duration in the earthquake's GCMT solution.

Fig. 2.4 Five models used to study the effect of the structure atop the outer core on differential traveltime TS3KS-TSKKS. PREM is shown in black. Gray, green, blue and red lines show models PREM-OC1, PREM-OC2, PREM-OC3 and PREM-OC4 which are the same as PREM except in the 400-km layer immediately below CMB, where the wave speed is reduced by about 0.23%, 0.11%, 0.06% and 0.09% on average.

(22)

Fig. 2.5 (Top) DSM radial-component synthetics based on PREM for the 2004 Sumatra earthquake at 54 stations in North America (source and station locations are shown in Fig.2.3). Red and blue lines are traveltime curves of SKKS and SmKS waves calculated by the TauP Toolkit. (Bottom) Same as the top plot but for Model PREM-OC1. The traveltime curves are for model PREM. All waveforms have been bandpass filtered between 0.005 Hz and 0.2 Hz.

(23)

Fig. 2.6 Zoom-in of the waveforms in Fig. 2.5 around the SKKS and S3KS arrivals. In the bottom plot, delays can be seen in the SKKS and S3KS waveforms relative to the PREM traveltime curves for corresponding waves.

Follow this procedure, it usually takes about 5 hours to compute the DSM synthetics up to 1 Hz using 8 computing nodes (12 cores on each node). The radial-component

(24)

synthetics are then filtered in the same way as the records to make S3KS-SKKS differential traveltime measurements. Using DSM it is very easy to replace the Earth model or the earthquake source and calculate their DSM synthetics. Fig.2.4 shows the model PREM as well as several PREM-like models we use in this study. Fig.2.5 displays the DSM synthetic seismograms for the same earthquake and stations as those in Fig.2.3 but for two different Earth models: PREM and PREM-OC1. PREM-OC1 is the same as PREM everywhere except for the top 400 km of the outer core, where the P-wave speed is reduced by ~0.23% on average. The synthetics clearly show that the S3KS-SKKS differential traveltimes are larger based on PREM-OC1, a model with a lower speed at the top of the outer core.

2.4 Differential traveltime measurements from records

Here we describe the process for measuring the differential travel times between the SKKS and the S3KS waves. Fig.2.3 shows the recorded radial-component seismograms from the 2004 Sumatra earthquake at 54 stations in North America in the distance range 120º-140º. We take the station FCC as an example to illustrate the measurement process, as shown in Fig.2.6.

First we choose an appropriate time window for the SKKS wave based on its TauP-predicted travel time in PREM. Usually the waveform of SKKS can be clearly identified, and we choose its window with a width of only half of its period in order to minimize the contamination of neighboring phases (Fig.2.6a). Then we apply a Hilbert transform to the recorded waveform with a focus on the S3KS wave:

S3KS

S3KS S3KS S3KS

1 ( )

( ) ( ) ( ) ( ) ( ) u ,

Hu t h t u t u h t d d

t

     

 

 

    

 

 (2.2)

where h(t) stands for the Heaviside-step function. The definition in eq.(2.2) shows that

(25)

the Hilbert transform adds a π/2 phase shift to the S3KS waveform. This is to account for the extra caustic along the S3KS ray path in comparison to SKKS. A time window is also determined in the Hilbert-transformed S3KS waveform (Fig.2.6b).

Finally we compute the cross-correlation between the SKKS waveform and the Hilbert-transformed S3KS waveform

SKKS S3KS SKKS S3KS

(u Hu )( )t ^ u ( ) Hu (t  ) ,d

 



  (2.3) and the differential time TS3KS-TSKKS is the lag time when the cross correlation reaches the maximum (Fig.2.6c). In this study, time series operations, such as deconvolution of instrument responses, filtering, Hilbert transform, and cross correlation are done in the software package SAC (Seismic Analysis Code, Helffrich et al, 2013).

Fig. 2.7 (a) Radial-component record of the 2004 Sumatra event to station FCC (yellow triangle in Fig. 2.3) showing the SKKS and S3KS. Red bars indicate the time window for SKKS used in making cross-correlation measurement; (b) Hilbert transform of the record in (a). Blue bars indicate the time window for Hilbert-transformed S3KS; and (c) Cross-corrlation of the windowed portions of the waveforms in (a) and (b). The dashed line shows the differential traveltime measurement. All waveforms have been band-pass filtered between 0.005 Hz and 0.2 Hz.

(a) (b)

(c)

(26)

This results in the observed differential traveltime between the SKKS and S3KS waves obtained from records:

 T T_S3KST_SKKS. (2.4) Following this process, we analyzed the waveforms of a total of 78 earthquakes and obtained high-quality measurements of differential traveltimes TS3KS-TSKKS at 374 stations.

2.5 Differential traveltime measurements from synthetics

The radial-component DSM synthetics in Fig.2.5 in the two models (PREM and PREM-OC1) from the Sumatra earthquake at 54 stations in North America in the distance range 120º-140º. We use the same station FCC to illustrate the measurements of differential traveltime TS3KS-TSKKS from synthetics, as shown in Fig.2.7. The process is exactly the same as for the records: A time window is chosen for the SKKS wave (Fig.2.7a). Then a Hilbert transform is applied to the waveform and a time window is chosen for the transformed S3KS waveform (Fig.2.7b). Finally the waveforms in the two windows are cross correlated, and the differential traveltime is determined from the maximum of the cross correlation (Fig.2.7c).

This results in the model-predicted differential traveltime between the SKKS and S3KS waves obtained from synthetic seismograms:

MODEL MODEL MODEL

S3KS SKKS .

T T T

     (2.5) Note that we use the same symbol in eq.(2.5) as in eq.(2.4), but with a "~" overhead to indicate model-predicted quantities. For every observed differential traveltime, we calculate the corresponding DSM synthetics for PREM and for several PREM-like models (Fig. 2.4), and measure SKKS and S3KS differential times in the same way.

(27)

Fig. 2.8 (a) Radial-component synthetic of the Sumatra event at station FCC showing on the SKKS and S3KS waves. Red bars indicate the time window for SKKS used in making cross-correlation measurement; (b) Hilbert transform of the synthetic in (a); and (c) Cross-corrlation of the windowed portions of the waveforms.

2.6 Quality control of differential traveltime measurements

The examples in Figs.2.6 and 2.7 demonstrated very robust measurements of differential traveltimes between SKKS and S3KS waves, as indicated by the similarity in the waveforms of SKKS and the Hilbert-transformed S3KS. However, as the epicentral distance increases, the S3KS wave becomes weaker, and the SNR decreases.

In addition, radiation pattern can also reduce the strength of SKKS and S3KS waves.

The effect of lateral heterogeneity introduces further complexities in the recorded waveforms. Therefore, the quality of the SKKS and S3KS waveforms are not always perfect, and there is a need to impose a quality control on both the observed and

(a) (b)

(c)

(28)

model-predicted differential traveltime measurements. In this study, we follow these steps and criteria:

1. For each data, we determine its measurement quality based on both the recorded and synthetic seismograms. In other words, we require the waveform of the Hilbert-transformed S3KS to be similar to SKKS to be true for both record and its corresponding synthetic in PREM. The requirement for synthetic ensures that the effect of lateral heterogeneity is small. Based on the visual inspection of the recorded waveforms of SKKS and S3KS waves, a quality score of 1 to 10 is assigned, with 10 being the highest quality. Similarly, a quality score is assigned to the synthetic. As a result, the highest quality differential traveltime has a score of 20, as in the case of FCC shown in Figs.2.6 and 2.7.

Fig. 2.9 (Top panels) Record (left) at station YBH from the Sumatra event and its Hilbert transform (right). A quality score of 10 is given to this record. (Bottom panels) Synthetic (left) and its Hilbert transform (right) for the same source-station pair. A quality score of 9 is given due to slight distortions of S3KS waveform.

(29)

2. In many cases, the SKKS waveforms are very clear, but the waveforms of S3KS are affected by other phases arriving slightly earlier or later, thus distorting the beginning and/or ending parts of the S3KS waveforms. This can also happen in records and/or synthetics. We assign quality scores of 19 and 18 based on the level of distortion of the S3KS waveforms. As shown in Figs.2.8 and 2.9.

Fig. 2.10 (Top panels) Record (left) at station TTW from the Sumatra event and its Hilbert transform (right). (Bottom panels) Synthetic (left) and its Hilbert transform (right) for the same source-station pair. Both record and synthetic are given a quality score of 9 due to slight distortions in the S3KS waveforms.

3. Sometimes there appears to be some anomaly in the relative amplitudes between the SKKS and S3KS waves, such as the record at station PNT from the Sumatra earthquake shown in Fig.2.10. This will also affect the quality score given to the measurement.

(30)

Fig. 2.11 Quality scoring for station PNT from the Sumatra event. The record is given a quality score of 9 due to anomalously large amplitude ratio between SKKS and S3KS, whereas the synthetic is given a score of 8.

Fig. 2.12 Quality scoring for station KNGO from the Sumatra event. Both record and synthetic are given a score of 8.

(31)

4. As the epicentral distance increases, especially beyond 135º, the relative strength of S3KS in the SmKS group decreases, and the S3KS waveforms have increasingly larger interference with the neighboring phases. As a result, the S3KS waveforms become less clear and the differential traveltimes obtained by waveform cross correlation are less robust, as shown in Figs.2.11 and 2.12. In this study, only measurements with quality scores of 17 and above are used for further analysis.

Fig. 2.13 Quality scoring for station LMQ from the Sumatra event. The record is given a quality score of 7, whereas the synthetic is given a score of 8.

(32)

Chapter 3 Modeling SmKS Differential Traveltimes

Follow processes described in Chapter 2, we measured the differential traveltimes between SKKS and S3KS from both recorded seismograms and the DSM synthetics in PREM as well as several PREM-like 1D models. We also assigned quality score to each observed differential traveltime measurement based on visual inspection of the record and its corresponding synthetic in PREM. In the end, a total of 606 measurements of quality score 17 and above are retained for further analysis. In this chapter, we compare the observed differential traveltimes with their corresponding predictions by PREM and PREM-like models to examine their discrepancies and their possible reasons.

3.1 Ray theory features of SKKS and S3KS waves

In order to identify the region of the Earth whose structure has the most influence on the S3KS-SKKS differential traveltimes, we first look at some ray theory features of these waves such as their ray parameters and depths of turning points in the outer core.

Fig.3.1 displays the values of the ray parameters and the turning depths of SKKS and S3KS waves in the distance range 120º-140º, calculated for Models PREM and PREM-OC1 by the TauP Toolkit. The ray parameters of the SKKS and S3KS waves in this distance range vary in 5.8~6.5 s/deg and 6.8~7.1 s/deg, respectively. For both waves, the ray parameters monotonically decrease with epicentral distance. The S3KS waves all turn in the top 200 km of the outer core, where the SKKS waves can reach down to

~500 km below the CMB. The S3KS-SKKS differential traveltimes for these two models displayed in Fig.3.2 clearly indicate that the model with a slower top outer core (PREM-OC1) has larger differential traveltimes at all distances.

(33)

Fig. 3.1 (Top) Ray parameter distributions for the SKKS and S3KS waves in epicentral distance 120°-140°, calculated for Models PREM (pink) and PREM-OC1 (gray) by the TauP Toolkit. Ray parameters of SKKS and S3KS decrease monotonically with distance almost linearly. (Bottom) Depths of the turning points of SKKS and S3KS measured from the CMB in the epicentral distance range 120º-140º, calculated for Models PREM (pink) and PREM-OC1 (gray) by the TauP Toolkit. For both waves, their turning depths decrease monotonically with ray parameter (or increase with distance). SKKS turns in the depth range 300-500 km; S3KS turns in the range 100-200 km below the CMB.

Fig. 3.2 S3KS-SKKS differential traveltimes in the epicentral distance range 120º-140º, calculated for Models PREM (pink) and PREM-OC1 (gray) by the TauP Toolkit.

(34)

3.2 Observed and model-predicted differential traveltimes

Here we discuss the observed differential traveltimes ΔT=TS3KS-TSKKS obtained from recorded seismograms. After a rigorous inspection of all the recorded and synthetic waveforms of SKKS and S3KS waves, we selected 606 high-quality (quality score of 17 and above) differential traveltime measurements for further analysis. These measurements come from 78 events (5.8 ≤ Mw ≤ 7.7) and they are plotted against the epicentral distance in Fig. 3.3.

Fig. 3.3 606 observed S3KS-SKKS differential traveltimes ΔT=TS3KS-TSKKS with quality scores of at least 17 plotted against the epicentral distance.

All the data in Fig.3.3 have been selected for their relatively high quality scores, and most of them are indeed clustered around a linear trend showing a clear increase in the differential traveltimes with the epicentral distance. However, the scattering of the measurements suggests a significant uncertainty in the data, and a number of them are quite far from the main cluster. Therefore, instead of using these raw measurements for analysis, we average them in 1º bins so that the main features in the data can be seen more clearly, and the comparisons with model-predicted results can also be done, as shown in Figs.3.4-3.8.

(35)

Fig. 3.4 Comparison between observed S3KS-SKKS differential traveltimes and predictions by Model PREM. The 606 measurements in Fig.3.3 have been averaged every 0.5º over 1º bins, and plotted in black dots with bars showing the standard deviation. Purple triangles are predictions obtained from DSM synthetics for PREM and averaged in the same way. Model predictions are clearly smaller than observations, suggesting that PREM is too fast beneath the CMB.

Fig. 3.5 Same as Fig.3.4 but the predictions (gray symbols) are for Model PREM-OC1.

With an average perturbation of about -0.23% below the CMB, Model PREM-OC1 is clearly too slow in the top of the outer core and will not be shown in later discussion.

(36)

Fig. 3.6 Same as Fig.3.4 but the predictions (green symbols) are for Model PREM-OC2.

With an average perturbation of about -0.11% below the CMB, Model PREM-OC2 performs better than PREM and PREM-OC1. Some of the model predictions are a bit larger than observations indicating that PREM-OC2 may be slightly too slow.

Fig. 3.7 Same as Fig.3.4 but the predictions (blue symbols) are for Model PREM-OC3.

With an average perturbation of about -0.06% below the CMB, Model PREM-OC3 seems to be similar to PREM. Many model predictions are still smaller than observations.

(37)

Fig. 3.8 Same as Fig.3.4 but the predictions (red symbols) are for Model PREM-OC4.

With an average perturbation of about -0.09% below the CMB, predictions by Model PREM-OC4 more closely follow the observations.

3.3 Double-differential travel times

The comparisons in Figs.3.4-3.8 between the observed S3KS-SKKS differential traveltimes with the predictions by several PREM-like models suggest that Model PREM is indeed too fast in the top of the outer core. The tested models PREM-OC1 is too slow below the CMB, while PREM-OC3 may still be too fast. Models PREM-OC2 and PREM-OC4 seem to perform better, with predictions by PREM-OC4 showing closest agreement with the observed differential traveltimes.

Another way to make a more quantitative assessment of the performance of different models is by examining the double-differential traveltimes

2TMODEL T TMODEL,

    ^ (3.1) where ΔT and T^^MODEL are the observed and model-predicted differential traveltimes defined in eqs.(2.4) and (2.5), respectively. An ideal model would make perfect

(38)

differential traveltime predictions so that ²T^MODEL 0. In Figs.3.9-3.12 we plot the values of the double-differential traveltimes for the tested models. They are plotted both as against the epicentral distances as in Figs.3.4-3.8 and in mapview.

Fig. 3.9 (Top) Plot of double differential traveltimes ²T^PREM against epicentral distance. The data have been plotted by the values averaged every 0.5º over 1º-bins with the bars showing the standard deviation. (Bottom) 606 source-receiver great-circle paths for which we obtained the observed differential traveltime data. The paths are plotted such that portions outside the core are in gray while portions inside the core are in color.

The color of a path indicates the value of ²T^PREM. Positive values (warm colors) mean PREM is too fast in the outer core.

(39)

In the mapview plots in Figs.3.9-3.12, each value is plotted along the great-circle path from which the observed differential traveltime ΔT is obtained. For each great-circle, the portion of the path outside the core is plotted in gray while the portion inside the outer core is in color according to the value of ²T^MODEL. Positive values are plotted with warm colors, indicating that the model is too fast in the outer core along those paths; negative values are plotted in cool colors, indicating that the model is too slow below the CMB along those paths. Therefore, the colors and their strengths can tell us how well a model performs in predicting the S3KS-SKKS differential traveltimes.

Fig. 3.10 Same as Fig.3.9 but for Model PREM-OC2 ²T^PREM-OC2.

(40)

Due to the location of deep earthquakes in subduction zones and the distribution of seismic stations on continents, the S3KS-SKKS paths are mostly concentrated in the Pacific Ocean, Northern Atlantic and Southern Indian Oceans, and under Eurasia.

However, the dominating warm-colored paths in Fig.3.9 strongly suggest that Model PREM is on average too fast in the top of the outer core. The same is true in the paths in Fig.3.11 for Model PREM-OC3. The paths for Models PREM-OC2 (Fig.3.10) and PREM-OC4 (Fig.3.12) are less dominated by warm colors, and the paths also have

(41)

weaker colors, suggesting a relatively better performance by PREM-OC2 and PREM-OC4 in predicting the S3KS-SKKS differential traveltimes.

3.4 Quantitative assessment of tested models

The double-differential traveltimes ²T^MODEL displayed in Fig. 3.9-3.12 suggest that Models PREM-OC2 (Fig.3.10) and PREM-OC4 (Fig.3.12) provide predictions of the S3KS-SKKS differential traveltimes closer to the observed values than PREM and

(42)

PREM-OC3. Here we use the histograms in Fig.3.13 to further quantify these observations. Values of the mean and standard deviation are given in each histogram. A positive mean value indicates the model is on average too fast relative to observation, and vice versa.

Fig. 3.13 Histograms of the double-differential traveltimes ²T^MODEL for the four models. (a) PREM, (b) PREM-OC3, (c) PREM-OC4, and (d) PREM-OC2. The dashed line in each plot shows the mean value. The values of the means and standard deviations are given in each plot. Also given in Panels (b), (c) and (d) are the average perturbations of the wave speed in the top of the outer core of the different models relative to PREM.

The histograms in Fig.3.13 quantitatively confirm our observations: Model PREM is too fast. It has the largest positive mean value of 0.36 s. For Model PREM-OC3 (Fig.3.13b), the reduction of 0.06% from PREM in the speed in the top of the outer core is too small, resulting in a still relatively large positive mean value of 0.13 s. On the

(43)

other hand, Models PREM-OC3 and PREM-OC4, with average wave speed reductions of 0.09% and 0.11%, respectively, both have small mean values. In particular, Model PREM-OC4 has the smallest mean value of -0.07 s. Therefore, based on our dataset, we conclude that the wave speed in the ~400 km zone beneath the CMB is on average about 0.1% lower than that in PREM. Although the average velocity decrease relatively to PREM seems to be small, it is clearly significant in the comparison of the observed S3KS-SKKS differential traveltimes with predictions by the models we have tested. In the next chapter, we will try to quantify the optimal wave speed variation in the top of the outer core by a Bayesian inversion of the observed S3KS-SKKS differential traveltimes.

(44)

Chapter 4 Inversion of SmKS Differential Traveltimes

4.1 Bayesian inversion

In the modeling results of the observed S3KS-SKKS differential traveltimes in Chapter 3, we have seen that the structure at the top region of the outer core is on average slower than PREM by about 0.1%. This suggests the possibility chemical stratification below the CMB. Our next goal is to conduct an inversion of the 606 observed S3KS-SKKS differential traveltimes between SKKS and S3KS to obtain the details of the 1D profile of the wave speed at the top of the outer core.

In our inversion problem, the data are the 606 differential traveltimes, and the unknowns are the wave speed in the uppermost layer of the outer core parameterized as wave speeds at a number of depth samples below the CMB. We can express this inverse problem as

( ) ,

g m d (4.1) where d is the data vector of length N (N=606 here), and m is the unknown model parameter vector of length M. The functions gi (i=1, 2, ..., N) that relate the data d to the model m is almost always non-linear (such as in the case of this study), and the inverse problem (4.1) is commonly solved either using a non-linear optimization scheme or by iteratively solving its linearized version. Here we adopt the former choice and solve (4.1) by a non-linear grid search method.

Theoretically speaking, a grid search can be easy to implement. There is no need to calculate the partial derivatives of the data with respect to model parameters, nor is there a need to compute inverse matrices. We only need to be able to carry out the forward modeling given in eq.(4.1) for any given model m, and drive the grid search

(45)

according to a pre-defined objective function such as



 d g m ( ) , (4.2)

where  represents a specific measure of the misfit between the data d and its model prediction g(m), such as the L2 norm:

2 2

1

[ ( )] .

N

L i i

i

d g







 ^m (4.3)

Under this definition, the inverse problem would be solved once the grid search has found the model m which yields the minimum value of the misfit _L₂. However, in practice the model space has a high dimension (M is large), the functions gi(m) is non-linear, and the data di have uncertainty. Due to these complications, there is no basis to regard the model m with the minimum misfit as the true solution to the inverse problem. We have to replace the deterministic view of the inverse problem with a probabilistic one, and here we adopt the Bayesian approach to the inversion of the S3KS-SKKS differential traveltime data.

The basic concept of the Bayesian inversion is to make inferences about the properties of the model parameters based on the a posteriori probability distribution (e.g.

Tarantola, 1987),

( | ) ( ) ( | ),L

 m d   m d m (4.4) where (m) is the a priori probability distribution of the model m, L d m( | ) is the likelihood function indicating the probability of predicting data d by the model m, and is defined as

2 2

1

| exp[ [ ( )] 2 ].

( ) ^N _i _i _i

i

d g

L 



 





d m m (4.5)

The a posteriori probability distribution ( | )m d represents the probability of the

(46)

solution of the inverse problem being the model m under the constrains of the observed data d with standard deviation τ and the specified prior information. If we can make an exhaustive sampling of the entire model space, we will have the complete knowledge of

( | )

 m d , which can be used to quantify the features of the model parameters that we are interested in, such as the model expectation and standard deviation. However, a complete sampling of the model space is an impossible task because of the extremely heavy computational cost.

4.2 Metropolis-Hastings Monte Carlo algorithm

As discussed in the previous section, the Bayesian inversion depends on making model inferences based on sufficient knowledge of the a posteriori probability distribution ( | )m d which is too costly if the entire model space must be sampled.

Therefore, we need an efficient way to sample the models so that we can acquire necessary knowledge of ( | )m d without the need for an exhaustive search in the entire model space. This is achieved by the so-called Metropolis-Hastings Monte Carlo (MHMC) algorithm (e.g. Metropolis et al., 1953; Hastings, 1970; Mosegaard &

Tarantola, 1995). The MHMC is based on a random (i.e. Monte Carlo) sampling of the model space. The Metropolis algorithm (Metropolis et al., 1953), later modified by Hastings (Hastings, 1970), provides a mechanism of rejection and acceptance of the sampling models. When the number of sampled models is large enough, the density distribution of the accepted models resembles ( | )m d , i.e. the MHMC samples the of the model space a posteriori probability distribution of the model space (Mosegaard &

Tarantola, 1995). Therefore, the MHMC has a natural connection with the Bayesian inversion and its efficiency in model sampling has made it a popular choice in solving

(47)

non-linear inverse problems.

The MHMC algorithm involves the following steps:

[1] Start with a so-called burn-in stage in which a series of random models are drawn from the model space. At the end of the burn-in stage, the current model sample is denoted as mi.

[2] Calculate L(d | mi), the likelihood function of model mi as defined in eq.(4.5).

[3] Draw a new random model mj from the model space.

[4] Calculate L(d | mj).

[5] Calculate the ratio between the values of the likelihood function of the two consecutive model samples

( | ) ( | ).

j i

r L

 L d m

d m (4.6) This ratio is used as the probability by which the new model mj is accepted or rejected, i.e., if r  1, then model mj is accepted with mi+1=mj; if r  1, then a random number s between 0 and 1 is generated, and if s ≤ r, mj is accepted with mi+1=mj, otherwise mj is rejected and mi+1=mi.

As Steps [2] to [5] are repeated and the number of accepted models increases, model expectations can be calculated frequently. Once the model expectation becomes stable, the sampling can be stopped.

4.3 Bayesian inversion of SmKS differential traveltimes

Although the MHMC algorithm provides an efficient scheme to conduct the Bayesian inversion, the number of model samples is still large. Since predicted data are needed to calculate the likelihood function, therefore for each sampled model we have to run the DSM simulations for all earthquakes and measure the S3KS-SKKS

(48)

differential traveltimes, which is still impractical for hundreds or thousands of models.

Therefore, we use the much more efficient Taup Toolkit (Crotwell et al., 1999) to predict the differential traveltimes. For a given model, there is a slight difference between the differential traveltimes calculated by the TauP Toolkit and those measured from DSM-calculated waveforms. The TauP result assumes high frequency and only SKKS and S3KS waves, whereas the DSM result is from full waveform which includes the finite-frequency effect as well as possible interference with neighboring phases near SKKS and S3KS. However, this most frequency-induced difference does not depend on the particular model. Therefore, for the sampled models in the inversion process, we calculate their predicted S3KS-SKKS differential traveltimes by TauP toolkit and then adjust them to the DSM results using the difference observed between the TauP and DSM results for model PREM, i.e.

PREM

DSM Taup

( ) ( ) ( ) ,

i i i i

g m  T^m  T^m T (4.7)

where (T_Taup^m )_i are the TauP-calculated differential traveltimes for model m, and

PREM

Ti

 is the difference between TauP and DSM differential traveltimes for PREM

PREM PREM PREM

Taup DSM

( ) .

i i

T T T

     (4.8)

Fig. 4.1 Difference in predicted S3KS-SKKS differential traveltimes obtained by TauP (green) and DSM (purple) for PREM in the distance range 120º-140º. The values are obtained by averaging every 0.5º over 1º bins.

(49)

The values of T_i^PREM shown in Fig.4.1 are estimated every 0.5º over 1º bins in the distance range 120°~140°. The differential traveltime data used in the final Bayesian inversion are selected from the 606 measurements displayed in Fig.3.3 after removing the outliers and the single data within a 1º bin (τi=0). In the end a total of 489 S3KS-SKKS differential traveltimes are used, as shown by the black dots in Fig.4.2a, and the standard deviation in Fig.4.2b are used for τi in eq.(4.5). Based on the modeling results for PREM and several PREM-like models in Chapter 3, we set the a priori probability distribution (m) to be the uniform distribution with a wave speed perturbation of -0.6% ≤ δ/ ≤ +0.2% with respect to  in PREM in the top 550 km layer of the outer core. We also impose a cosine taper so that the perturbation δ/

decreases smoothly to 0 at the depth of 550 km below the CMB. The model is parameterized by  at 20 depth grids over the 550-km layer (i.e. M=20).

Fig. 4.2 (a) 489 differential traveltimes measurement (block dots) used in the final inversion. (b) The measurements in (a) averaged every 0.5º over 1º bins. The bars show the standard deviations used for τi in eq.(4.5) to evaluate the likelihood.

(50)

4.4 A slower outer core model from Bayesian inversion

We carried out the Bayesian inversion of the 489 observed S3KS-SKKS differential traveltimes on an average PC. The sampling is still very time consuming even with the efficient MHMC algorithm, taking about 2 minutes for each model on average. After continuously running for about 7 days, 5092 models were sampled, with 109 accepted models. Repeated estimation of the model expectation shows little change after sampling more than 3600 models with 50 accepted models, suggesting that the MHMC sampling has stabilized. Fig.4.3 shows the result of Bayesian inversion from the 5092 model samples. The model expectation (red solid line) is computed by the assemble average of the 109 accepted models weighted by their respective a posteriori probability values.

Fig. 4.3 109 accepted models (gray) by the MHMC algorithm in Bayesian inversion.

The black, green and red solid lines are PREM, the maximum likelihood, and the expectation models. The red dash lines show the standard deviation of accepted models.

The black dashed lines show the boundaries of the prior model: the a priori probability distribution, i.e. (m) in eq.(4.4), has a uniform value of 1 between the black dashed lines and 0 otherwise.

(51)

This expectation model in Fig.4.3, named as PREM-SOC, has an average velocity perturbation of about -0.11% from PREM in the top 550 km of the outer core, and the perturbation is variable with depth. The maximum perturbation is about -0.31% at 60 km below the CMB, while there is a narrow depth range around 150 km below the CMB where the velocity is close to PREM.

The histogram in Fig.4.4 shows the comparison between the observed S3KS-SKKS differential traveltimes and predictions by model PREM-SOC. The values of the mean and standard deviation are not as small as those for Models PREM-OC2 and PREM-OC3 (Fig.3.13), although their average perturbations of the wave speed are all about -0.1%. As discussed earlier, due to the inherent uncertainty in the data and the uneven and incomplete coverage, fitness to the data is not the only criteria in determining the outcome of the Bayesian inversion, and the model expectation and its standard deviation are more meaningful results from the inversion.

Fig. 4.4 Histogram of the double-differential times for the PREM-SOC ²T^PREM-SOC. In this comparison the predicted differential traveltimes are obtained from DSM-calculated waveforms. The values of the means and standard deviations are given in the plot.

(52)

In Fig.4.5 we compare the expectation model from our Bayesian inversion with the Model KHOCQ of Helffrich & Kaneshima (2010). In the top 300 km of the outer core, our expectation model PREM-SOC has a good agreement with Model KHOCQ at long wavelength. The low-velocity zone in our model extends to a greater depth, and our model has a high-velocity layer ~150 km below the CMB. Our inversion result has a larger standard deviation than the error bounds of KHOCQ, perhaps due to different definitions between our standard deviation and error bounds in Helffrich & Kaneshima (2010). However, both results clearly indicate a lowered wave speed at the top of the outer core than that in PREM.

Fig. 4.5 Comparison of our inversion result with Helffrich & Kaneshima (2010). Solid red line is our expectation model PREM-SOC and red dash lines show the standard deviation. Solid blue line is Model KHOCQ in Helffrich & Kaneshima (2010), and the blue dashed lines are its error bounds. Black dashed lines show the prior in this study.