RECEIVER FUNCTION ANALYSIS FOR CRUSTALSTRUCTURE BENEATH PAKISTAN
Riaz Ahmed Soomro∗
MEE08161

Supervisor: Tatsuhiko HARA∗∗

ABSTRACT
We performed a receiver function analysis to study crust and uppermost upper mantle structure beneathPakistan. We used data recorded between June 2008 and June 2009 by a newly installed broadband
network of Micro Seismic Studies Programme. We selected events whose magnitudes are greater than 5in the epicentral distance range from 30 to 90 degrees. We found that amplitudes of tangential
components are comparable to those of radial components for some stations. We investigated theorientations of horizontal components by varying rotation angles. The results suggest that the deviationof the orientations of the horizontal components from the North-South and East-West directions is minor,
and that the large observed tangential components are likely to be caused by structural effects such as
dipping of interfaces.
We calculated radial receiver functions using extended time multitaper frequency domaincross-correlation receiver function estimation method. We selected stations for which receiver functions
with high signal-to-noise ratios are obtained among the 17 stations. The obtained receiver functions show
a significant azimuthal dependence, which suggests the complexity of the structures beneath thesestations. This is consistent with the large tangential components mentioned above. We performed
inversion of the obtained receiver functions to determine velocity structures of crust and upper mantle
using a genetic algorithm inversion technique. The results suggest that the Moho beneath the stations in
the Potwar region is thick (around a depth of 50 km), while it becomes thinner in the center of Pakistan.
This is qualitatively consistent with previous studies.
Keywords: Receiver function, Crustal Structure, Genetic algorithm.

1. INTRODUCTION
Pakistan lies in a seismically active region. The creation of Himalayas is due to subduction of the Indian
plate beneath the Eurasian plate. The relative speed between them is around 40 mm/y at present. Due tothis movement there are many active faults in and around Pakistan. It is essential to study crustal
structure and obtain reliable models for accurate hypocenter determination and understanding
geophysical and geological processes in crust and upper mantle.
There are several previous studies that investigated seismic structures of crust and upper
mantle in and around Pakistan by such methods as seismic tomography, relocation of hypocenters,
receiver function method, and so on.
Micro Seismic Studies Programme has upgraded its network from that consisting of shortperiod sensors to broadband sensors in 2006. There are 17 stations running at present; the number ofstations will increase in the near future.

∗

Micro Seismic Studies Programme, Islamabad, Pakistan.
International Institute of Seismology and Earthquake Engineering, Building Research Institute, Japan.

∗∗

12. THEORY AND METHODOLOGY
2.1 Receiver Function method
Seismic waves arrive at station carry information about earthquake source, propagation path and earth
structure. To determine velocity structure beneath a station it is necessary to isolate earth structure
information from the other factors. Langston (1979) provided such an isolation procedure. The receiver
function method is one of the techniques to study crustal structure. It is a time series computed fromthree component seismograms and provides information on earth structure beneath a seismic station(receiver). P to S converted phases which arrive after an arrival of a direct P-wave contains informationon the earth structure. Timings and amplitudes of these waves provide information on structure beneathseismic stations. Park and Levin (2000) developed a RF computation algorithm using Multiple - Taper
spectral correlation. Helffrich (2006) improved this method so that it is applicable to a longer time series.
However, it preserves phase information by using a sequence of short multiple tapers. He used a seriesof windows with 50% window overlap. Shibutani et al. (2008) improved this technique and used 75 %window overlap. They showed that the resulting taper has a flat level.
2.2 Receiver Function Inversions
Inversion of receiver functions for structure is non linear problem. There are two approaches to solve
non linear inverse problems; one is local optimization, and the other is global optimization technique.
In the former, the non linear problem is linearized first and (partial) derivative of observables isnecessary. In this approach, local minima can be a problem. An initial model should be close to thetrue model and a final model depends upon this initial model. Genetic algorithm is a method among
global optimization techniques. It is a stochastic search technique; it has resemblance with biologicalevolution systems. Genetic algorithms do not require (partial) derivatives of observables for inversion.
It is much freer from the local minimum problems; it is possible to treat complex problems which havemultiple minima. Therefore, it is suitable for inversion of receiver functions for crust and uppermost
upper mantle structure. The genetic algorithm consists of the following three steps:
1. Selection and Reproduction step.
2. Crossover step.
3. Mutation step.
The procedure is started with randomly generated initial models. Sambridge and Gallagher (1993)
used a direct ‘misfit’ based scheme in the reproduction step, where interim population of parents are
generated by selecting models with probability based on misfit calculations. We have used tournament
based selection scheme following Shibutani et al (1996). Initially there are Q model strings. Afterreproduction the parent population is paired to form Q/2 model strings and then these pairs are again
selected at random and their probability is calculated and if the probability is less than the probabilityof cross over (Pc) then both parents are crossed over to form the off spring population, otherwise they
are passed to next stage without any change (See Sambridge and Drijkoningen (1992) for details).
After crossover step the final step is mutation. In biology mutation brings diversity in new generations,and in genetic algorithm this is achieved by changing the parity of one or more chosen bits. The
mutation probability is kept small (in our case it is 0.01). After all three steps a new population of Q
bit model strings is produced, which will be the parent population and this process is repeated until aconvergence criterion is met.

3. DATA
3.1 Data Preparation
Micro Seismic Studies Programme has replaced the old analogue stations to the new broadband digital
network; at present 17 stations are running, and the number of stations will increase in the near future.

2

The seismic stations are located mainly in the northern part of Pakistan. The sensors are CMG-3T,
CMG-3ESP, and CMG-5T from Guralp Systems Limited. The observation period of the data used inthis study is that from June 2008 to June 2009. We selected 100 teleseismic events. We selected
teleseismic events whose magnitudes are greater than 5 with the epicentral distance in the range between
30 and 90 degrees. We used the newly installed broadband network data. The data format is GCF(Guralp Compressed Format). We converted data format from GCF to SAC using the utility “gcf2sac”,
which is available at Guralp website. After conversion we changed file names; the extension “*.sac” isreplaced by “*.Z”, “*.N”, “*.E” for vertical, north-south, and east-west components, respectively.
Header information for event and station locations in SAC headers is set.
3.2 Data CheckIn calculation of receiver functions, it is assumed that structure beneath a station can be expressed byhomogeneous horizontal layers. Thus, radial receiver functions are expected to be dominant. While
processing the data, we noticed that the tangential motions are comparable to radial motions for the data
recorded at the MCET station. There are two possibilities for these large tangential motions. The one isthat misorientation of horizontal components causes them. The other is that these tangential motions arecaused by structural effects such as scattering, dipping interfaces, and anisotropy. In the former case,radial motions will be dominant when the orientation is corrected. Therefore, we checked the orientationof the horizontal components of MCET. First, we performed rotation of horizontal components with anangle from backazimuth -90° to backazimuth +90° with an interval of one degree. Then we calculated
the amplitude ratio between the radial and tangential components and found an angle which maximizes
this ratio, which we call θmax hereafter. Even we changed the rotation angle from the back azimuth; theamplitude ratio was not large for many events. This suggests that the observed large tangential motions
are not due to misorientation but due to complexity of the structure beneath this region. Figures 1 and 2
show comparison of θmax and back azimuths, and frequencies of differences between them, respectively.
Although there is a scatter, we found the linear relationship between them (Fig. 1), and that theirdifferences scatter around 0, which also suggests that the observed large tangential motions are notcaused by misorientation. The mean and median of their differences are -4.51 and -5.02, respectively.
This suggests small deviation (around -5 degree) of the orientations of horizontal components from theNorth-South and East-West directions.
350

300

Amplitude Ratio

250

200
Series1

θmax

150

100

50

0
0

50

100

150

200

250

300

350

Back Azimuth

θmax.

Figure 1. Comparison of θmax andbackazimuth.

Figure 2. Frequencies of the
Differences of backazimuth and θmax.

3

4. RESULTS
4.1 Receiver functions calculation
We calculated radial receiver functions for the 17 seismic stations using extended time multitaper
frequency domain cross-correlation receiver function estimation method. We selected five stations,
MNIL, MSAG, MFAG, MNPR, and MCET, for which receiver functions have relatively highsignal-to-noise ratios among the 17 stations. Here we show receiver functions for MNIL, and MCET
for which a relatively large number of receiver functions are obtained. Receiver functions are aligned
according to their back azimuths. The red and blue colors represent the positive and negative
amplitudes, respectively. Figure 3 shows the receiver functions obtained for MNIL, although thenumber of the obtained receiver functions is small They have shown azimuthal dependence. Therefore,
for further analyses such as inversion for structures, it is necessary to divide receiver functions intosome groups based on their similarity, back azimuths, and epicentral distances. For inversion using thegenetic algorithm, we selected two receiver functions with high signal-to-noise ratios for the events
that occurred close to each other near the Kuril Island on 15/01/2009, and 07/04/2009. We stacked
Receiver functions for these two events in order to achieve better signal to noise ratio. Figure 4 shows
the receiver functions results obtained for MCET. We found that the receiver functions to be modeledby preliminary inversion using the genetic algorithm for this station are difficult; we do not show theirinversion results.

Figure 3. Receiver functions calculated for
the MNIL station. The receiver functions
are ordered by their backazimuths.
Positive and negative amplitudes areshown in red and blue color, respectively.

Figure 4. Receiver functions calculated for the
MCET station. The receiver functions areordered by their backazimuths. Positive andnegative amplitudes are shown in red and blue
color, respectively.

We have applied genetic algorithm developed by Shibutani et al (1996) which was originally based on
Sambridge and Gallagher (1993). The model parameters are layer thickness, shear wave velocity,
Vp/Vs ratio and density of each layer. We model the crust and upper mantle using six layers; asediment layer, basement layer, upper crust, middle crust, lower crust and uppermost upper mantle.
Based on the studies by Johnson and Vincent (2002) and Vinnik et al. (2007), the crustal thickness inthe study region is expected to be in the range between 30 to 60 km. The model space to be searchedcovers this range. The radial receiver functions for five stations MNIL, MSAG, MFAG, MNPR andMCET were inverted using this technique. We found that the receiver functions of the MCET and
MNPR stations were difficult to be modeled.

4

(a)

(b)

(c)

(d)

Figure 5. Inversion result for the MNIL station. (a) The best model. (b) Comparison between theobserved receiver function and the synthetic receiver function computed for the model shown in (a).
The dotted and solid curves correspond to the observed and synthetic receiver functions, respectively.
(c) The models for which the misfits are comparable (up to 5 percent difference) to that for the bestmodel. (d) The Moho depths beneath three stations determined in this study.
Figure 5 shows the result of receiver function inversion for the MNIL station. The best model (Fig. 5a)shows that the Moho depth is about 52 km. Figure 5(b) shows comparison between the observed
receiver function and the synthetic receiver function computed for the model shown in fig. 5a. Somesignals with large amplitudes are reproduced by this model. Figure 5(c) shows the models for which
the misfits are comparable (up to 5 percent difference) to that for the best model. Although there is large
uncertainty for velocity structure, the Moho depth is relatively well constrained, and their estimates arearound 52 km. In fig. 5d the Moho depths obtained for MNIL, MFAG and MSAG station are shown.

5. DISCUSSION
We applied the genetic algorithm using receiver functions to determine the crustal structure beneaththree stations. Our results suggest a thicker crust beneath MNIL and MFAG stations with the thicknessof 52 and 44 km., respectively. The crust is relatively thinner beneath MSAG (~37km.). This is firsttime that structure beneath Pakistan is studied using local seismic data. Our results are consistent withprevious studies that crust is getting thinner from the north to the south in Pakistan.
Johnson and Vincent (2002) developed a 3-D Velocity model for India and Pakistan region by
integrating more than 60 previous studies. They developed a regional 1D model based on their results,
in which the average crustal thickness is 50km. They also showed the lateral variation of crustal

5thickness in the region MCET is located in an area where there is a rapid lateral variation in crustal
thickness. This may be a reason for the observed large tangential component amplitudes and for that it
was difficult to model the observed receiver functions. The crust beneath MNIL is a little thicker, as
compared to Johnson and Vincent (2002). Our results are consistent with their results considering a
relatively small data set used in this study.
Vinnik et al. (2007) performed joint inversion of P and S receiver functions using data recorded at
the NIL station of the Global Seismograph Network. In their model, the crustal thickness beneath the
NIL station is 58±2 km, and the S wave velocity in the intermediate and lower crust is around 4.0 km/s.
Our results are qualitatively consistent with their results, although their thickness and S wave velocity
are larger than our estimates. One of the possible reasons for this difference is that they used both P
and S receiver functions. It is necessary to enlarge our dataset to efficiently perform stacking to
improve our models for further discussions. Also, it may be necessary to use S receiver functions.
They also determined deeper structures up to 300 km including anisotropy as a parameter to be
determined. Such a study is a possibility for extension of the present study.

6. CONCLUSIONS
We performed a receiver function analysis to study crust and uppermost upper mantle structure beneathPakistan. We found that amplitudes of tangential components are comparable to those of radial
components for certain stations. We investigated the orientations of horizontal components by varying
rotation angles. The results suggest that the deviation of the orientations of the horizontal components
from the North-South and East-West directions is minor, and that the large observed tangential
components are likely to be caused by structural effects such as scattering, dipping of interfaces andanisotropy. We also noticed azimuthal dependence from the receiver functions beneath these stations,which suggests the complexity of the structures beneath these stations.
In conclusion we suggest that the Moho beneath the seismic stations in the Potwar region is
deep at a depth of around 50 km. and that the crust is thinner in the central part of Pakistan. They arebasically consistent with previous studies. These results are preliminary attempt to determine structuresof crust and upper mantle using data from the newly installed broadband network. Accumulation ofdata in the near future, enhancement of signal-to-noise ratios and addition of new broadband stations
which cover the whole region in Pakistan will make it possible to determine structures beneathPakistan and we will be able to propose a better velocity model for our region by application of thedata analysis procedure of this study. This study provides a basis for such further studies in our region.

ACKNOWLEDGEMENT
I express my sincere gratitude to Dr. Shibutani for allowing me to use his codes and suggestions andvaluable comments and also to Dr. N. Hurukawa for his valuable comments to improve this manuscript.

REFERENCES
Helffrich, G., 2006, Bull. Seism. Soc. Am., 96, 344-347.
Johnson, M., and Vincent, C., 2002, Bull. Seism. Soc. Am., 92, 2893-2910.
Langston, C. A., 1979, J. Geophys Res., 84, 4749-4762.
Park, J. and Levin, V., 2000, Bull. Seism. Soc. Am., 90, 1507-1520.
Sambridge, M., and Drijkoningen, G., 1992, Geophys. J. Int., 109, 323–342.
Sambridge, M., and Gallagher, K., 1993, Bull. Seism. Soc. Am., 83, 1467–1491.
Shibutani, T., Sambridge, M., and Kennett, B., 1996, Geophys. Res. Lett., 23, 1829–1832.
Shibutani, T., Ueno, T., and Hirahara, K., 2008, Bull Seism. Soc. Am., 98, 812-816.
Vinnik, L., Singh, A., Kiselev, S., Kumar, M. R., 2007, Geophys. J. Int., 2007, 171, 1162-1171.

6