Crustal thickness and velocity structure beneath Singapore’s seismic network

Crustal thickness and velocity structure beneath Singapore’s seismic network

Journal of Asian Earth Sciences 64 (2013) 245–255 Contents lists available at SciVerse ScienceDirect Journal of Asian Earth Sciences journal homepag...

2MB Sizes 3 Downloads 131 Views

Journal of Asian Earth Sciences 64 (2013) 245–255

Contents lists available at SciVerse ScienceDirect

Journal of Asian Earth Sciences journal homepage: www.elsevier.com/locate/jseaes

Crustal thickness and velocity structure beneath Singapore’s seismic network Kenneth A. Macpherson a,⇑, Dannie Hidayat a, Lujia Feng a, Siang Huat Goh a,b a b

Earth Observatory of Singapore, Nanyang Technological University, Singapore Department of Civil Engineering, National University of Singapore, Singapore

a r t i c l e

i n f o

Article history: Received 1 August 2012 Received in revised form 12 December 2012 Accepted 19 December 2012 Available online 4 January 2013 Keywords: Receiver functions Crustal thickness Singapore Velocity structure

a b s t r a c t We estimated the crustal thickness and velocity structure beneath the five stations comprising the Republic of Singapore’s seismic network. Our data set was composed of 697 teleseismic receiver functions and 7 months of broad-band data that was cross-correlated to produce inter-station Green’s functions. Surface wave group velocities were extracted from the Green’s functions to obtain dispersion data for a path from central Sumatra to Singapore in order to provide a complimentary data set to the receiver functions. Crustal thickness was estimated via an H  k stacking technique, and high-resolution 1D Pwave velocity profiles were generated beneath each station by jointly inverting receiver function stacks and the group velocity data using a linearised time-domain inversion scheme. Crustal thickness beneath four stations was found to be between 28.0 km and 32.0 km, while one station in the northeast of Singapore indicates 24.0 km thick crust. This implies a significant crustal thinning beneath Singapore over the lateral extent of 50.0 km. Inversion results exhibit several crustal features that are observable in the derived models at all five stations, indicating that they exist across Singapore as a whole. There appears to be an upper-crustal high-velocity zone beneath Singapore, underlain by a velocity inversion. Station NTU shows slower near-surface velocities than the other stations, consistent with its situation above the sedimentary Jurong formation. These results expand the available global velocity data set, as well as being useful for assessing the seismic hazard in Singapore. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction The Republic of Singapore, a highly-urbanised city-state with a population of 5.2 million, is located at the southern tip of the Malay peninsula. The nation occupies 63 islands with a combined area of roughly 700.0 square kilometres, of which some 77% contains significant built infrastructure. The seismic hazard in the southern Malay peninsula is low; the peak horizontal acceleration with 10% probability of exceedance in 50 years being 5% g (Petersen et al., 2004). However, although Singapore lies some 700.0 kilometres (km) from the nearest tectonic margin, long-period ground motions induced by large events on the Sumatran subduction interface have historically been felt in the nation (Pan and Sun, 1996). More recently, the Mw = 8.6 April 11, 2012 earthquake off of the west coast of north Sumatra induced swaying of taller structures in the city that was clearly noticeable to occupants, and prompted some 126 reports of Modified Mercalli Intensity of 2 (USGS, 2012). Also, simulation studies indicate that high-rise structures on soft-soil sites in Singapore may be at risk from long-period ground-motions if they are of considerable duration (Megawati and Pan, 2009). Such ground-motions are likely in the event of a ⇑ Corresponding author. Tel.: +65 9383 9509. E-mail address: [email protected] (K.A. Macpherson). 1367-9120/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jseaes.2012.12.027

large (Mw > 8.0) earthquake on the Mentawai segment of the Sumatran megathrust. This segment, which is the closest part of the megathrust to Singapore, has not fully ruptured since 1797 (McCloskey et al., 2010), and super-cycle patterns indicate that it could produce a large event in the coming decades (Sieh et al., 2008). Because of the extant hazard of long-period ground motions to Singapore, it is important to understand the velocity structure beneath the city and how it might respond to seismic waves. Singapore’s Meteorological Service Division operates a seismic network consisting of one very-broad-band station (BTDF) housing a Streckeisen STS-2/VBB seismometer and four teleseismic stations (BESC, KAPK, NTU, PTK) housing Kinemetrics WR-1 seismometers located in a rough east–west line across the island (see Fig. 1). The frequency response of the STS-2 instrument is flat between 0.002 Hertz (Hz) and 10.0 Hz, while the WR-1 instruments are flat between 0.1 Hz and 20.0 Hz. Data from this network are shared with other regional networks and with the Global Seismic Network, where they are available through the IRIS consortium. Despite being operational since 1996, little work has been done to characterise the crustal thickness and velocity structure beneath the network (Walling et al., 2010). Such analysis is needed in order to develop velocity models for Singapore that will be useful for seismic hazard studies such as ground-motion simulations. Additionally, such velocity models will add to the existing global

246

K.A. Macpherson et al. / Journal of Asian Earth Sciences 64 (2013) 245–255

At right

° 60

30 °

Event location Location of Singapore (enlarged below)

Below Network MS seismic station

5 km piercing points BESC NTU

BTDF PTK

KAPK

60 °

PTK

1.4˚

3.6 mm/yr

1.4˚

BTDF BESC

KAPK

NTU

10 km 1.2˚ 103.6˚

1.2˚ 103.8˚

104˚

Fig. 1. Study area showing the locations of the five stations that comprise Singapore’s seismic network. The globe shows the locations of the teleseismic events that were used to compute receiver functions. The five kilometre piercing points of teleseismic observations, calculated by assuming a 6.0 km/s crustal velocity, are plotted, colour-coded for each station. The motion of the continuous GPS station located at NTU relative to the Sunda block is plotted in magenta, with the 95% confidence error ellipse. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

database of available crustal thicknesses and velocity profiles. In this study, we investigate the crustal thickness and velocity structure beneath Singapore’s seismic network by analysing teleseismic receiver functions and regional surface waves. 1.1. Geologic and seismotectonic setting The near-surface geology of the main island of Singapore is characterised by three broad regions. The south and southwest is primarily composed of sedimentary rocks of the Jurong formation, the central region consists of the Bukit Timah granite and the Gombak norite, while the east of the island is primarily covered by the Old Alluvium formation of Quaternary sediments (Rahardjo et al., 2004). Recent alluvial deposits of the Kallang formation are distributed throughout the island. Because of the regions humid tropical environment, significant weathering has occurred to these formations, resulting in extensive deposits of residual soils (Rahardjo et al., 2004). On a more regional scale, Singapore is located, along with Peninsular east Malaysia, within the Indochina terrane to the south east of the Bentong-Raub suture zone (Metcalfe, 2000). This terrane likely separated from Gondwana during the Devonian. Singapore is located within the tectonically stable Sunda shelf. This is the continental shelf of southeast Asia, and is considered to be distinct from the Eurasian plate due to significant relative motion between the two induced by the collision of India with Asia (Bock et al., 2003). A continuous GPS station located in the west of Singapore (see Fig. 1) indicates 3.6 mm per year of motion to the

northeast. The closest likely source zone for significant seismic events is the strike-slip Great Sumatran fault, whose nearest segment is more than 300.0 km away from Singapore. The largest earthquakes produced by the Sumatran fault are likely to be less than a moment magnitude of 8.0 (Natawidjaja and Triyoso, 2007). Therefore, the most salient seismic hazard to Singapore is a great earthquake (Mw > 8.0) on the Mentawai segment of the Sumatran megathrust.

2. Receiver functions The analysis of teleseismic receiver functions has become a routine tool to investigate the depth of discontinuities beneath 3-component seismic stations, and has been undertaken in diverse tectonic environments, e.g., Langston (1977), Ammon and Zandt (1993), Searcy et al. (1996), Dugda et al. (2007), Park et al. (2009), Lloyd et al. (2010), and Macpherson et al. (2012). Receiver functions are time series that have been derived by deconvolving the vertical components from radial components of 3-component recordings. The deconvolution removes the source factor, and because P-to-S converted waves have higher amplitudes in the radial component than the vertical, results in a time series that is primarily a function of the discontinuities beneath the station (Ammon, 1991). Receiver functions may be modelled to provide constraints on crustal thickness and depth to discontinuities and may be combined with complimentary data sets, such as surface-wave dispersion curves, in order to estimate sub-station velocity structure.

K.A. Macpherson et al. / Journal of Asian Earth Sciences 64 (2013) 245–255

2.1. Teleseismic data In order to compute receiver functions for Singapore’s seismic network, we chose events at teleseismic distances of between 30.0° and 90.0° in the National Earthquake Information Center catalogue. We also restricted the magnitude range to 5.0 6 Mw 6 8.0 so that the events would be large enough to help increase the signal to noise ratio of the recordings, but not so large as to exhibit considerable source complexity. Because the network consists of five stations that have been operational since 1996, these selection criteria resulted in a data set of several thousand recordings. In order to manage such a large data set, we computed receiver functions via the iterative deconvolution technique of Ligorria and Ammon (1999), which is well suited to automated processing of large data sets due to the method’s stability (Crotwell and Owens, 2005). We used trial and error to choose a Gaussian width parameter of 1.5 for use in the deconvolutions. This value attenuates frequencies greater than about 0.75 Hz, and allows a reasonable compromise between minimising noise while still providing sufficient detail. While the four stations housing WR-1 seismometers are not true broad-band stations, they do have a flat response down to 0.1 Hz, providing sufficient bandwidth to compute receiver functions given our choice of Gaussian width. Receiver functions computed from shorter-period instruments may not be as stable as those computed from broad-band data, and may exhibit negative side-lobes around positive peaks (Shaw-Champion et al., 2006). However, because we use the time-domain deconvolution, we do not experience stability issues, and we note that receiver functions obtained from short-period data have been used to successfully model crustal structure in previous studies (Tomlinson et al., 2006). Because of the large quantity of teleseismic recording available, we were able to be highly selective with our data. We examined the quality of each recording by convolving the computed receiver function with the corresponding vertical component and comparing the result to the original radial trace. Deconvolutions that resulted in a larger than 90.0% fit were retained. This criteria resulted in the retention of 148 observations for station BESC, 271 for BTDF, 30 for KAPK, 188 for NTU, and 60 for PTK. See Fig. 1 for the source locations and the 5.0 km deep piercing points of the teleseismic events used to compute receiver functions for this study. The radial receiver functions that were retained for analysis are shown in Fig. 2. The data are shown as a move-out record section with the plots organised by ray parameter. We indicate the theoretical arrival times of the converted Ps phase and the PpPs and PpSs + PsPs multiple phases, computed based on a crustal thickness of 30.0 km with an average crustal velocity of 6.0 km/s. Observations at station BESC (Fig. 2a) show a clear Ps arrival, with the PpPs and PpSs + PsPs phases less obvious. The presence of an intra-crustal discontinuity is indicated by a clear arrival between the direct P and the Ps. This phase is more distinct for events closer to Singapore that arrive at a shallower angle. Waveforms for larger ray parameter observations appear more complex in general for this station, and this indicates that they may be sampling more complex structure. The receiver functions for station BTDF (Fig. 2b) appear relatively simple compared to station BESC. The Ps phase as well as the multiples are easily visible. The crust appears to be relatively homogeneous beneath BTDF, as Ps is the first significant arrival following the P arrival. The simple waveforms observed at this station may be due to its situation above the Bukit Timah granite. Station KAPK (Fig. 2c) exhibits complex waveforms, and individual phases are difficult to identify. If we assume that the phase nearest to the Ps theoretical arrival time is indeed Ps, then there is a distinct intra-crustal arrival. Overall, station KAPK exhibits the

247

most complicated waveforms and seems to be sampling more complex structure than the other stations. Receiver functions at station NTU (Fig. 2d) show clear Ps and PpPs arrivals, while the PpSs + PsPs phase is less distinct. While there is no clear arrival from an intra-crustal converter, the Ps arrival does exhibit a side-lobe in many of the traces, indicating lateral velocity variations within the crust. Station PTK (Fig. 2e) receiver functions show clear Ps phases but the multiple phases are difficult to identify. There is not much evidence of intra-crustal converters, although a few Ps phases have side-lobes. 3. Crustal thickness and Vp/Vs ratio beneath Singapore The depth of the Mohorovicic discontinuity (Moho), or crustal thickness, is a property that is fundamental to understanding the crustal structure of a region. The topography of the Moho is also of interest from the standpoint of seismic hazard, as it can have a significant impact on ground-motions via the crustal wave-guide effect (Zhu and Kanamori, 2000). Teleseismic receiver functions are routinely used to constrain crustal thickness, the simplest method being to examine the delay in the Ps arrival after the direct P, denoted by tPs, via the following equation:

tP H ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 2   p  p2 V2 V2 s

ð1Þ

p

where H is thickness, p is the ray parameter, and Vs and Vp are the shear and compressional wave velocities, respectively. The limitation of this method is that there is a strong trade-off between crustal velocity and thickness, specifically, the thickness (H) is highly dependent on k = Vp/Vs. To avoid this problem, we use the H  k stacking method of Zhu and Kanamori (2000). This method, rather than consider the Ps phase only, also uses the timing of the multiple phases PpPs and PpSs + PsPs to transform the time-domain receiver function into H  k space for stacking, where H is crustal thickness and k is Vp/Vs. The stacking is defined as

sðH; kÞ ¼ w1 rðt 1 Þ þ w2 rðt 2 Þ  w3 rðt 3 Þ

ð2Þ

where the wi are weighting factors to account for the variable signal to noise ratios between the phases, r is the radial receiver function, and the ti are the predicted travel times for the three phases given a particular value for H and k. The best estimate of H and k is achieved when the three phases stack most coherently and s is at a maximum. For each of the five stations comprising Singapore’s seismic network, we applied the H  k stacking algorithm using all available receiver functions. Stacking observations with many different ray parameters and back azimuths is an advantage as it helps to reduce the effects of lateral variations, providing an average crustal thickness for the station. We used an interpolated CRUST2.0 model to determine an average crustal velocity for Singapore of Vp = 6.0 km/s, and this value was used for all stations. Because the signal to noise ratio is higher for the Ps phase than the later phases, we weighted it the most heavily, and decreased the weighting of the multiples. We experimented with many different weights, but found the results were unaffected as long as we chose a reasonable weighting scheme. We used a weighting of 0.7 for the Ps phase, 0.2 for the PpPs multiple, and 0.1 for the PpSs + PsPs multiple. The results of this analysis are shown in Fig. 3. The method yielded a crustal thickness estimate beneath station BESC of 28.0 ± 0.1 km, beneath station BTDF of 29.5 ± 0.0 km, beneath KAPK of 32.3 ± 0.7 km, beneath NTU of 32.0 ± 0.2 km, and beneath PTK of 24.0 ± 0.1 km. We have higher uncertainties for stations with fewer observations, so that we have a larger range of values for station KAPK, with only

K.A. Macpherson et al. / Journal of Asian Earth Sciences 64 (2013) 245–255

radial RF for BTDF, rel ampl

0.08

0.08

0.07

0.07

0.06

0.06

0.05

0.05

0.04

0.04

(c) Ray Parameter (s/km)

(b)

radial RF for BESC, rel ampl

(d)

radial RF for KAPK, rel ampl

Ray Parameter (s/km)

Ray Parameter (s/km)

(a)

radial RF for NTU, rel ampl

0.08

0.08

0.07

0.07

0.06

0.06

0.05

0.05

0.04

Ray Parameter (s/km)

248

0.04 0

10

20

30

40

0

10

time (s)

20

30

40

time (s)

(e)

radial RF for PTK, rel ampl

Ray Parameter (s/km)

0.08

0.07

0.06

0.05

0.04 0

10

20

30

40

time (s) Fig. 2. Plots of all radial teleseismic receiver functions that were analysed in this study. Traces are organised by ray parameter, and the theoretical arrival time of the Ps, PpPs, and PpSs + PsPs phases are indicated. Plots are shown individually for each station: BESC (a), BTDF (b), KAPK (c), NTU (d), and PTK (e).

30 receiver functions stacked, while station BTDF, with 271 receiver functions stacked, has a range that is less than j0.1j. Although our H  k analysis provides depth to the Moho at only five points across Singapore, we created an isopach map of crustal thickness in order to better visualise the topography of the Moho (see Fig. 4). This was done simply by interpolating the thickness values between the five stations using a variable tension spline. The amount of variation in crustal thickness over such a small (1000.0 km2) area is surprising. While the four stations on the

main island of Singapore have thicknesses within a kilometre or two of 30.0 km, station PTK on Ubin island differs significantly, with a thickness of 24.0 km. While PTK does have the second fewest number of observations out of the five stations, we feel that a 60 receiver functions stack is enough to provide a degree of confidence in the estimate. We conclude that the Moho beneath Singapore has an average dip from northeast to southwest on the order of 10°. Vp/Vs ratios vary between 1.68 beneath KAPK and 1.99 beneath PTK, where the crust is thinner. The average k for all five

K.A. Macpherson et al. / Journal of Asian Earth Sciences 64 (2013) 245–255

Station: MS-BESC

(a)

249

Station: MS-BTDF

(b)

Station: MS-NTU

Station: MS-KAPK

(c)

(d)

Station: MS-PTK

(e)

Fig. 3. Plots of results of the H  k stacking analysis for station BESC (a), BTDF (b), KAPK (c), NTU (d), and PTK (e). The greyscale colourbar shows the coherence of the Ps, PpPs, and PpSs + PsPs phases achieved by stacking with different crustal thicknesses (H) and Vp/Vs(k) ratios. The preferred values for H and k are those that result in the largest coherence.

250

K.A. Macpherson et al. / Journal of Asian Earth Sciences 64 (2013) 245–255

crustal thickness (km) 20

30

40

PTK

1.4˚

1.4˚ NTU

BTDF BESC KAPK

10 km 1.2˚ 103.6˚

1.2˚ 103.8˚

104˚

Fig. 4. Crustal isopach map for Singapore. The surface was created by interpolating the crustal thickness values estimated at each station from the H  k analysis with a regularised spline with tension.

Singapore stations is 1.83, with is higher than the global average of 1.78 (Zhu and Kanamori, 2000). 4. Velocity profiles beneath Singapore’s seismic network We developed 1D velocity profiles beneath each station in Singapore’s network by jointly inverting stacks of teleseismic receiver functions and surface wave dispersion data. Fig. 1 shows the piercing points at a depth of 5.0 km for all of the observations at each station. There is good azimuthal coverage for all stations with the exception of relatively fewer observations from back azimuths to the southwest. However, because little azimuthal variation was observed in the waveforms (Fig. 2), we stacked data from all back azimuths in order to maximise signal to noise ratios and obtain a representative profile beneath each station. P-wave receiver functions are sensitive to material interfaces with sharp velocity contrasts, and determination of the depths of these interfaces trades-off with layer velocity (Ammon et al., 1990). In order to constrain both interface depths and layer velocities, it is necessary to combine receiver functions with another independent, complimentary data set. Surface wave dispersion data provide constraints on absolute velocities in different depth ranges along great-circle paths between source and receiver. They provide a complimentary data set for receiver functions, allowing more robust estimates of velocity structure than either data set may achieve individually (Juli et al., 2000). 4.1. Surface wave observations There are several techniques available for determining group velocities between source and receiver or between pairs of receivers. The single station method can find the average group velocities along the great circle path between the source and receiver. However, if the source-receiver offset is large, perhaps thousands of kilometres, the technique can only provide crustal information averaged over the long great-circle path. This method may also be affected by source complexity and off-great-circle path propagation (Mokhtar et al., 2001). The two station method and the ambient noise cross correlation method are able to derive group velocities along the great circle path between two stations. Both

methods eliminate source complexity issues and can be used on a highly localised path, the only constraint on the minimum station offset being the maximum period of interest. In order to develop dispersion curves for input into our inversions, we computed inter-station Green’s functions between station BTDF in Singapore and GEOFON station BKNI, a broad-band station equipped with an STS-2 seismometer on the island of Sumatra, by using the ambient noise cross-correlation method (see Fig. 5). The great circle path between these stations is 325 km (see Fig. 5a), a path sufficiently short to provide general crustal information in the vicinity of Singapore. Reliable group velocity measurement requires an inter-station distance of three wavelengths (Ritzwoller et al., 2007). If we assume an average Rayleigh wave velocity of 3.0 km/s, then we can have confidence in dispersion measurements at periods of less than (325 km)/(9 km/s)  36 s. This provides us with a relatively short path, but long enough periods to sample the entire crust. Seven months of data, recorded between January and July of 2010 were used to produce the inter-station Green’s functions. We followed the general processing steps suggested by Ritzwoller et al. (2007), including deconvolving the instrument response, demeaning, de-trending, and spectral whitening before performing the daily cross-correlations. Finally, all daily cross-correlation were stacked in the time-domain. Fig. 5b shows the positive-lag interstation Green’s functions between station BKNI and BTDF, filtered between 0.01 Hz and 0.5 Hz. Group velocities were extracted from the waveforms via a multiple filter analysis (Herrmann, 1987). This technique applies a series of narrow bandpass filters in order to compute the travel time of surface wave groups. We extracted velocities for periods between 4.0 s and 38.0 s. The results of the multiple filter analysis are shown in Fig. 5c. We interpret these results as representative of average group velocities in the vicinity of Singapore, and consider these data as complimentary to our teleseismic receiver function data. 4.2. Inversions For the five stations comprising Singapore’s seismic network, we computed 1D velocity profiles by jointly inverting stacks of teleseismic receiver functions with the surface-wave dispersion

251

K.A. Macpherson et al. / Journal of Asian Earth Sciences 64 (2013) 245–255

(a)

100˚

BTDF

BKNI 0˚



50 km 100˚

(b) RAD

VERT

TAN

0

100

200

time (s)

(c)

4.8

Dispersion between BKNI and BTDF

group velocity (km/s)

4.4 4.0 3.6 3.2 2.8 2.4 2.0 10

20

30

40

period (s) Fig. 5. Station geometry and recordings used in our surface wave group velocity analysis. (a) Map showing the locations of the two-stations used in the cross-correlation analysis. The great-circle path between the two stations is approximately 325.0 km. (b) Positive-lag inter-station Green’s functions derived by cross-correlating 7 months of broad-band data at the two stations. The waveforms have been filtered between 0.01 Hz and 0.5 Hz. (c) Surface wave group velocities extracted from the recordings in (b) by multiple filter analysis. The error bars assume that travel time may be mis-measured by one half filter period (Herrmann, 1987).

data for S-wave velocity (Vs). In order to ease interpretation of the results, here we report P-wave velocities (Vp), simply by assuming pffiffiffiffiffiffiffiffi that V p ¼ 3V s . The inversions were accomplished using the

linearised joint-inversion scheme of Juli et al. (2000). This scheme uses a damped least-squares technique to solve for the model space. For receiver functions, this technique is sensitive to the

252

K.A. Macpherson et al. / Journal of Asian Earth Sciences 64 (2013) 245–255

choice of the initial model, and ideally, the input model will be similar to the actual velocity structure. Because limited velocity model studies have been undertaken for Singapore, we were forced to use the global crustal velocity model CRUST2.0 as the initial model for our inversions (Bassin et al., 2000). An initial 1D model for Singapore was derived by interpolating the CRUST2.0 grid points using a variable tension spline. Initial models for each station were then constructed by modifying the depth to the Moho in order to be consistent with the results of our H  k stacking analysis. We explored the sensitivity of the inversions to the input velocity model by performing a suite of 20 inversions for each stack with slightly different input models. This was accomplished by randomly perturbing the CRUST2.0 profile with a maximum perturbation of 10%. This allows us to gain a degree of confidence in features that persist in many different output models. In order to determine how best to weight the two data sets, we performed many inversion with many different weights and an attempt was made to find a compromise between good waveform fits and relatively simple output models. We found that a 20.0% weight assigned to the surface wave data to be a good compromise. All of the receiver function stacks at each station were inverted with the dispersion data assuming a 20-layer model with 2.0 km thick, flat-lying layers. Using a large number of layers allows the inversion code some freedom in where it places interfaces. 4.3. Inversion results The inversion results are shown in Fig. 6. Each panel shows the output from inverting the receiver function stack for each station jointly with the dispersion data. The input Vp model, modified from CRUST2.0, is shown at the top of each panel in blue, while the Vp models derived from the suite of 20 inversion with slightly different input models are plotted in grey. The best fitting output model in plotted in red. Each panel shows a plot of the stack along with a synthetic computed using the Vp output from the best-fitting inversion at the bottom, and the RMS fit is indicated. The inversion results for station BESC are shown in Fig. 6a. The best fitting inversion produces a waveform with a better than 93% fit. Although the different input models produce some differences in the inversions, the general character of the output velocity profiles is similar. The top-most layer in all models is significantly faster, at around 4.8 km/s, than the input CRUST2.0 model. This top layer is underlain by a fast (Vp  6.1 km/s) layer, beneath which there is a velocity inversion. Below this inversion, velocities gradually increase to a depth of around 28.0 km, at which point the velocity increases relatively sharply to a velocity of Vp  7.9 km/s. We interpret this velocity increase as the Moho, and note that the depth is consistent with the results from our H  k stacking analysis. The inversion results for station BTDF are shown in Fig. 6b. The output models are very similar to those obtained for station BESC, but with an inter-station distance of just over 9.0 km, this is not unexpected. Again, waveform fits are good, with the best fitting model achieving almost a 95% fit. Although there are difference imposed by different input models, persistent features emerge in the output models. Models for BTDF exhibit a pattern similar to BESC of a near-surface layer much faster than CRUST2.0 (Vp  4.5 km/s), beneath which there exists a low velocity zone. The increase to Vp  6.0 km/s beneath the low-velocity zone is sharper for this station than for BESC. The Vp increases slowly from a depth of 10.0 km until around 30.0 km, at which point the velocity increases to Vp  7.6. Although the increase is not as sharp as for station BESC, we interpret the velocity increase at a depth of 30.0 km as the Moho, consistent with our H  k stacking results. Station KAPK exhibits the most complicated waveform of the five stations considered, and consequently, the models resulting

from the inversions are more complex. Unfortunately, this station had the fewest number of teleseismic observations available (30), and we expect this stack to contain more noise than the others. Waveform fits are not as good, with the best-fitting output model achieving just less than an 80% fit. However, some features persist in all of the output models, and also show similarities with other stations. Near-surface velocities are again higher than the CRUST2.0 input model, and there exists a low-velocity zone in the upper crust for all models. Between a depth of 10.0 km and 30 km there are two velocity inversions, separated by a high-velocity (Vp  7.0 km/s) layer. The velocity increases sharply at a depth of 30.0 km to reach Vp  8.0 km/s at a depth of 32.0 km. This appears to be the Moho, and the depth agrees with our H  k stacking results. The output velocity models for station NTU are similar in character to stations BESC and BTDF. The best-fitting output model achieves a greater than 93% fit. Again, the top-most layer in all models is significantly faster, at around 4.4 km/s, than the input CRUST2.0 model. The top layer in underlain by a fast (Vp  5.8 km/s) layer, more pronounced than in the other stations, beneath which there is a velocity inversion. Below this inversion, velocities gradually increase to a depth of around 32.0 km, at which point the velocity increases to Vp  8.0, indicating the presence of the Moho and providing general agreement with our H  k stacking results. The inversion results for station PTK exhibit some similarities with other stations, but some distinct differences at depth. The best-fitting output model achieves a waveform fit of over 96%. In the upper crust, we observe the now familiar pattern of a top-most layer that is significantly faster than CRUST2.0, underlain by a faster layer above a low-velocity zone at a depth of between 4.0 km and 8.0 km. Beneath a depth of 8.0 km, velocities increase sharply at first, and then more slowly to a depth of 20.0 km. Here we observe significant differences from the other stations. There is a low-velocity zone between a depth of 20.0 km and 24.0 km, beneath which velocities increase rapidly, eventually reaching Vp  7.8 km/s. The sharp increase in velocity observed at a depth of 24.0 km is consistent with our H  k stacking results. We speculate that beneath station PTK, there exists a low-velocity zone a the base of the crust, and that the Moho is gradational at this location. This observation is supported by the fact that the Vp/Vs ratio is relatively high for this station, and this may indicate the presence of a low-rigidity layer.

5. Discussion This study inverted stacks of teleseismic receiver functions jointly with surface wave dispersion data for five different station, the greatest inter-station distance being approximately 40.0 km. As might have been expected, many of the crustal features observed in the inversions are visible at all five stations. In particular, there exists a high-velocity layer in nearly every model between a depth of 2.0 km and 4.0 km. Beneath this layer, most models exhibit a velocity inversion, before Vp gradually increases with depth. The minimal amount of azimuthal variation observed in teleseismic waveforms indicates that crustal layers beneath Singapore are flat-lying and have little horizontal variation. However, we do observe differences in near-surface velocities at station NTU relative to the other stations in the network. The lower near-surface velocities observed at NTU are likely a result of the station’s situation above the Jurong formation, which is composed primarily of sandstone, siltstone, shale, and conglomerate, in contrast to the crystalline rock present in central Singapore. Also, while waveforms indicate that crustal layers are flat, the H  k analysis shows that the crust itself thins quite dramatically to the northeast of the

253

K.A. Macpherson et al. / Journal of Asian Earth Sciences 64 (2013) 245–255

Vp (km/s)

(a)

0

2

4

Vp (km/s)

6

(b) 0

8

4

6

depth (km)

20

20

30

30 Station: BESC

Station: BTDF

Input model Output models Synthetic shown below

Input model Output models Synthetic shown below

40

40

stack synthetic fit: 94.8%

stack synthetic fit: 93.5%

−5

0

5

10

15

20

−5

0

5

time (s)

4

15

20

Vp (km/s)

6

8

0

(d) 0

10

10

depth (km)

depth (km)

2

10

time (s)

Vp (km/s)

(c)

8

10

10

depth (km)

2

20

2

4

6

8

20

30

30 Station: KAPK

Station: NTU

Input model Output models Synthetic shown below

Input model Output models Synthetic shown below

40

40

stack synthetic fit: 93.2%

stack synthetic fit: 79.3%

−5

0

5

10

time (s)

15

20

−5

0

5

10

15

20

time (s)

Fig. 6. Inversion results for station (a) BESC, (b) BTDF, (c) KAPK, (d) NTU, and (e) PTK. Each panel shows the results from inverting the receiver function stack for the station jointly with the dispersion data. The top of each panel shows the output velocity models from a suite of 20 inversions obtained by randomly perturbing the starting velocity model with a maximum perturbation of 10%. The unperturbed input velocity model from an interpolated CRUST2.0 is plotted in blue. The best fitting output velocity model is plotted in red, and the synthetic from this model is plotted in red along with the stack in the bottom panel. The fit of the best-fitting synthetic is indicated in the lower panels. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

254

K.A. Macpherson et al. / Journal of Asian Earth Sciences 64 (2013) 245–255

Singapore and GFZ-Potsdam. Deconvolution codes were provided by Professor Charles Ammon, whose on-line receiver function guide was also very helpful. The Computer Programs in Seismology codes and manuals, provided by Saint Louis University Earthquake Center, were used extensively in this study. We used the FuncLab code to do the H  k stacking analysis. Two anonymous reviewers provided helpful comments that improved the final manuscript. This study benefited from helpful discussions with Paul Tapponnier, Emma Hill, Kusno Megawatti, Yanger Walling, and Christina Widiwijayanti. The authors thank Amy Macpherson for providing mapping expertise. Many of the figures were rendered using GMT (Wessel and Smith, 1998). Numerous other free and opensource software packages were utilised in this study, particularly, GNU/Linux.

Vp (km/s)

(e)

2

4

6

8

0

depth (km)

10

20

References

30 Station: PTK Input model Output models Synthetic shown below

40 stack synthetic fit: 96.3%

−5

0

5

10

15

20

time (s) Fig. 6. (continued)

island. The increase in velocity apparent in the inversions for stations PTK at a depth of 24.0 km supports the conclusion of thinner crust in this area. 6. Conclusions This study analysed 697 teleseismic receiver functions and surface wave data derived from the ambient noise cross-correlation of 7 months of broad-band data in order to estimate crustal thickness and velocity structure beneath the five station comprising Singapore’s seismic network. Crustal thickness was constrained beneath each station via an H  k stacking technique. Crustal thickness was found to vary across Singapore from around 32.0 km to 24.0 km, with a pronounced thinning to the northeast. High resolution (twenty 2.0 km-thick layers) 1D velocity models were estimated for each station by jointly inverting stacked teleseismic receiver functions with surface wave group velocity data. These velocity models differ significantly from the CRUST2.0 global velocity model, interpolated for Singapore. We note that non-uniqueness is a problem when inverting geophysical data, and the output shown in this study is just one of many possible estimates. However, by jointly inverting independent data sets, receiver functions and group velocities, we gain a degree of confidence in the results. The velocity models developed in this study add to the available global database of crustal thickness and seismic velocity. At a local level, these models can be combined with site-condition studies and Vs,30 estimates, and should be useful for seismic hazard analysis in Singapore. Acknowledgements This is Earth Observatory of Singapore contribution number 49. Data for this study were provided by the Meteorological Service of

Ammon, C.J., 1991. The isolation of receiver effects from teleseismic P waveforms. Bull. Seismol. Soc. Am. 81, 2504–2510. Ammon, C.J., Zandt, G., 1993. Receiver structure beneath the southern Mojave block, California. Bull. Seismol. Soc. Am. 83 (3), 737–755. Ammon, C.J., Randall, G.E., Zandt, G., 1990. On the non-uniqueness of receiver function inversions. J. Geophys. Res. 95 (B10), 15303–15318. Bassin, C., Laske, G., Masters, G., 2000. The current limits of resolution for surface wave tomography in north America. EOS Trans. AGU 81, F897. Bock, Y., Prawirodirdjo, L., Genrich, J.F., Stevens, C.W., McCaffrey, R., Subarya, C., Puntodewo, S.S.O., Calais, E., 2003. Crustal motion in Indonesia from global positioning system measurements. J. Geophys. Res. 108 (B8), 2367. Crotwell, H.P., Owens, T.J., 2005. Automated receiver function processing. Seismol. Res. Lett. 76 (6), 702–709. Dugda, M.T., Nyblade, A.A., Julia, J., 2007. Thin lithosphere beneath the Ethiopian plateau revealed by a joint inversion of Rayleigh wave group velocities and receiver functions. J. Geophys. Res. 112 (B8), B08305. Herrmann, R.B. (Ed.), 1987. Computer Programs in Seismology, vols. I–VII. Saint Louis University, Saint Louis, Missouri, USA. Juli, J., Ammon, C.J., Herrmann, R.B., Correig, A.M., 2000. Joint inversion of receiver function and surface wave dispersion observations. Geophys. J. Int. 143 (1), 99– 112. Langston, C.A., 1977. Corvallis, Oregon, crustal and upper mantle receiver structure from teleseismic P and S waves. Bull. Seismol. Soc. Am. 67 (3), 713–724. Ligorria, J.P., Ammon, C.J., 1999. Iterative deconvolution and receiver-function estimation. Bull. Seismol. Soc. Am. 89 (5), 1395–1400. Lloyd, S., van der Lee, S., Franca, G.S., Assumpcao, M., Feng, M., 2010. Moho map of south America from receiver functions and surface waves. J. Geophys. Res. 115 (B11), B11315. Macpherson, K.A., Hidayat, D., Goh, S.H., 2012. Receiver function structure beneath four seismic stations in the Sumatra region. J. Asian Earth Sci. 46 (0), 161–176. McCloskey, J., Lange, D., Tilmann, F., Nalbant, S.S., Bell, A.F., Natawidjaja, D.H., Rietbrock, A., 2010. The September 2009 Padang earthquake. Nat. Geosci. 3 (2), 70–71. Megawati, K., Pan, T.-C., 2009. Regional seismic hazard posed by the Mentawai segment of the Sumatran megathrust. Bull. Seismol. Soc. Am. 99 (2A), 566–584. Metcalfe, I., 2000. The BentongRaub suture zone. J. Asian Earth Sci. 18 (6), 691–712. Mokhtar, T.A., Ammon, C.J., Herrmann, R.B., Ghalib, H.A.A., 2001. Surface wave velocities across Arabia. Pure Appl. Geophys. 158, 1425–1444. Natawidjaja, D.H., Triyoso, W., 2007. The Sumatran fault zone – from source to hazard. J. Earthq. Tsunami 1 (1), 21–47. Pan, T.-C., Sun, J., 1996. Historical earthquakes felt in Singapore. Bull. Seismol. Soc. Am. 86 (4), 1173–1178. Park, S.J., Lee, J.M., Ryu, I.-C., 2009. 1d velocity structure beneath broadband seismic stations in the Cretaceous Gyeongsang basin of Korea by receiver function analyses. Tectonophysics 472 (1–4), 158–168 (deep seismic profiling of the continents and their margins). Petersen, M.D., Dewey, J., Hartzell, S., Mueller, C., Harmsen, S., Frankel, A., Rukstales, K., 2004. Probabilistic seismic hazard analysis for Sumatra, Indonesia and across the southern Malaysian Peninsula. Tectonophysics 390 (1–4), 141–158 (strong Ground Motion, Earthquake Hazard and Risk in Alpine-Himalayan and Pacific Region). Rahardjo, H., Aung, K., Leong, E., Rezaur, R., 2004. Characteristics of residual soils in Singapore as formed by weathering. Eng. Geol. 73 (12), 157–169. Ritzwoller, M.H., Barmin, M.P., Levshin, A.L., Lin, F., Moschetti, M.P., Shapiro, N.M., Yang, 2007. Processing seismic ambient noise data to obtain reliable broadband surface wave dispersion measurements. Geophys. J. Int. 169, 1239–1260. Searcy, C.K., Christensen, D.H., Zandt, G., 1996. Velocity structure beneath college station Alaska from receiver functions. Bull. Seismol. Soc. Am. 86 (1A), 232–241. Shaw-Champion, M.E., White, N.J., Jones, S., Priestley, K.F., 2006. Crustal velocity structure of the British Isles: a comparison of receiver functions and wide-angle data. Geophys. J. Int. 166, 795–813. Sieh, K., Natawidjaja, D.H., Meltzner, A.J., Shen, C.-C., Cheng, H., Li, K.-S., Suwargadi, B.W., Galetzka, J., Philibosian, B., Edwards, R.L., 2008. Earthquake supercycles

K.A. Macpherson et al. / Journal of Asian Earth Sciences 64 (2013) 245–255 inferred from sea-level changes recorded in the corals of west Sumatra. Science 322 (5908), 1674–1678. Tomlinson, J.P., Denton, P., Maguire, P.K., Booth, D.C., 2006. Analysis of the crustal velocity structure of the British Isles using teleseismic receiver functions. Geophys. J. Int. 167, 233–237. USGS, 2012. ‘‘Did You Feel It’’. .

255

Walling, M.Y., MacPherson, K.A., Hidayat, D., Megawati, K., 2010. Crustal velocity structure under Singapore inferred from receiver functions study. In: AGU, Fall Meeting. Wessel, P., Smith, W., 1998. New version of the generic mapping tools released. EOS Trans. AGU 79, 579. Zhu, L., Kanamori, H., 2000. Moho depth variation in southern California from teleseismic receiver functions. J. Geophys. Res. 105 (B2), 2969–2980.