Crustal structure of Borneo, Makassar Strait and Sulawesi from ambient noise tomography

SUMMARY

Borneo and Sulawesi are two large islands separated by the Makassar Strait that lie within the complex tectonic setting of central Indonesia. The seismic structure beneath this region is poorly understood due to the limited data availability. In this study, we present Rayleigh wave tomography results that illuminate the underlying crustal structure. Group velocity is retrieved from dispersion analysis of Rayleigh waves extracted from the ambient noise field by cross-correlating long-term recordings from 108 seismic stations over a period of 8 months. We then produce a 3-D shear wave velocity model via a two-stage process in which group velocity maps are computed across a range of periods and then sampled over a dense grid of points to produce pseudo-dispersion curves; these dispersion curves are then separately inverted for 1-D shear wave velocity (Vs), with the resultant models combined and interpolated to form a 3-D model. In this model, we observed up to ± 1.2 km s⁻¹ lateral Vs heterogeneities as a function of depth. Our models illuminate a strong low shear wave velocity (Vs) anomaly at shallow depth (≤ 14 km) and a strong high Vs anomaly at depths of 20–30 km beneath the North Makassar Strait. We inferred the sediment basement and Moho depth from our 3-D Vs model based on iso-velocity constrained by the positive vertical gradient of the Vs models. The broad and deep sedimentary basement at ∼14 ± 2 km depth beneath the North Makassar Strait is floored by a shallow Moho at ∼22 ± 2 km depth, which is the thinnest crust in the study area. To the east of this region, our model reveals a Moho depth of ∼45 ± 2 km beneath Central Sulawesi, the thickest crust in our study area, which suggests crustal thickening since the late Oligocene. Moreover, the presence of high near-surface Vs anomalies with only slight changes of velocity with increasing depth in southwest Borneo close to Schwaner Mountain confirm the existence of a crustal root beneath this region.

1 INTRODUCTION

Borneo and Sulawesi are located in the middle of the Indonesian Archipelago, where they are separated by the relatively narrow (∼200–300 km wide) Makassar Strait. The ongoing collision of the Indo-Australian plate with the Eurasian and Philippine plates has led to a complex tectonic setting in this region, which to date has only seen limited studies into its deeper structure. Borneo is divided among three countries, namely Brunei and Malaysia in the north and Indonesia in the south. The Indonesian territory is called Kalimantan (Fig. 1), which occupies nearly 75 per cent of the total land area. Borneo is located at the front of a laterally extruded block of the Indo-Australian-SE Asian collision zone (Tapponnier et al. 1982; Hall et al. 2008) and is relatively stable, with few earthquakes and a small number of inactive Quaternary volcanoes (Hall et al. 2008). Sulawesi, located to the east of Borneo, has four arms that together form its characteristic K-shape. Sulawesi was assembled by the collision of multiple continental fragments. For example, Central and Southeast Sulawesi, Banggai-Sula, and Buton are believed to have originated in the northern part of the Australian continent (Metcalfe 1988, 1990; Audley-Charles 1991; Surono and Bachri 2002). Borneo and Sulawesi were part of a single land mass in the Late Mesozoic, but were separated during the Cenozoic by rifting in the Makassar Strait (Katili 1978; Hamilton 1979; Hall et al 2009). The Makassar Strait has two basins, the North and South Makassar Basins, which are separated by the Paternoster Platform (PP). The North Basin lies beneath a water depth of ∼2.5 km and is connected to the Celebes Sea (CBS) to the north. The South Basin sits beneath ∼2 km of water and is connected to the shallow shelf area of the Java Sea to the south. On the western side of the North Makassar Basin is the Kutai Basin (KB), the largest and deepest basin in Indonesia with a depth to basement of 15 km (Rose & Hartono 1978); it is also one of the richest hydrocarbon reserves in Indonesia (Hall et al. 2009). The opening of Makassar Strait resulted from rifting during the Middle Eocene (∼48 Ma, Guntoro 1999; Moss & Chambers 1999; Calvert & Hall 2007), which field observations and seismic imaging have shown, produced a highly asymmetrical and extremely wide rift in the Late Eocene (Hall et al. 2009).

On the eastern side of North Makassar Basin is the Palu-Koro Fault, a major sinistral strike-slip fault-oriented north-northwest (NNW)—south-southeast (SSE) across onshore Central Sulawesi (CS) and into the CBS. This fault caused the M_w 7.5 earthquake on 2018 September 28, which unexpectedly produced a large tsunami that inflicted major damage and casualties in the city of Palu and surrounding areas (Omira et al. 2019).

During the last few decades, the increased availability of high-quality seismic data has been instrumental for the application of seismic tomography, for example, Zheng et al. (2021), Greenfield et al. (2022) and Chen et al. (2023). Seismic tomography is ideal for mapping subsurface velocity changes, which can be used to make inferences about the structure and tectonic history of a study region (Widiyantoro & Van Der Hilst 1997; Curtis et al. 1998; Laske et al. 2013; Greenfield et al. 2022). In this study, we use ambient seismic noise tomography (ANT) to image the lithosphere beneath Borneo, Makassar Strait and Sulawesi. The basis of ANT is the extraction of empirical Green's Functions (EGF) from long-term recordings of ambient noise via cross-correlation (CCF, Shapiro & Campilo 2004). The primary advantage of this method is that we can obtain a high-resolution image of crustal structure using a reasonably dense seismic network across a large region, for example, beneath the Australian continent—see Saygin & Kennett (2010, 2012) and Chen et al. (2023), without relying on earthquakes or active sources (Suemoto et al. 2020). Body wave information is typically challenging to extract from ambient noise, so most studies are restricted to the use of surface waves, which can provide good horizontal resolution, but depth resolution is limited. Moreover, the exploitation of group and phase velocity dispersion means that S-wave velocity tends to be the final product of ANT, with P-wave velocity being too poorly constrained by the data to permit retrieval. Ambient noise and other surface wave imaging methods have been widely used in the Indonesian region, for example, Zulfakriza et al. (2014, 2020), Martha et al. (2017), Yudistira et al. (2021), Rosalia et al. (2022) and Pranata et al. (2020), but its application remains limited in Borneo, Makassar Strait and Sulawesi, with the exception of a focused study in northern Borneo (Greenfield et al. 2022). The aim of our study is to recover broad scale crustal structure through the determination of 3-D shear wave velocity in order to contribute to a better understanding of the complex tectonic setting of central Indonesia.

2 DATA

We use continuous data collected from a dense seismic network, with an average station spacing of 58 km, distributed throughout Borneo and Sulawesi (Fig. 2a). This network consists of broadband STS-2 instruments from Indonesia's permanent seismic network (network code IA) operated by the Meteorological, Climatological and Geophysical Agency (BMKG), four of Malaysia's permanent broad-band STS-2 stations (network code MY) operated by the Malaysian Meteorological Service, and two temporary seismic networks deployed in Borneo and Sulawesi between 2019 and 2020 (Greenfield et al. 2018). Together, these networks provide 108 three-component seismic stations distributed across Borneo and Sulawesi. The temporary seismic networks, which feature intermediate broad-band stations (Guralp 6TDs and Nanometric Trillium Compacts) are especially important for improving the ray coverage of Borneo, which has fewer permanent seismic stations. For the CCF, we used eight months of continuous data recorded on the vertical component, which allowed us to retrieve good quality Rayleigh wave Green's functions.

3 MATERIALS AND METHODS

ANT proceeds by first calculating the ambient noise CCFs between station pairs. These can be considered to be the EGFs corresponding to Earth structure between the stations (Weaver & Lobkis 2001; Fichtner et al. 2008). To improve the EGF signal-to-noise ratio (SNR), we stack daily CCFs and then measure the Rayleigh wave group velocity dispersion by applying wavelet domain analysis (see Section 3.1). Using the dispersion measurements, we then construct a group velocity map for each period and extract pseudo-dispersion curves on a 1° × 1° grid. These pseudo-dispersion curves are then inverted for 1-D shear wave velocity (Vs) at each gridpoint, with a final 3-D shear wave velocity model of the crust obtained by stitching all the separate 1-D models together via interpolation.

3.1 Dispersion curve measurement

We use a modified version of the Python code NoisePy (Jiang et al. 2020) to compute the dispersion curve between station pairs. The code was designed specifically for large-scale ambient noise seismology and provides routines for downloading data, processing and dispersion analysis (Jiang et al. 2020). NoisePy implements both the wavelet domain scheme as in Fichtner et al. (2008) and time–frequency analysis (FTAN; e.g. Bensen et al. 2007) for the group velocity estimation (Jiang et al. 2020). As discussed in Fichtner et al. (2008), the wavelet domain approach is equivalent to the FTAN in the case that the width of the sliding window does not depend on time or frequency. We applied a standard ambient noise processing workflow (Bensen et al. 2007) as follows: the data were resampled into daily segments at a sampling rate of 2 Hz, before the instrument response, mean and trend were removed from the data; a bandpass filter between 0.02 and 1 Hz was then applied. This sampling rate was chosen to reduce computation time, since we only consider the ambient noise signal in the period band 8–50 s. We used a combination of spectral whitening with a running mean average (rma) across 10 sample points and waveform CCF between station S and R in the Fourier domain defined by

where |${U}_\mathrm{ S}$| (f) is the Fourier transform of the seismogram at virtual source station S, * denotes the complex conjugate, |${U}_\mathrm{ R}$| (f) is the Fourier transform of the seismogram at virtual receiver station R, | | denotes the absolute value and {} represents a smoothing of the spectrum using rma (Jiang et al. 2020). The CCFs were computed for each station pair with a lag-time of ± 800 s. After being transformed back into the time domain, we combined all the daily CCFs from each station pair using a linear stacking method. Fig. 2(b) shows the stacked CCFs extracted from all interstation pairs, which reveals very clear EGFs across the entire distance range.

The dispersion plot for each interstations pair was obtained using the wavelet domain analysis tool included in NoisePy. It uses a Morlet wavelet with default parameters for frequency spacing, scale of the wavelet and number of scales as defined in the pycwt package (https://pycwt.readthedocs.io, Jiang et al. 2020). To provide high-quality dispersion data, we first identified and picked the best EGFs by using the SNR calculated as the ratio of the maximum absolute amplitude of the signal window and the root mean square (RMS) of a noise window as suggested by Ritzwoller & Feng (2019). The signal window is taken between the minimum time, |${t}_{\mathrm{ min}}$|⁠, and the maximum time, |${t}_{\mathrm{ max}}$|⁠, where |${t}_{\mathrm{ min}}$| and |${t}_{\mathrm{ max}}$| are defined by 4.5 and 1.5 km s⁻¹, respectively, for a given interstation distance. The noise window is taken between |${t}_{\mathrm{ max}}$| and the end of the time-series. We used an SNR threshold of five, as used in a previous ANT study (Zheng et al. 2021). This threshold of five was chosen to satisfy the cut-off between the number and the quality of EGF data. Using this threshold, the number of EGFs is reduced by about 40 per cent from the total EGFs obtained from all the station pairs; however, this step leaves only high-quality EGFs, with in this case 1700 subsequently used in the dispersion curve analysis. Finally, the dispersion curves were extracted from the maximum amplitude of the group velocity measurements at each period for each station pair. The quality of dispersion data obtained in these measurements were then assessed, as discussed in Section 3.2.

3.2 Dispersion measurement uncertainties

We perform quality checks on the interstation dispersion measurements prior to group velocity analysis in order to weight the data in the inversion step. Here, we examine the spatial and temporal variability of the dispersion curves, for example, Bensen et al. (2007), as summarized in Figs 3 and 4. Fig. 3 shows dispersion curves extracted from four interstation paths, all with similar geometries, suggestive of similar subsurface properties and standard deviations ranging from 0.01 to 0.1 km s⁻¹ (shown by the blue-shaded area) across all periods. However, this assessment is only valid when a cluster of stations subtends a small angle to a relatively distant station, which only occurs for a subset of measurements (Bensen et al. 2007). Thus, we estimate the dispersion uncertainties quantitatively using the temporal variability of dispersion measurements. To do so, we create three separate stacks that comprise the first and second four months of stacked CCFs and the full eight months of stacked CCFs (Fig. 4a). Fig. 4(b) shows that the dispersion curve generated from the second four months of stacked CCFs, which has lower SNR (with an SNR of 4.6, as shown by the red line), leads to a significant deviation from the full eight month dispersion curve. However, the dispersion curve obtained from higher SNR data recorded in the first four months (with an SNR of 7.2, as shown by the blue line) is observed to be quite similar to the curve obtained from the eight months of stacked CCFs. This shows that the SNR check is suitable for assessing the quality of the stacked CCFs or EGFs prior to dispersion curve analysis. To obtain quantitative uncertainty estimates of the dispersion measurements, we compute the standard deviation of each period extending from 8 to 50 s with a sampling of 1 s. Fig. 4(c) shows an example of the dispersion data uncertainties calculated from the station air of PKKI and TTSI. We can observe that the standard deviation tends to increase with increasing period. Using this approach, we can calculate the dispersion data uncertainties for all station pairs. Fig. 4(d) shows the average standard deviations obtained from all interstation pairs.

3.3 Group velocity anomalies

We used the Fast-Marching Surface Tomography code (Rawlinson & Sambridge 2005) to produce group velocity maps as a function of period across the study area. In this step, we only consider the dispersion data ranging from 8 to 45 s with an increment of 1 s to produce group velocity maps, since our data have larger uncertainties at periods above 45 s. We estimated the traveltime between station pairs by picking the corresponding interstation dispersion curve as a function of period and converting path-average group velocity to time. To obtain a starting model for each period, we computed the uniform velocity field that best satisfied the data set. The resulting starting model exhibits an approximately linear increase in velocity with period, from 2.7 km s⁻¹ at 8 s period to 3.6 km s⁻¹ at 45 s period. The inversion was conducted for the target area spanning Borneo, the Makassar Strait and Sulawesi. To test the ability of our data set to recover features of a certain size, evaluate which regions are well-recovered and choose optimal grid spacing, we perform checkerboard tests at all periods before running the inversion using real data. Fig. 5 shows example checkerboard tests at 10, 20, 30 and 40 s period using the station-pair distribution shown in Fig. 5 (bottom). We test checkerboard anomalies with 0.5° (Fig. 3, top) and 1° grid spacing (Fig. 5, middle). From the results of these tests, the 1° checkerboard anomalies are recovered across most of the study area, while the 0.5° checkerboard anomalies were only recovered across Sulawesi. Consequently, we used a 1° grid for the inversion of our observed data set, while acknowledging that some finer scale structure in Sulawesi that may be constrained by the data will not be recovered. We perform trade-off tests to determine appropriate damping and smoothing parameters for the real-data inversion. A damping value of 15 was found to provide a suitable balance between satisfying the data and producing models that do not differ significantly from the starting model. A smoothing value of 3 was chosen to provide smooth models that still exhibit a satisfactory reduction in data variance (see Fig. S3, Supporting Information). Although there is some variation in path coverage between different periods, the chosen values of damping and smoothing were found to be suitable across the period range that was used.

Following the procedure outlined by Barmin et al. (2001), we identify and remove the traveltime outliers with an absolute residual between the observed and predicted traveltime greater than two times the standard deviation. Our final Rayleigh wave group velocity model for each period was obtained after 12 iterations and the average RMS traveltime misfit reduction is about 90 per cent. Fig. 6(a) shows the number of paths at each period after removing the traveltime outliers. Examples of the resultant group velocity map at each period are shown in Fig. 7. From these figures, prominent low-velocity anomalies are clearly revealed at periods between 10 to 20 s. Furthermore, prominent high-velocity anomalies at longer periods (periods ≥ 20 s) are clearly shown, that is, in the North Makassar Strait (NMS), PP and south Borneo.

In addition, we also performed the inversion using a finer grid of 0.5°. The group velocity maps produced using the finer grid (0.5°) are shown in Fig. S4 (Supporting Information), which reveals some of the extra details that can be recovered. However, our primary interest is to interpret similar-scale features across the entire study region, and therefore prefer to use the coarser grid in the inversion to produce a single model that is everywhere well resolved.

3.4 Shear wave velocity inversion

To obtain a 1-D depth-dependent Vs model, we extract pseudo-dispersion curves on a 1° grid spacing covering the study area. 116 of 340 gridpoints were selected based on well-recovered areas from checkerboard tests (see Fig. 5). The pseudo-dispersion curves at the selected gridpoints (Fig. 6b) were then used as inputs for shear wave velocity inversions. Strong low-velocity anomalies are clearly shown in the period ranges of 10–15 s (Fig. 6b). These anomalies are robust because the group velocity maps (Fig. 7) are well constrained across these periods, as explained in Section 3.3, and also contain distinct low-velocity anomalies. These curves are then inverted using an iterative least-squares technique contained in the package surf96 (Herrmann & Ammon 2002; Herrmann 2013). According to Herrmann & Ammon (2002), the difference in observed, |${U}_{\mathrm{ obs}}$| and calculated group velocity, |${U}_{\mathrm{ cal}}$|⁠, is given by:

where N is the number of layers, |$\Delta {V}_{s_i}$| is updated model from the prior and the partial derivative of group velocity with respect to shear wave velocity are derived by following Rodi et al. (1975):

where c is the phase velocity and T is the period. When group velocities are computed, the phase velocities are also computed at the two periods (1 ± h)T. The required phase velocity partial derivative can be calculated as follows:

The parameter h is given a value of 0.005. By introducing the regularization parameters of weight, W, and damping, |$\sigma $|, to eq. (2), the updated model from the prior, |$\Delta {V}_s$|⁠, can be retrieved using the generalized inverse as follows:

where the U, |$\lambda $| and V are decomposition matrices of group velocity partial derivatives, |$\partial U/\partial {V}_s$|⁠.

The prior model in each case was parametrized as a depth dependent seismic structure from the surface to 60 km depth, with homogeneous layers 2 km thick from the surface to 40 km depth and 5 km thick below 40 km depth, since our data do not have much sensitivity below 40 km. We use a monotonic Vs constraint as a prior model by setting Vs to linearly increase with increasing depth (see blue lines in Figs 8b and d) as also done in, for example, Feng & Ritzwoller (2019), Feng (2021a, b) and Feng & Diaz (2023). The Vs in the crust increases from 3 km s⁻¹ at 2 km depth to 4.0 km s⁻¹ at 36 km depth, which we assume to be the Moho. Below the Moho, Vs increases from 4.0 km s⁻¹ at 36 km depth to 4.2 km s⁻¹ at 60 km depth. We use Brocher (2005) regression to estimate the prior Vp model and the density of each layer.

To estimate the uncertainties and to obtain an appropriate damping value, we perform an inversion of the synthetic group velocity (plotted as black dots with error bars in Fig. 8a) computed from the synthetic Vs model (plotted as green line in Fig. 8b). The synthetic Vs model was created using a simple 1-D Vs model obtained from ak135 (Kennet et al. 1995). We chose an optimum damping value to obtain the best 1-D Vs model (see Fig. S5, Supporting Information). We also perturb the synthetic data 1000 times by adding Gaussian noise with a standard deviation ranging from 0.14 km s⁻¹ at 8 s period to 0.24 km s⁻¹ at 45 s period, and then invert each synthetic pseudo-dispersion curve using the same parametrization as before to obtain the Vs model. As shown in Fig. 8(b), the final 1-D Vs model (red line) is similar to the average Vs model taken from the ensemble of 1000 models, indicating that our results are robust. In addition, the grey area shows the standard error of these models, which varies from 0.03 km s⁻¹ at shallow depths to 0.1 km s⁻¹ at greater depths and below the Moho. However, as shown in Fig. 8(b), the velocity jump across the mid-crust (Conrad discontinuity) and the Moho, as shown in the synthetic model, cannot be resolved by the inversions, even though the dispersion misfits are relatively small, as shown in Fig. 8(a). Previous studies of the Moho using seismic surface waves reported that Vs inversions suffer from a trade-off between the Moho depth and the uppermost mantle Vs (Lebedev et al. 2013). Here, we perform an inversion of the synthetic group velocity that explicitly includes the Moho by applying a small weight in the neighbourhood of the Moho to permit significant change in the Vs model. Both the depth and velocity contrast at the discontinuity can be subsequently modified during the inversion. Similar to Figs 8(a)–(c) show the synthetic (black dot with error bar) and calculated dispersion curves (red line) while Fig. 8(d) shows the Vs model (red line) produced by the inversion of synthetic group velocity with the reduced regularization constraint at 30 km depth. The velocity jump at the given Moho depth can therefore be recovered by this inversion scheme and the error reduced (see the grey lines in Fig. 8d).

For the real data inversions, we use Moho depth information from the CRUST1.0 model (Laske et al. 2013) at 116 selected gridpoints covering the study area. Fig. 9 shows an example of the Vs models obtained from the inversion using Moho depth constraints of CRUST1.0 and Vs gradients. The gradient is calculated from each Vs solution model. The impulsive gradients at the Moho are clearly shown at points A and C, indicating velocity jumps at this depth. However, at coordinate point B, a weak gradient occurs at the Moho, while a rapid velocity change occurs at shallower depths (around 15 km), which may indicate that the real Moho could be shallower than the CRUST1.0 Moho at this location.

3.5 Crustal interface measurements

Once the 1-D Vs models have been obtained—as described above—they are combined into a pseudo-3-D model using a smooth interpolator. To delineate the primary interfaces beneath a study area, Chen et al. (2023) suggested using the gradient maxima to determine the sediment–basement interface and the Moho depth. This method assumes a ‘typical’ crust that features a distinct sediment–basement interface and Moho, as seen in Fig. 9 at points A and C. Since we use Moho depths from the CRUST1.0 model, the location of velocity jumps is variable in depth. However, it is difficult to find the gradient maxima for gradients that are not pronounced, for example, at point B. As a consequence, we use a different approach based on prescribed velocity contours—in this case, a value of 3.0 ± 0.1 km s⁻¹ for the sediment–basement interface and 3.9 ± 0.1 km s⁻¹ for the Moho provided the velocity gradient is positive (velocity increases with depth). Mooney et al. (1998) produced a P-wave velocity (Vp) classification for sediment layers: 2.0–3.0 km s⁻¹ for soft sedimentary layers on continents and 4.0–5.3 km s⁻¹ for hard sedimentary layers. Using Brocher (2005) regression, the equivalent for i-wave velocity (Vs) is 0.6–1.4 km s⁻¹ for soft sediments, and 2.3–3.2 km s⁻¹ for hard sediments. Based on these results, we thus assume that a Vs of 3 km s⁻¹ is the velocity of the sediment–basement interface. For the Moho, we assume that a Vs of 3.9 km s⁻¹ is the highest velocity in the crust—based on the ak135 model of Kennet et al. (1995)—and can therefore be used as a proxy for the Moho. Both the basement interface and Moho can be clearly recovered in ∼98 per cent of our 1-D Vs models.

4 RESULTS

We show horizontal slices through the recovered Vs model between 8 and 40 km depth in Fig. 10 and the corresponding errors. From our horizontal slices, we observe up to ± 1.2 km s⁻¹ of lateral Vs heterogeneities as a function of depth. The results illuminate significant low and high Vs anomalies across the region. We observe a strong low near-surface Vs anomaly zone (LSVZ) at shallow depths (≤ 14 km) beneath NMS as delineated by the red dashed lines in Fig. 10. Other LSVZs are also observed beneath Tomini Bay (TB), Bone Gulf (BG, including the South Arm of Sulawesi, SS), the eastern part of the Southeast Arm of Sulawesi (SES) and Sarawak (northwest Borneo). The LSVZ beneath the North Makassar Strait is strongest at 8 km depth but can be observed down to a depth of 14 km (Fig. 10). This prominent LSVZ extends from KB, East Borneo, to the west coast of Sulawesi in an east–west direction and spans a distance of ∼300 km; it also extends in a north–south direction from PP to Mangkalihat Peninsula (MP), a distance of ∼400 km. This is the broadest and deepest LSVZ observed in this region and also has the lowest Vs of ∼2.2 km s⁻¹ observed at 8 km depth. In contrast, we observe a high near-surface velocity zone (HSVZ) at 8 km depth (Fig. 10) in Southwest Borneo and CS. Moreover, we also observe high deep Vs anomaly zones (HDVZs) at depths of 30 km (as indicated by the blue dashed lines in Fig. 10). These HDVZs tend to be co-located with LSVZs at shallower depth (∼8 km), that is, in the North Makassar Strait, TB, SS and the SES, indicating rapid velocity changes with increasing depth. The errors in the model are highest in the North Makassar Strait where strong Vs anomalies are found.

Fig. 11 shows six vertical profiles of our Vs model through Borneo, the Makassar Strait and Sulawesi, along with their associated error. The black solid lines denote the iso-velocity obtained in this study, while the red and black dashed lines reveal the sediment basement and Moho respectively. The iso-velocities of 3.8–4.0 km s⁻¹ are largely in the vicinity of the CRUST1.0 Moho model (plotted as black dashed lines); hence, it appears reasonable to assume that the iso-velocity contour of 3.9 ± 0.1 km s⁻¹ is appropriate for approximately locating the Moho from our Vs models. However, in some regions, our Moho model deviates significantly from CRUST1.0, for example, in the NMS and in CS to the southwest of TB—see Fig. 11 cross-section B-B′ and E-E′. To delineate the geometry of the crustal interfaces, we compute the 2-D and 3-D depth contour profiles that correspond to the sediment basement and Moho depth as shown in Fig. 12. From this figure, we can infer the depth variation of the sediment basement (Fig. 12a) and the Moho (Fig. 12b) across the study area. The average depth-to-basement estimated from our Vs model is ∼6 km, with a minimum of 2 km and maximum of 14 km, while the average Moho depth is ∼36 km, with a minimum of 24 km and maximum of 45 km.

From the vertical Vs profiles, the crustal structure beneath the study area can be readily observed. The sediment, the crust and the lithospheric mantle layer are clearly depicted in our model, which illuminates strong LSVZs in the upper crustal layer and in the sedimentary layer, that is, in the neighbourhood of the NMS, including the KB in the west and the MP in the north, as well as TB and BG. This suggests the likely presence of sedimentary material with a Vs ≤ 3.0 km s⁻¹ (Mooney et al. 1998) beneath the NMS, the SS, the BG and TB.

5 DISCUSSION

The widespread LSVZs beneath Sarawak Basin (SB), NMS, SS, BG, TB and the eastern margin of the SES indicates the presence of extensive sedimentary basins across the study area. The broad and thick LSVZs observed by this study shows good agreement with previous geology and geophysical studies. For example, the LSVZ observed beneath Sarawak confirms the presence of SB, which has been discussed previously by, for example, Hamilton (1979) and Moss (1998). The SB incorporates the Rajang–Crocker Group of Sarawak, which is interpreted as remnant oceanic basin fill created by deposition of an extensive Santonian–Palaeocene (∼86–59 Ma) submarine fan on a passive margin or subduction accretionary complexes (Hamilton 1979; Moss 1998). The other example, the LSVZs in SS and BG, confirms the findings of previous geology studies that reported the BG as a Neogene interarm basin (∼23–4 Ma) dominated by extension and underlain by a variety of pre-Neogene rocks. Proven petroleum plays reported in the BG include gas fields in the East Sengkang Basin in the SS (Camplin & Hall 2014). Our models also show widespread LSVZs in the eastern margin of the SES. Previous geology studies reported oil shows and discoveries in and around the island of Buton at the tip of the SE Arm (Camplin & Hall 2014). Oil seeps are reported offshore near Kolaka (Thompson et al. 1991). Furthermore, the LSVZ in the TB remains enigmatic and less studied; however, based on marine magnetic studies, the TB is underlain by an oceanic-like crust with suspected thermal extension, thinning and differential subsidence in the basement (Kusnida et al. 2009). According to Hall (2011), the additional extension which contributed to subsidence in this region was caused by subduction rollback at the North Sulawesi Trench.

The most broad and deep sedimentary basement extending to a depth of ∼14 ± 2 km can be observed in the NMS (Fig. 12a). The NMS Basin, which includes the KB, has an asymmetrical shape and extends laterally by ∼300 km in the east–west direction and ∼400 km in the north–south direction. This basin is bounded by Central Borneo Mountain (CBM) in the west, CS in the east, PP in the south, and the southwestern margin of the CBS in the northwest. The widespread LSVZ (Vs < 3.0 km s⁻¹) is underlain by a deep sediment–basement interface at ∼14 ± 2 km depth (Fig. 12a) and is co-located with a shallowing of the Moho (Fig. 12b). This indicates the presence of thin crust beneath the NMS, likely resulting from rifting during the Middle Eocene (∼48 Ma). The thinnest crust beneath the NMS is underlain by a Moho at ∼22 ± 2 km depth, which is slightly thicker than the 20 km elastic thickness estimated by Cloke et al. (1999), who assumed that the crust beneath the NMS was 48 Myr-old oceanic crust. Previous gravity modelling (Guntoro 1999) showed that the lowest free-air anomaly of −40 mGal is in the middle of the NMS basin and increases to + 70 mGal towards the continental shelves of Borneo and Sulawesi. Guntoro (1999) postulated that there is a change in the thickness of the crust, due to deformation by extensional thinning. Furthermore, most of Sulawesi is characterized by crust rarely thicker than 35 km (Fig. 12b). However, our models show a significant increase in crustal thickness from ∼22 ± 2 km beneath the NMS basin depocentre to ∼45 ± 2 km beneath CS. The crust to the west of the Makassar Strait is not as thick (30–35 km); this would suggest that the change in crustal thickness is not related to the extension of the NMS but is instead a pre-existing feature.

As shown in Fig. 12, the deep Moho beneath CS is caused by low deeper Vs anomalies (LDVZ) that extend into the lithospheric mantle. To understand what this region of low velocities could represent, we briefly describe how Sulawesi was assembled. As discussed earlier, the island of Sulawesi is a juxtaposition of continental fragments from both Laurasia and Gondwana. The process that formed Sulawesi started in the Oligocene (∼30 Mya) with the continental fragment of Banggai–Sula drifting westwards toward Sulawesi, driven by westward subduction beneath the East Sulawesi Ophiolite (which comprises ES and SES). At this time, East Sulawesi and West Sulawesi (comprising SS and the western part of CS) were not connected (Wilson & Moss 1999). During the late Oligocene, the East Sulawesi Ophiolite was accreted onto West Sulawesi. This terminated the rifting of the Makassar Stait (Guntoro 1999) and sediments began to accumulate in the subsiding basins. The supply of sediments was increased by the collision of East and West Sulawesi, which also caused uplift. During the Miocene or earliest Pliocene, the Banggai–Sula continental fragment was accreted onto Eastern Sulawesi. A transpressional regime between 5 and 2 Ma then generated further uplift in CS estimated to be 200–700 m Ma⁻¹ (Wilson & Moss 1999). The vertical extent of LDVZ explains the thickened crust beneath CS and the rapid uplift the region has experienced since the late Oligocene. The rapid uplift of CS and the tropical latitudes have caused very high erosion rates and by extension the thick sedimentary package infilling the NMS and imaged in our velocity model. Due to the limit of vertical resolution at greater depths in our model, we cannot say whether the subducted oceanic crust is still attached. However, there is no indication in global, for example, LLNL-G3Dv3 (Simmons et al. 2012), or regional (Zenonos et al. 2019; Wehner et al. 2022a,b) tomography models of a subducting slab in this region. It is likely therefore that the slab has detached, further increasing uplift rates in CS (van Hunen & Allen 2011). The uplift and gravitational collapse of Sulawesi is still ongoing and driven by the transpressional stress regime and internal deformation along major strike-slip faults. The broad westward vergent fold and thrust belt in offshore West Sulawesi has been documented by Sandwell & Smith (2009) and Puspita et al. (2005). The fold and thrust belt is evidence of shelf failures caused by major uplift and high exhumation (Hall 2011) and are still active and capable of producing large magnitude earthquakes (Greenfield et al. 2021; Meilano et al. 2023). Another example of a vertically extensive LDVZ, can be found beneath the CBM, which explains the thickened crust in the region. Previous geology studies, for example, Hall et al. (2008) and Hall (2011), suggested that CBM crust is mechanically weak and has experienced shortening and uplift in the Early Miocene as a result of compressional forces applied by adjacent regions of stronger crust.

The presence of a shallow Moho or thin crust suggests the lack of a crustal root as discussed by Szanyi et al. (2021) and Simonová et al. (2019). A previous ANT study in the Eastern Alps presented by Szanyi et al. (2021) found that higher velocities were observed at shallower depths beneath basin areas, suggesting the presence of thin crust and hence the absence of a crustal root. Our results beneath NMS, TB and BG show similar crustal characteristics, suggesting the lack of a crustal root beneath these locations (see Fig. 11). In particular, we observe low near-surface Vs anomalies (Vs < 3.0 km s⁻¹) pointing to the presence of extensive sedimentary basins, which are underlain by a strong positive Vs gradient, suggesting the presence of a shallow Moho and hence a thin crust. Moreover, Szanyi et al. (2021) also suggest that the ∼homogeneous structure embedded in the higher velocity mantle material, which is characterized by slight changes of velocity with depth, indicates the presence of a crustal root. From our results, the presence of HSVZs in southwest Borneo, close to Schwaner Mountain (SM), confirms the existence of a crustal root, since this region is also characterized by slight changes in velocity with increasing depth (see Fig. 11). Borneo has a Palaeozoic continental core in its west to southwest regions, surrounded by ophiolitic, island arc and microcontinental crust accreted during the Mesozoic between 252 and 72 Ma (Hamilton 1979), which is interpreted as a block derived from the Gondwana margin of NW Australia (Hall et al. 2009; Hall 2012). The SM itself is characterized by Cretaceous granitoid plutons intruding the older metamorphic rocks (Hall et al. 2008).

6 CONCLUSIONS

We determined the crustal shear wave velocity structure of Borneo, Makassar Strait and Sulawesi using Rayleigh wave group velocity ANT. The primary conclusions of our study are as follows:

1. We observed lateral Vs heterogeneities of up to ± 1.2 km s⁻¹ in the crust and uppermost mantle. The major sedimentary basins in the study area are well depicted from our Vs models, for example, in the NMS, TB, BG (including the SS), the eastern part of the SES and in Sarawak (northwest Borneo).

2. Our models illuminate a strong low shear wave velocity (Vs) anomaly at shallow depth (≤ 14 km) and a strong high Vs anomaly at depths of 22–30 km beneath NMS.

3. We inferred the depth of the sediment basement and Moho from our 3-D Vs model using iso-velocity contours along with a requirement for a positive vertical gradient in Vs. The broad and deep sedimentary basement at ∼14 ± 2 km depth beneath the NMS is floored by a shallow Moho at ∼22 ± 2 km depth, which is the thinnest crust in the study area.

4. Our models reveal a significant change in crustal thickness to the east, from ∼22 ± 2 km beneath the NMS basin depocentre to ∼45 ± 2 km beneath CS, which likely reflects crustal thickening since the late Oligocene.

5. The presence of high near-surface Vs anomalies and only slight changes in velocity with increasing depth in southwest Borneo close to SM confirms the presence of a crustal root beneath this region.

Acknowledgement

We thank two anonymous-peer-reviewers and the editor Prof Huajian Yao, for constructive comments that greatly improved the quality of this manuscript. We are grateful to Meteorological Climatological and Geophysical Agency (BMKG) and Malaysian Meteorological Service (Ministry of Science, Tech and Innovation, MOSTI-MMD-MET) for providing raw seismic data from Borneo and Sulawesi that are available at https://geof.bmkg.go.id/webdc3 and FDSN Web Services at IRIS (http://service.iris.edu). This work is supported by BMKG and Hibah Riset Institut Teknologi Bandung (ITB, https://www.itb.ac.id) 2022, awarded to ADN, KEMENDIKBUDRISTEK (Penelitian Dasar Kemitraan), contract no: 083/E5/PG.02.00.PT/2022, and matching funds from NFIS awarded to SW and British Council Newton Impact fund G107511. Most of the figures were plotted using the Generic Mapping Tools (Wessel et al. 2013) and most of the color figures were plotted using scientific colour maps from Crameri (2018).

DATA AVAILABILITY

The data sets used and analysed during the current study are :

1. The waveform data of the Malaysian Seismic Network (https://www.fdsn.org/networks/detail/MY) with station code; KKM, KSM, LDM and SBM are available from the Incorporated Research Institutions for Seismology (IRIS) Data Services (http://ds.iris.edu/mda/MY).

2. The waveform data from Indonesian Seismic Network in Borneo and Sulawesi is restricted access due to the institutional regulations of BMKG. The data are available at https://geof.bmkg.go.id/webdc3, with permission from BMKG required to access the data.

3. Several IA stations are also available in the GEOFON Program (GEOFON), (https://www.fdsn.org/networks/detail/IA)

4. The waveform data from the temporary network deployed in Borneo (https://doi.org/10.7914/SN/9G_2018) and Sulawesi during 2019–2020 are restricted; requests for these data should be sent to TG or NR for Borneo data and to ZZ or SW for Sulawesi data.

5. The 3-D shear wave velocity (3-D Vs) model produced in this study and the data underlying this paper are available at https://doi.org/10.5281/zenodo.10662743.

References