Large‐Scale Crustal Structure Beneath Singapore Using Receiver Functions From a Dense Urban Nodal Array

Geophysics has a role to play in the development of “smart cities”. This is particularly true for Singapore, one of the world's most densely populated countries. Imaging of Singapore's subsurface is required to identify geological faults, model shaking from future earthquakes, and provide a framework for underground development. A noninvasive geophysical technique that is well suited for urban areas is passive seismic surveys using compact seismic instruments called nodes. Here we image Singapore's crustal structure using receiver functions generated by a 40‐day deployment of a 88 station nodal array. We generate high resolution receiver functions, despite the noisy environment and short recording time and also create common‐conversion point images. Our results reveal a complex crustal structure, containing multiple discontinuities. Azimuthal variations indicate a distinct change in crustal structure on either side of the postulated Bukit Timah fault, which has implications for seismic hazard.


Introduction
Increasing urbanization creates an ever pressing need to understand the subsurface of our cities. "Smart" cities aim to use technology to improve their liveability, economic, and environmental stability. Geophysics, and seismology in particular, has a role to play in smart city development, for example, through geohazard mitigation and subsurface imaging for underground construction. Technological advancements are allowing us to image the subsurface beneath cities in ways that were previously difficult. In particular, passive seismic surveys using nodes are a noninvasive acquisition advancement that are well suited to urban environments. Here we demonstrate one application of seismic nodes in the world's smartest city (IMD, 2019).
Singapore is a densely populated city nation of 5.6 million people, and better understanding of its seismic hazard and subsurface structure is prescient. Singapore experiences shaking from earthquakes in Sumatra, which is only a few hundred kilometers away (Pan & Sun, 1996). The Mentawai segment of the Sumatran megathrust is the closest to Singapore and is predicted to rupture as a large earthquake in the near future (Sieh et al., 2008). Simulations of future earthquake scenarios in Singapore are needed to ensure it is properly prepared. Imaging the crustal structure beneath Singapore is also required for underground construction, resource development, and identification of geological faults. Local geological faults pose an as yet unassessed seismic hazard. Although Singapore lies on the relatively stable Sunda continental shelf, old fault zones nearby have been reactivated due to stress transfer from the Sumatran megathrust (Shuib et al., 2017).
Several earthquakes, up to M l 5, have occurred within the Malay peninsula during the instrumental time period (ISC, 2019). For example a series of earthquakes up to M l 3.7 in 2007-2009, occurred on a previously unidentified fault 20 km from Kuala Lumpur (Shuib, 2009). It is therefore important to identify significant geological faults, especially in areas of high population density.
Conventional seismic surveys to image crustal structure are not well suited to urban environments. Such surveys typically use a network of large semi-permanent seismometers that use up valuable space and are deployed in relatively quiet environments. On the other hand, dense arrays of nodes are far better suited to urban environments. Seismic nodes can be deployed rapidly, directly into the ground and without any bulky equipment. Nodes-short period autonomous geophones-were originally developed in the energy industry for shallow land imaging with active sources (e.g., Manning et al., 2019, Ourabah et al., 2019. However, several recent studies have demonstrated the utility of nodes for passive seismic imaging on a basin or lithospheric scale (Lin et al., 2013;Liu et al., 2018;Ward et al., 2018;Wang et al., 2019). One drawback of nodes compared to conventional instruments is their short battery life of approximately 1 month. However, the considerable cost saving of nodes compared to conventional instruments permits the acquisition of very dense surveys and receiver density is the principle factor controlling subsurface image quality.
A key passive seismic technique to image crustal structure is the receiver function method (Langston, 1979;Vinnik, 1977). The technique employs body waves generated by teleseismic earthquakes to isolate mode conversions at discontinuities below a receiver. In order to increase the signal above the noise, receiver function imaging typically requires stacking of more than 1 year of teleseismic body waves at a single station (e.g., Macpherson et al., 2013). The long-time needed normally limits the application of the technique to seismic stations that have been active for several years. However, nodal arrays can generate high resolution receiver functions despite their short recording time (Ward et al., 2018). This is because the dense station distribution permits stacking in space, in addition to stacking receiver functions over time. The close proximity of stations also permits higher frequency receiver functions to be used, since coherent mode conversions can be continuously traced between stations.
In this study, we image Singapore's lithospheric structure using teleseismic receiver functions generated from a short-period nodal array. The motivation behind the nodal array was to (1) illuminate the large scale crustal structure beneath Singapore to provide a skeletal framework for future smart city refinement and (2) assess the performance of a nodal array in a noisy urban environment. We show that nodes are fully capable of generating high frequency receiver functions, even in a noisy environment with only 40 days of data. Our results reveal a complex crustal structure, reflecting Singapore's complicated geological history. Azimuthal variations in receiver functions indicate a distinct change in crustal structure on either side of the postulated Bukit Timah fault, which has implications for seismic hazard in the world's third most densely populated country (World Population Review, 2020).

Geology of Singapore
The near-surface geology of Singapore is relatively well known due to an abundance of shallow boreholes and building works (e.g., Leslie et al., 2019). However the geological structure below several hundred meters depth remains almost entirely unknown. The near surface is composed of three principle geological units-variegated arc-derived plutonic rocks and dykes of the Bukit Timah Center granite , very low metamorphic grade sediments of the Jurong Group, and Quaternary shallow marine and terrestrial sediments ; Figure 1a). The Bukit Timah Center granite underlies central Singapore and is the most extensive geological unit, representing repeated igneous intrusions into existing Carboniferous metamorphic basement between ∼285 and 230 Ma (Oliver et al., 2014). Western Singapore and the southern islands archipelago are underlain by the Middle Triassic (242 ± 13 Ma) Jurong Group and Upper Triassic (< 209 ± 12 Ma) Sentosa Group . These groups consist of lightly metamorphosed marine, fluvial, and volumetrically minor pyroclastic rocks that are intensely folded and faulted . Quaternary sediments up to a few hundred meters thick cover much of eastern Singapore and only a small area of northwest Singapore; they exclusively overlie rocks of the Bukit Timah Center (Woon & Yingxin, 2009). Finally, over 20% of the present-day area of Singapore is newly reclaimed land constructed from ∼5 to 10 m of quarried, dredged, or imported sediment.
The nature of the original relationship between the Jurong Group and the Bukit Timah Centre granite is not clear, as the contact does not crop out. While age constraints suggest the Jurong Group could have originally been derived from and deposited on top of exposures of the oldest, westernmost units of the Bukit Timah Center granite, the nearly linear contact between the igneous and sedimentary rocks appears to be a subvertical fault, often referred to as the Bukit Timah Fault Zone (Woon & Yingxin, 2009). The similarities in age between the surficial Jurong Group and the Bukit Timah Center granitoids (which were likely emplaced at significant crustal depth) suggest that at least several kilometers of dip-slip displacement must have occurred along the Bukit Timah Fault Zone, although the time period of such displacement is unknown.

Nodal Array
A nodal array comprising 88, 5-Hz, three-component geophones (Z-land nodes) was deployed across Singapore for a continuous 40-day period, from February to April 2019 ( Figure 1a). Site locations were distributed across Singapore, with a denser profile of 35 nodes located across the boundary between Bukit Timah granite and the Jurong Group. The largest station spacing was 8 km, for a node deployed on a nearby island. The smallest station spacing was 100 m for the nodes deployed across the geological boundary. Site locations included parks, schools, nature reserves, weather stations, roadsides, and industrial sites and so had a range of ambient noise levels. In addition to the nodal array, four permanent seismic station operated by Singapore's Meteorological Service were used.

Methodology
Receiver functions are time series composed of P-to-S converted waves generated at structural boundaries in the Earth beneath a seismometer. Teleseismic waveforms are used to compute receiver functions as they approximate a near-vertically incident plane wave. We construct receiver functions for earthquakes larger than M w 5.5 with an epicentral distance from 30 • to 90 • -a total of 22 earthquakes fitting this criteria occurred during the acquisition time (Figure 1b). We process the raw waveforms by deconvolving the instrument response, removing the mean and linear trend, bandpass filtering from 0.05 to 10 Hz, rotating the horizontal components into the earthquake reference frame and windowing based on the estimated P-wave onset. Removing the instrument response has the effect of boosting low frequencies relative to higher frequencies (supporting information, Figure S1).
Receiver functions are generated by a time-domain iterative deconvolution with a Gaussian low pass filter (Ligorría & Ammon, 1999). The Gaussian low pass filter removes high-frequency noise in the receiver function at the expense of resolution. A Gaussian width of 5 was chosen, which corresponds to a cutoff frequency of approximately 2.5 Hz. This value reduces noise while keeping as much high frequency signal as possible (supporting information, Figure S2). Using a dense nodal array in the United States, Ward et al. (2018) generated receiver functions for a Gaussian value up to 10 (approximately 4.8 Hz). However, in our case, such high Gaussian values generate noisy receiver functions, because there is abundant man-made noise above approximately 3 Hz (Figures S1 and S2). However, our receiver functions still have over one octave more signal compared to conventional surveys (e.g., Macpherson et al., 2013). The dense acquisition means that coherent mode conversions can be traced between stations, which allows more high frequency signal to be utilised, even in an urban environment.
Receiver functions are quality controlled automatically based on their signal-to-noise ratio and subsequently by visual inspection. From the 22 earthquakes that were used to generate receiver functions for the node data, nine earthquakes produced receiver functions of sufficiently high quality. The earthquakes used to compute receiver functions are distributed in three azimuthal bins: 280 • for earthquakes near Fiji, 240 • for earthquake near Japan, and 75 • for one earthquake that occurred in the Indian Ocean. Receiver functions were also generated for four permanent stations in operation in Singapore, composed of one broadband and three short-period instruments. On average, 50 high quality receiver functions were generated per station from 229 suitable earthquakes (M w greater than 5.5 in years 2012, 2013, and 2018).

Results
A stacked receiver function profile across Singapore is shown in Figure 2. Stations within 7 km are projected on to the profile, which strikes across the geological units. Receiver functions with an azimuth between 200 • to 300 • are stacked at each station-this is approximately 90% of all high quality receiver functions and variations within this azimuth range are small. An average of five receiver functions (a minimum of one and maximum of eight) are stacked at each node station and an average of 32 receiver functions are stacked at each permanent station. Receiver functions are also stacked in space. We stack all stations within a radius of 4 km, with the center station doubly weighted. A radius of 4 km is chosen in this case to best increase signal while not smoothing across geological variations. For this cutoff radius, an average of 19 neighboring stations (with a minimum of 0 neighbors and maximum of 39 neighbors) are included in the weighted stack to create the final receiver function at each station. Stacking over time and space increases the signal to noise ratio and coherency of the receiver functions (supporting information, Figure S3).
A clear Ps phase from the Moho is seen at 4 s, matching the previous receiver functions of Macpherson et al. (2013). Additionally, intracrustal discontinuities are present. In eastern Singapore there is a clear positive arrival between the direct P and Ps phases, with negative arrivals preceding and following the peak. This peak arrives at later times towards the center of Singapore and then disappears completely. Conversely, western Singapore appears relatively homogeneous, with negative amplitudes that increase further west.
In order to place the discontinuities at the correct location in depth, we utilize a migration method known as common-conversion point (CCP) stacking (Kosarev et al., 1999;Zhu et al., 2006). The amplitude at each point along the receiver function is back-projected to the conversion point using a background velocity model. We use a 1D velocity model generated by joint inversion of receiver functions and surface waves at a broadband station in the center of Singapore (Macpherson et al., 2013). After back-projection of all receiver functions, the crustal volume is divided into bins, and all amplitudes within a bin are stacked, such that the amplitude of each bin represents the impedance contrast at that location. We use bins of 2 km in length and width and 0.25 km in depth (later interpolated to equal bins of 0.5 km) and smooth the model by a moving average with a width of two bin sizes in length and one in depth. Figure 3 shows the resulting CCP image at several profiles. The Moho is a high amplitude peak at 32 ± 1 km depth. Macpherson et al. (2013) proposed that the depth to the Moho varied at a short wavelength across Singapore (32 to 24 km depth change over ∼25 km); however, we show that the Moho depth is near constant. There are significant intracrustal variations from west to east across Singapore. The negative arrival in western Singapore corresponds to the top of a low velocity zone at 15 km depth. This trough appears to disappear towards the east, where it is replaced by a peak at 19 km depth, marking the top of a high velocity layer. In the east there is an additional low amplitude peak at 7 km depth. It is possible that this peak continues to the west, but it is not clearly separated from the main P arrival. The profiles C-D and E-F appear more uniform, since they are parallel to expected structural trends. The largest change appears to occur in the middle of Singapore, away from surface contacts; we suggest that faults with no surface exposure may be present here. The CCP images show remarkably complex intracrustal lateral variations over a distance of only 40 km.
The majority of teleseismic earthquakes occur to the east of Singapore and so imaging is dominated by energy arriving from the east (azimuths of 240-280 • , where azimuth is defined as the direction from the source to receiver, Figure 1). Fortuitously, however, one suitable earthquake occurred in the Indian ocean during the 10.1029/2020GL087233 nodal deployment and produced high quality receiver functions. Seismic waves from the Indian ocean event travel to Singapore at an azimuth of 75 • , illuminating the far western side of Singapore (Figure 4). This is particularly useful for investigating variations in crustal structure from the Jurong Group to the Bukit Timah granite, across the postulated Bukit Timah fault. Figure 4 shows receiver function profiles generated in the three azimuths groups, using a stacking radius of 1 km. The piercepoints of rays at the base of the crust are shown in Figure 4a for each azimuth group. Waves arriving from an azimuth of 75 • penetrate the crust of the Jurong Group, while the other azimuths sample the crust below the Bukit Timah granite. The receiver functions from an azimuth of 75 • are dramatically different from the other azimuth groups. A peak at 2 s due to an intracrustal boundary appears, with sharp negative arrivals surrounding the peak. On the other hand, the profiles from azimuths 240 • and 280 • show a relatively simple structure with one intracrustal negative arrival as shown in Figure 2. Transverse component receiver functions also show variation with azimuth (supporting information, Figure S4). We propose two explanations for the distinct azimuthal variation in receiver functions shown in Figure 4. The first is a smooth variation in velocity from northeast to west, which would explain the apparent splitting of the negative intracrustal arrival for the westernmost part of the profile at azimuths of 280 • . The second is that there is a near-vertically dipping fault as a sharp velocity boundary, such that the receiver functions from azimuths of 75 • sample the western side, whereas azimuths of 240 • and 280 • sample the eastern side of the fault. The surface expression of this fault would coincide with the postulated Bukit Timah fault. The presence of a significant fault also agrees with geological information showing that the ages of the two units are very similar (Oliver et al., 2014) and therefore that significant tectonic movement must have occurred in order for sedimentary rocks to be adjacent to a granite pluton that formed at the same time.

Discussion and Conclusions
We have generated robust receiver functions at high frequencies using a nodal array in a noisy urban environment. The method relies on (1) high receiver density allowing coherent signals to be identified and stacked, and (2) man-made noise predominantly occurring at frequencies greater than 3 Hz. To create receiver function profiles across Singapore, we first de-noise the data by stacking. We stack both receiver functions from different earthquakes and receiver functions from nearby stations. Assuming a 1D velocity model for the region, we also show for the first time that a temporary nodal array can be used for receiver function imaging using the popular common conversion point method. The CCP method is a ray-based method, and we suggest that future studies might profitably focus on wave-based migration methods (such as reverse time migration, Chen et al., 2005), which may prove superior with high receiver density (Shang et al., 2017). In this study, we do not invert for seismic velocity, since receiver functions have low sensitivity to absolute seismic wave speed. Instead, the receiver functions presented here will be combined with surface waves in a joint inversion for 3D velocity structure.
Our results indicate a complex crustal structure beneath Singapore, reflecting Singapore's complex geological history. In contrast to the previous work of Macpherson et al. (2013), which used around 10 years of data from five permanent stations, we produce more detailed results using only 40 days of data but for 88 stations. For example, we show that the Moho depth is constant across Singapore and does not vary at short wavelength. Our higher frequency results (cutoff frequency of 2.5 Hz vs. 0.75 Hz for Macpherson et al., 2013) also show significant lateral intracrustal variations across Singapore that could not previously be interpreted due to low resolution and low receiver density. In particular, azimuthal variations in receiver functions show distinct crustal structure below the Jurong Group and Bukit Timah granite, indicating the presence of the Bukit Timah fault. Old faults in this area can be reactivated by postseismic deformation following large earthquakes at the nearby Sumatran subduction zone (Yong et al., 2017). Therefore, geological faults in and close-by to Singapore pose a seismic hazard to this densely populated area and warrant further characterization.