DETERMINATION OF CRUSTAL STRUCTURE
BENEATH A BROADBAND STATION IN MALAYSIAUSING RECEIVER FUNCTION ANALYSIS.
Zahratul Ain Din ∗MEE09180

Supervisor: Tatsuhiko HARA ∗∗ABSTRACT
We performed receiver function analyses in order to investigate crust structure beneath a broadband
seismic station in Malaysia. We selected 11 teleseismic events from the Malaysian National Seismic
Network Database whose magnitudes are in the range between 5.8 and 7.1 with good signal-to-noise
ratios. We compared the observed receiver functions to the synthetic receiver functions computed for amodel for the site of station IPM taken from the global crust model, CRUST 2.0. There is a significant
difference between them. This is likely to be due to the thin sedimentary layer in the crust model, andour comparison suggests that such a sedimentary layer does not exist beneath the station site. Then, weapplied a genetic algorithm to perform inversion. We modeled the crust and uppermost upper mantle
with six major layers: sediment layer, basement layer, upper crust, middle crust, lower crust anduppermost mantle. The Moho depth is relatively well constrained and is estimated to be about 35km.
The thin sedimentary layer is not obtained in the inversion result, although a thin layer with relatively
low S wave velocity is obtained. This study is a preliminary attempt to determine velocity structuresand obtain an appropriate velocity model for Malaysia. Further accumulation of data and stacking
using a larger dataset will be useful to construct such a model by application of the data analysis
procedure of this study.
Keywords: Receiver Functions, Crustal Structure, Genetic Algorithm.

1. INTRODUCTION
Malaysia is tectonically situated within relatively stable Sunda-land, and considered as acountry with a very low seismicity. It is bordered to the west and to the east by the Indonesia(Indo-Australian Plate) and the Philippines (Philippines Sea Plate), which are two of the mostseismically active countries in this region with frequent earthquakes. The Great Sumatra Earthquakewith a magnitude of 9.3 on December 26th, 2004 caused a huge tsunami in the Indian Ocean that hit
many countries including Malaysia. From that historical moment, Malaysian Government has decided
to establish its own monitoring center. The main purpose of this study is to determine the crustal
structure beneath a broadband station in Malaysia using a seismic data from the newly establishedseismic network.

2. DATATo determine crustal structure, we chose a broadband station installed at Ipoh (its station code is IPM)
in Malaysian Peninsula. We retrieved waveform data for 11 events from the Malaysian National
Seismic Network Database whose magnitudes are in the range between 5.8 and 7.1 that occurred in the∗

Malaysian Meteorological Department (MMD), Malaysia.
International Institute of Seismology and Earthquake Engineering, Building Research Institute, Japan.

∗∗

1teleseismic distance range (Table 1). In order to obtain an effective receiver function modeling, wechose data with good S/N (signal to noise ratio).
Table 1. Teleseismic events used in this study.
Mag
(Mw)
Date
2/18/2010 6.9
2/07/2010 6.3
2/26/2010 7.0
3/08/2010 6.1
1/27/2010 5.9
3/20/2010 6.6
2/01/2010 6.2
1/03/2010 6.6
1/05/2010 6.8
1/03/2010 7.1
4/11/2010 5.8

Depth
(km)
573.8
16.5
22.0
446.5
24.7
423.5
33.0
10.0
35.0
30.5
35.0

Epicentral
Azimuth distance(degree) (degree)
30.15
46.50
46.91
28.98
47.97
34.03
66.99
44.42
67.46
25.89
97.88
51.72
100.8054.62
103.4055.02
103.4358.17
103.5056.99
106.1067.25

Slowness
0.071
0.081
0.079
0.072
0.082
0.068
0.066
0.064
0.064
0.064
0.058

Location
China-Russia-North Korea
Southwestern Ryukyu Island
Ryukyu Islands, Japan
Maug Islands Reg, N.Mar
Philippine Islands Region
New Ireland Region, P.N
Bougainville Region, P.N
Solomon Islands
Solomon Islands
Solomon Islands
Santa Cruz Islands

3. THEORY AND METHODOLOGY
3.1. Receiver Functions
A waveform is a recorded time series of the earth movement beneath a seismic station (receiver). The
relative response can be computed using three-component seismograms, vertical (Z), East-West (E-W),
and North-South (N-S). Teleseismic waves arriving at a seismic station from different parts of theearth contains different information on the earth structure. Receiver function is one of the techniques
to study the crustal structure. In receiver function analysis, information regarding earthquake source
and the propagation through the mantle is removed from observed waveform data. We can extract
information on structure of the earth’s crust and uppermost upper mantle beneath a seismic station by
analyzing the vertical and radial component records.
3.2. Data Preparation for Receiver Function Analysis
For receiver function analysis, we need three-component observations with a wide bandwidth period.
We use SAC (Seismic Analysis Code) program to set several header variables (event latitude,
longitude & depth and station latitude & longitude) in each waveform. The information is necessary
for SAC program to rotate horizontal seismograms into the radial and tangential components by anangle of the back azimuth. The final data-preparation consists of windowing the P waveform from thepre-signal noise and the rest of the seismic signal. We are using signal of 60s before and after thearrival of the P wave. To avoid signal processing artifacts later in the processing, we remove the mean
and taper the ends of signal. We used a program developed by Dr. Charles J. Ammon which isavailable at http://eqseis.geosc.psu.edu/~cammon/HTML/RftnDocs/prep01.html
3.3. Calculation of Receiver Function
The receiver function analysis uses the converted phases and multiples recorded on the horizontal
seismograms (e.g., Langston 1977, 1979, 1981). It is necessary to isolate earth structure informationfrom other factors. Teleseismic waves include conversions and reverberation phases generated at
discontinuities beneath each station (Langston, 1979). Langston (1979) developed a source
equalization procedure to remove the effects of near-source structure and source time functions. The
impulse response of the earth structure for vertical components is assumed to be the Dirac delta
function. We calculate receiver functions following Owens et al. (1984). To compute receiver
functions, the vertical component signal is deconvolved from the radial component to produce source
equalized radial seismogram (Ammon, 1991). For deconvolution, we used a water-level stabilization
method and a low-pass Gaussian filter to remove high-frequency noise. A time series computed by this2

deconvolution is called a radial receiver function. The final expression for the radial receiver function
in frequency domain is as follows

ER (ω ) =DR (ω ) DV (ω )G (ω )ϕ (ω )

Where ER(ω) is Fourier transform of the radial receiver function, DR(ω) is Fourier transform of theradial component of motion and DV(ω) is Fourier transform of the vertical component of motion.

[]ϕ (ω ) = max{DV (ω ) DV (ω ), c max DV (ω ) DV (ω ) } and G (ω ) = e

−

ω2
4α 2

(Gaussian Filter)c = water-level
3.4. Genetic Algorithm in Receiver Functions Inversion
Genetic algorithm is one Table 2: Model parameters in genetic algorithm receiver function inversion.
of global search method
C rustS ediment B ase ment
M antleupper
m iddlelower(e.g., Sambridge and
0.0
0.0
5.0
5.0
0.0
0.0
Thi ckne ss (km ) Lower
Drijkoningen,1992).
Upper
2.1
3.5
20.0
20.0
15.0
0.0
n
3
3
4
4
4
0
Genetic
algorithm
Incre ment
0.3
0.5
1.0
1.0
1.0
0.0
consists of reproduction
Lower
1.5
2.2
2.9
3.1
3.3
4
Vs(km/s)
Upper
2.2
2.9
3.6
4.0
4.0
4.7
step, crossover step, and
n
3
3
3
3
3
3
mutation step. Genetic
Incre ment
0.1
0.1
0.1
0.1
0.1
0.1
1.70
2.30
2.67
2.18
3.18
3.25
Density(g/cc) Lower
algorithms work with a
Upper
2.40
2.60
2.67
2.18
3.18
3.25
group of Q (size of
n
3
2
0
0
0
0
Incre ment
0.1
0.1
0.0
0.0
0.0
0.0
working
population)
Lower
1.80
1.80
1.73
1.73
1.73
1.73
Vp/Vs
models simultaneously,
Upper
2.50
2.10
1.73
1.73
1.73
1.73
n
3.0
2.0
0.0
0.0
0.0
0
initially chosen at random,
Incre ment
0.1
0.1
0.0
0.0
0.0
0.0
and code each into a Qα
100
675
1450
1450
1450
1450
Qβ
25
300
600
600
600
600
binary string.
We have modeled the crust and uppermost mantle with six major layers: sediment layer,
basement layer, upper crust, middle crust, lower crust and uppermost mantle. The model parametersin each layer are thickness, S wave velocity, density, and the velocity ratio between Vp/Vs. Thesemodel parameters are represented by a binary string. The length of the binary string is 46. For eachmodel parameter, upper limit, lower limit and 2n possible values are specified. The best model issearched within these limits. The total number of model parameters is 15. The ranges of the model
parameters, the n values and the incremental values are shown in Table 2. The total size of model
space to be searched is 246. The inversion procedure was iterated for 200 generations; as a result
200,000 models were generated. The radial receiver functions for IPM station were inverted using
this technique. We used the code developed by Shibutani et al. (1996) to perform inversion ofreceiver functions for crustal structure.

4. RESULTS AND DISCUSSION
4.1. Observed and Synthetic Receiver Function
Velocity structure in crust and uppermost upper mantle beneath seismic stations can be determined byinversion of receiver functions with time duration about 30-40 second (Owen et al., 1984). Before
performing receiver function inversion, we compared the observed receiver functions to synthetic
receiver functions computed for an existing crust model for the station site. In this analysis, we used
CRUST 2.0 (http://igppweb.ucsd.edu/~gabi/crust2.html) to compute synthetic receiver functions.
Figure 2 shows P and S wave velocity structure for IPM station. There is a thin sedimentary layer exist
with S wave velocity of 1.2km/s in this model. The Moho thickness beneath the IPM station is about
32.5km.

3

Figure 3 shows a comparison of the observed and synthetic receiver functions. There are
significant differences between them. In the synthetic receiver function there are several peaks in theearly part (say, 0-2 s), while there is one peak in the observed receiver function. These peaks in thesynthetic receiver function may be caused by the thin sedimentary layer in the model taken from
CRUST2.0. To investigate this possibility, we replace the sedimentary layer by the layer with theunderlying upper crust property to make another synthetic crust model. Then we recalculated thesynthetic receiver function for this modified model. The computed receiver function is shown in
Figure 3(c). The several peaks found in the receiver function shown in Figure 3(b) disappear. This
result suggests that these peaks are produced by the thin sedimentary layer. Since the receiver function
computed for the modified model (i.e., there is no thin sedimentary layer) is more similar to theobserved receiver function, the thin sedimentary layer is not likely to exist beneath station IPM.

Figure 2. P and S wave velocity models for IPM
station taken from CRUST 2.0 model are shownin the left and right panels, respectively.

Figure 3. Radial receiver functions (a)observed, (b) synthetic computed for the model
shown in Figure 2, and (c) the modified model,
for the 18/2/2010 event that occurred in
China-Russia-North Korea.

4.2 Receiver Functions Inversion by Genetic Algorithm
Among the 11 receiver functions that we calculated, we chose those for the events that
occurred on 2/07/2010 and 2/26/2010 to perform inversion based on signal-to-noise ratios, theirsimilarity, back azimuths and epicentral distances. Then, following the suggestions by Dr. Ammon
the receiver functions are down-sampled from 20 Hz to 10 Hz, and the signals up to 20 second fromthe P onset are used for inversion. Figure 4 shows the best model for which the residual between theobserved receiver function and the synthetic receiver function is the smallest. There is no sedimentarylayer in this model, while there is a thin layer with a relatively low velocity. Figure 5 shows
comparison between the observed receiver function and the synthetic receiver function computed forthe model shown in Figure 4. Some significant features are reproduced, although there are several
differences. Figure 6 shows the models for which the residuals between the observed receiver
functions and synthetic receiver functions are comparable (up to 2 per cent difference) to that for thebest model shown in Figure 4. The Moho depth seems relatively well constrained to be around 35 km,
while the density and the Vp/Vs ratios of the shallow layers are not well constrained.
To quantitatively evaluate the uncertainty of the inversion result, we calculated marginal
posteriori probability density functions for each model parameter using all of the models generated in
the inversion. Figure 7 shows the probability of thickness for each layer. The values and shapes ofprobabilities indicate whether the parameter in each layer is well constrained or not. We can see thethickness in each layer is well constrained with a value of Moho depth about 35km. This is consistent
with the results show in Figure 6.

4

Figure 4: The best model obtainedby inversion for IPM station.

Figure 6. The models
shows
the
residuals
between the observed
receiver functions andsynthetic
receiverfunctions are comparable
to that for the best model.

Figure 5. Comparison between thesynthetic receiver function and theobserved receiver function.

Figure 7. Marginal posteriori probability density functions for thickness in each layer.

55. CONCLUSIONSIn this study, we performed receiver function analysis to determine the crustal structure beneath thebroadband station IPM in Malaysia. We used data from the 11 teleseismic events with magnitudes in
the range between 5.8 and 7.1 from the Malaysian National Seismic Network Database. We calculated
receiver functions by deconvolving vertical components from radial components with water levels to
make computations stable. Before performing inversion for structure, we computed synthetic receiver
functions for IPM station. The crust model is chosen from the global crust model, CRUST 2.0; themodel contains a thin sedimentary layer with S-wave velocity of 1.2km/s. There is a significant
difference when we compared synthetic receiver functions to the observed receiver functions. The
differences suggest that such a sedimentary layer does not exist beneath the station.
Then, we applied a genetic algorithm to perform inversion of the observed receiver
functions for crustal structure beneath IPM station. We chose the receiver functions from two events
for inversion based on signal-to-noise ratios, their similarity, back azimuths and epicentral distances.
There is no sedimentary layer in the obtained model, while there is a thin layer with a relatively low
velocity. Some significant features of the receiver functions are reproduced by the obtained model. We
calculated marginal posteriori probability density functions for each model parameter using all of themodels generated in the inversion. The thickness for each layer, the S-wave velocities in the crust andthe Moho depth is relatively well constrained, which suggests that the Moho depth beneath IPM
station is around 35km. Constraints for Vp/Vs ratios and densities in the sedimentary and basement
layer seem weak.
This is the first receiver function analysis done using data from the Malaysian National
Seismic Network. When we compare our results to the model from CRUST2.0 a significant difference
is found for the very shallow structure beneath the IPM station. Our results suggest that there is a thinlayer with a relatively low S-wave velocity and that the sedimentary layer in the CRUST2.0 model
does not exist. Further accumulation of data and application of the technique used in this study will be
useful for construction of a crust model appropriate for Malaysia.

ACKNOWLEDGEMENT
I would like to present my sincere appreciation to Director of IISEE Dr. Shoichi Ando, to Seismology
Course Leader Dr. Toshiaki Yokoi and Dr. Bunichiro Shibazaki for their valuable comments toimprove my dissertation. I express my sincere gratitude to Dr. Shibutani for allowing me to use hiscodes and Dr Charles J. Ammon for the useful instruction in his webpage which I used as a reference
to perform receiver function analysis.

REFERENCES
Ammon, C. J., 1991, Bull. Seism. Soc. Am., 81, 2504-2510
Langston, A. C., 1977, Bull. Seism. Soc. Am., 67, 1029-1050
Langston, A. C., 1979, J. Geophys. Res., 84, 4749-4762
Langston, A. C., 1981, J. Geophys. Res., 86, 3857-3866
Owens, T. J., Zandt, G., and Taylor S. R., 1984, J. Geophys. Res., 89, 7783-7795
Sambridge, M., and Drijkoningen, G., 1992, Geophys. J. Int., 109, 323-342
Shibutani, T., Sambridge, K., and Kennett, B. L. N, 1996, Geophys. Res. Lett., 23, 829-1832
Web page: http://eqseis.geosc.psu.edu/~cammon/HTML/RftnDocs/rftn01.html
Webpage: CRUST2.0; http://igppweb.ucsd.edu/~gabi/crust2.html

6