Post‐Subduction Tectonics of Sabah, Northern Borneo, Inferred From Surface Wave Tomography

We use two‐plane‐wave tomography with a dense network of seismic stations across Sabah, northern Borneo, to image the shear wave velocity structure of the crust and upper mantle. Our model is used to estimate crustal thickness and the depth of the lithosphere‐asthenosphere boundary (LAB) beneath the region. Calculated crustal thickness ranges between 25 and 55 km and suggests extension in a NW‐SE direction, presumably due to back‐arc processes associated with subduction of the Celebes Sea. We estimate the β‐factor to be 1.3–2, well below the initiation of seafloor spreading. The LAB is, on average, at a depth of 100 km, which is inconsistent with models that ascribe Neogene uplift to wholescale removal of the mantle lithosphere. Instead, beneath a region of Plio‐Pleistocene volcanism in the southeast, we image a region 50–100 km across where the lithosphere has thinned to <50 km, supporting recent suggestions of lower lithospheric removal through a Rayleigh‐Taylor instability.

advocated the Celebes Sea (CS) slab breaking-off and/or lithospheric delamination beneath Mt. Kinabalu as causing uplift. An alternative model focused on results from global tomography, inversion of river profiles, and Neogene magmatism over much of Borneo (Roberts et al., 2018). The authors concluded that much of the uplift through the Neogene could be attributed to the emplacement of a hot asthenosphere layer and subsequent melting and intrusion.
A common feature of all these models is that they rely on field observations of surface features, but deeper, geophysical images of the subsurface, especially of the thicknesses of the crust and lithosphere is lacking. In this study, we address this issue by using a dense network of seismic stations deployed across Sabah to construct a shear wave velocity (Vsv) model using two-plane-wave tomography. We show that lithospheric thinning has occurred, consistent with extension across north and central Sabah, yet the presence of mantle lithosphere is still present beneath much of Sabah, which cannot be reconciled with the idea of wholesale replacement by hot asthenosphere. Instead, we image a localized region where the lithosphere appears to have been removed beneath the youngest volcanic province.

Data Set
We use seismic waveforms recorded by the nBOSS array in Sabah (network YC, Figure 1; Rawlinson, 2018) as our base data set. This is supplemented by stations from the permanent network of seismometers operated by the Malaysian Meteorological Service (network MY, Figure 1) and two stations in Kalimantan (network 9G, Figure 1; Greenfield et al., 2018). In total, there are 72 instruments, which recorded high-quality waveforms from 2018 to 2020, with an average station spacing of ∼35 km. All the instruments are broadband seismometers but vary in their corner period from 120 to 30 s.
For subsequent surface wave modeling, we exploit a good azimuth coverage of earthquakes that have an epicentral distance between 30° and 150°, a depth <200 km and a magnitude >5.7 (Forsyth & Li, 2005;Harmon & Rychert, 2016; Figures S11 and S12 in Supporting Information S1). The instrument response is removed before using frequency-time analysis with phase-matched filtering (Goforth & Herrin, 1979;Landisman & Dziewonski, 1969) to isolate the Rayleigh Wave fundamental mode on the vertical component. At each period of interest, we measure the phase and amplitude of the Fourier transformed wave. For each measurement, we require a signal-to-noise ratio of at least 10 and observations to be made by at least eight stations (Forsyth & Li, 2005;Harmon & Rychert, 2016).

Two Plane Wave Tomography
We use two-plane-wave tomography (Forsyth & Li, 2005;Yang & Forsyth, 2006) to model the phase velocities across Sabah. The two-plane-wave method approximates the incoming wavefield as the superposition of two incoming plane waves with the effect of scattering accommodated using analytical sensitivity kernels (Zhou et al., 2004). We parametrize the velocity structure onto a grid of nodes with a grid-spacing of 0.5° across Sabah ( Figure S1 in Supporting Information S1). An outer set of nodes assigned a larger a priori uncertainty and with a spacing of 1° is included to allow for unmodeled structure outside of Sabah. We describe the continuous velocity model using a Gaussian weighting function with a characteristic length scale of 65 km. This length scale parameter imposes implicit regularization in the inversion and sets the minimum recoverable size of velocity anomaly in the final phase velocity models.

Inversion
We first proceed by inverting for a single 1D phase velocity curve for Sabah ( Figure 2a) using a standard, iterative, and linearized technique (Tarantola & Valette, 1982) alternating with a simulated annealing adjustment of the two-plane-wave parameters (Forsyth & Li, 2005;Yang & Forsyth, 2006). This step allows us to remove observations that are clearly degrading the result. As a starting model, we use a dispersion curve from the Global Dispersion Model (Ekström, 2011) for a point in the center of Sabah (black line, Figure 2a). After each inversion we remove observations with high amplitude (>0.3) or phase (>4 s) residuals and update the starting model.
Following the 1D inversion, we invert for the 2D phase velocity structure using the same inversion technique, the groomed data set and the 1D average phase velocity as the starting model ( Figure 2a). The resulting phase velocity maps ( Figure 2) show significant deviations from the 1D average (up to 0.5 km/s) and fit the data ∼50% better than the 1D model indicating such deviations are required by the data. We estimate the posterior covariance of our model and isolate regions with a standard deviation of 60 m/s or less as being well constrained by the data (Harmon & Rychert, 2016). Recovery is good between 25 and 80 s period; at longer periods, the well constrained region diminishes in size until no part of the 130, 150, and 200 s models have standard deviations less than 60 m/s.

Synthetic Tests
We perform a suite of synthetic recovery tests to investigate the ability of our data set to recover velocity heterogeneities (Figures S2-S7 in Supporting Information S1). These tests include checkerboard tests, which, despite their Mt. Kinabalu C r o c k e r R a n g e The recovery of features observed in our final phase velocity model ( Figure 2) is investigated using synthetic recovery tests. Low-velocity regions, such as those seen at short periods are well recovered, with only a slight decrease in amplitude and low-amplitude artefacts over the model area. A low-velocity anomaly to the southeast of Sabah, beneath the Semporna Peninsula (Figures 1 and 2) is recovered to periods of 55 s, but loses its definition at longer periods, while still being recognizable. Since, the observed feature is also present up to this period, we cannot say whether such a feature extends to longer periods, due to a lack of sensitivity.

Anisotropy
Azimuthal anisotropy is modeled by allowing the phase velocity to vary with back azimuth using the standard assumption of π-periodic (or 2θ) anisotropy (Forsyth & Li, 2005). We divide our modeled region into three regions: outside of Sabah, SE Sabah, and NW Sabah based on a NW-SE change in the surface geology and to avoid the structure beneath Sabah being affected by artefacts produced by the anisotropic structure in the surrounding ocean basins.
At short periods (<37 s), the orientation of the fast axis of anisotropy (fast direction) is NW-SE (Figures 2c and 2d). In the mantle (Figures 2e and 2f), our model shows a difference in the fast direction between SE and NW Sabah. Under the Crocker Range, in NW Sabah, the fast directions are NE-SW, parallel to the dominant structural trend at the surface. In SE Sabah the fast direction varies slightly but is predominantly in a N-S orientation.

Shear Wave Inversion
To calculate the shear wave velocity structure, we first extract pseudo-dispersion curves from the phase velocity maps every 0.2°. Each pseudo-dispersion curve is inverted using a transdimensional, reversible-jump, Markov-chain Monte-Carlo algorithm implemented by the Bayhunter package (Dreiling & Tilmann, 2019;Green et al., 2020). We run 80 different chains, to explore the parameter space, each with a total of 500,000 model realizations. The first 250,000 models are discarded as "burn-in" models. We discard chains that have a median misfit greater than 5% of the average misfit and take a random sample of 100,000 models to assemble the final ensemble. The final model is calculated as the mean of the ensemble at 10 km intervals ( Figure 3f). Our shear wave model shares a similar horizontal resolution and sensitivity to that of the phase velocities. Due to the variation in the size and shape of the area of good recovery at each period, calculating a region which is well-recovered at all depths is not trivial. We conservatively show a region which is covered by at least 10 phase velocity observations. This seems to avoid complicated structure near the boundary of the model caused by insufficient data to constrain the Vsv structure. In depth, the calculated error suggests that we can successfully constrain structure between 20 and 300 km ( Figure S10 in Supporting Information S1). However, because of the shape of the sensitivity kernels, the resolution of our model decreases with depth such that at depths of 150 km we cannot image anomalies <50 km thick.
At shallow depths (20 km, Figure 3a), the Vsv structure is dominated by a large low-velocity region extending through the center of Sabah in a NE-SW orientation. It is flanked on either side by high velocities beneath the Crocker Range and the Semporna Peninsula. At approximately the depth of the Moho (40 km, Figure 3c), the Vsv structure features an anomaly pattern that is opposite in polarity to that observed at 20 km depth, with high velocities across the middle of Sabah and low velocities beneath the Crocker Range and across a NE-SW oriented band in eastern Sabah. Between 60 and 80 km depth, low velocities are found in the SE corner of Sabah. At greater depths, the velocity structure is more uniform to ∼150 km and due to the reduction in vertical resolution we elect not to interpret the 2D structure below these depths.

Anisotropy
The observed anisotropy of Rayleigh Waves in NW and SE Sabah at periods from 25 to 130 s are inverted for depth-dependent anisotropy parameters G c,s , B c,s , and H c,s in each region to a depth of 300 km using the method outlined by Montagner and Nataf (1986). Using this framework, the Rayleigh-2θ cosine and sine components can be written as follows: where, c is the phase velocity, U the group velocity and K L , K A , and K F are the eigenfrequency Fréchet derivatives of the corresponding Love parameters, which depend on frequency, ω, and depth, z.  price, 1998;Russell et al., 2019). The inversion is then performed using a damped-weighted least-square approach with first and second degree smoothing.
Results (Figure 3h) for both regions show a reduction in the strength of anisotropy through the crust, reaching a minimum at Moho depths (40-50 km). Another peak in anisotropy is observed at ∼100 km depth. If we assume the depth to the lithosphere-asthenosphere boundary (LAB) is ∼100 km (Section 4.3), the mean strength of anisotropy in the lithosphere is 5%-8%, at the upper end of estimates from teleseismic shear-wave-splitting observations (Pilia et al., 2021). In both regions the direction of anisotropy is −050° at shallow depths. However, at greater depths (>50 km), the two models diverge. The fast direction is modeled as 000° in SE Sabah, while in NW Sabah, the fast direction is oriented at 045°, subparallel to the Crocker Range.
The range-parallel fast directions we observe at upper mantle depths are similar to that observed beneath young orogens (e.g., Hirn et al., 1995;Meissner et al., 2002;Petrescu et al., 2020). However, the slow present-day deformation rate suggests this would be a frozen fabric from the collision between Sabah and Dangerous Grounds approximately 21 Ma.
In the crust, the range of possible mineral assemblages is much greater than in the mantle. Thus, relating anisotropic structure to deformation is more difficult (Lloyd et al., 2009). Given the prominent, along-strike metamorphic fabric (Tongkul, 1991) and a lack of any candidate collision (Hall, 2012), it seems unlikely that the range-perpendicular fast directions (NW-SE) we observe at shallow depths could be the result of a "frozen-in" fabric. Instead, we interpret the shallow anisotropy due to crustal extension, similar to that observed beneath the Aegean (Endrun et al., 2011).

Crustal Thickness
We interpret the depth of a sharp velocity gradient between 20 and 50 km as the base of the crust (Figure 3f) and use the 4.1 km/s contour as a proxy for Moho depth throughout Sabah ( Figure S8 in Supporting Information S1). Crustal thicknesses using this contour vary between 20 and 25 km in the north of Sabah to >50 km beneath the Crocker Range and a segment of eastern Sabah (Figure 3b). It is important to note here that because we attribute the low velocities found in the SE of Sabah (Figure 3e) to slow mantle velocities, we cannot estimate the thickness of the crust in this region.
The thickest crust in our model is greater than the maximum crustal thickness of the CRUST1 model in our study region (∼35 km; Laske et al., 2013), although the average thickness is similar. To the northwest, our crustal thickness is similar to that obtained from seismic refraction lines in this region (Franke et al., 2008), suggesting that our choice of 4.1 km/s is appropriate. Between the two regions of thick crust, the crust thins to 20-25 km in a region that extends into the Sulu Sea (SS). Assuming the 40-50 km thick crust beneath the Crocker Range and Eastern Sabah represents the pre-stretching crustal thickness, we calculate a β-factor (degree of extension) between 1.3 and 2.

Temperature and Thickness of the Lithospheric
To estimate the lithospheric thickness we convert our Vsv model to temperature using the following equation (Priestley & McKenzie, 2006;Roberts et al., 2018): where, z is the depth (in km), T is the temperature (in K) and b v , c, Α, E, m, and V a are empirically fitted constants with the values 3.84 × 10 −4 km −1 , 4.72 km s −1 , −1.8 × 10 16 m s −1 , 409 KJ mol −1 , −2.8 × 10 −1 m s −1 , and 10 × 10 −6 m 3 mol −1 , respectively. The pressure is assumed to be lithostatic and is calculated as P = ρgz, where, ρ is the density (3,330 kg m −3 ) and g is the acceleration due to gravity (9.81 ms −1 ). Using Equation 2 we convert our Vsv model into temperature and estimate the depth of the LAB by assuming the ambient temperature of the asthenosphere is 1333°C ( Figure S14 in Supporting Information S1; i.e., a potential temperature of 1300°C, an adiabatic temperature gradient of ∼0.3°C/km and a depth of ∼100 km; Parsons & Sclater, 1977); thus the LAB depth is the depth to the 1333°C contour (Figure 3f).
The results show that the depth to the LAB is approximately 100 km, much thicker than the ∼50 km estimated from global tomography models (Priestley et al., 2018;Roberts et al., 2018;Schaeffer & Lebedev, 2013). In addition, there are significant variations in the depth to the LAB across Sabah. In the SE corner of Sabah, beneath the Semporna Peninsular and Sempatik Island, the depth to the LAB is ≤50 km. In the north, a shallow LAB region extends into Sabah in the NE-SW direction from the SS. To the east and west of this region, there are regions of increased LAB depth.

Extension
The Vsv model presented here allows us to investigate the crustal and upper mantle structure beneath Sabah. The region has been the site of both southeastward subduction of the proto-South China Sea and northwestwards subduction of the CS (Hall & Blundell, 1996). Despite this post-subduction setting, in our model we see no evidence for any high-velocity dipping slabs beneath Sabah. However, this is to be expected since: (a) surface waves are relatively insensitive to these dipping structures (although, notably they can image thick subducting plates, e.g., Wang & Tape, 2014); and (b) it is likely that the slabs have broken off and are freely descending through the mantle (Hall, 2013;Pilia et al., 2021). While little is known about the current depth of the slabs, it is likely that they are deeper than the maximum depth imaged in this study (200-300 km).
Our model reveals that the crustal and shallow mantle structure follows that expected from extension acting upon the region, agreeing with estimates derived from receiver function analysis (Pilia et al., 2021). Low velocities at shallow depths (Figure 3a) are likely to be due to Neogene sedimentary rocks in these regions; while our data are insensitive to the thickness of such sediments, estimates from sedimentary logging suggest that in places, it could be up to 7.5 km thick (Balaguru & Nichols, 2004;Tongkul & Chang, 2003). The spatial pattern of the low-velocities suggest that these sediments extend across the central region of Sabah in a broad NE-SW oriented swath. The lowest velocities are found on the northern coast. Further inland, the relative reduction of the velocities decreases, likely due to a change in the depth extent of the low-velocity sediments, with the depth sensitivity of our data modulating the observed velocities.
Further evidence for crustal extension can be found in the crustal thickness estimates, depth slice at 40 km (Figures 3b and 3c) and seismic anisotropy. High-velocity regions indicate the locations of thinned crust relative to the thick crust beneath the Crocker Range and Eastern Sabah. The region of thinnest crust is in the north of Sabah and extends into the SS, suggesting NW-SE orientated extension, matching the observed fast directions in the crust. Extension in the SS dates to the early Miocene and proceeded to seafloor spreading Silver & Rangin, 1991), before the construction of a subaerial island arc on top of the oceanic basement (Hall, 2013). The SS reaches its maximum depth closer to the Philippines, but seafloor >3 km depth extends to at least the continental shelf off the north-east coast of Sabah. Our results support a model in which extension continued into Sabah, causing thinning of the crust and the formation of an asymmetric sedimentary basin (Hall, 2013;Lai et al., 2021;Tsikouras et al., 2021). Other back-arc basins show similar features (e.g., Lau Basin, Parson & Wright, 1996).
Our LAB depth estimates also reveal thinning of the lithosphere in a similar location, although the spatial pattern is not the same as at crustal depths. We measure a region of thinned lithosphere (relative to the average, 100 km) extending into the SS approximately 50 km northwest of the thinnest crust. Accompanying the thinned lithosphere is a region of thicker than average lithosphere to the south-east, while to the north-east an extended region of slightly thicker lithosphere underlies the Crocker Range.

Removal of the Lithosphere
Our LAB depth map ( Figure 3d) has two regions of extremely thin lithosphere, to the south, beneath the Semporna Peninsula and to the west. To the south, our results show that the lithospheric thickness is less than 50 km, and shear velocities are ∼3.9 km/s, much lower than expected for ambient mantle ( Figure S13 in Supporting Information S1). Removal of the lithosphere and the emplacement of a hot layer beneath Borneo was proposed by Roberts et al. (2018) as an explanation for the Neogene uplift of Borneo. This model was partly guided by the shear wave velocities from global tomography models (Priestley et al., 2018;Schaeffer & Lebedev, 2013), which show low-velocities beneath most of Borneo ( Figure 3f). However, global velocity models have limited resolution beneath Borneo, because of a restricted number of seismic stations. Our Vsv model (Figure 3f) indicates that beneath Sabah, these global models are underestimating the velocity, which is instead closer to ak135/ typical continental lithosphere (Kennett et al., 1995). The lack of a continuous low-velocity layer beneath Sabah precludes widespread Neogene uplift being driven by the whole-scale removal of the lithosphere. Instead, our results suggest that any lithospheric removal occurs in localized patches possibly caused by convective instabilities and "lithospheric drips" (Houseman & McKenzie, 1982). We are aware, however, that there is no evidence in our model for a high-velocity region underlying the low-velocity layer beneath Tawau which would conclusively argue for a drip.
Removal of the lithosphere beneath Tawau explains the presence of recent (Lim, 1985) volcanic units distributed over the Semporna Peninsula, as suggested by Pilia et al. (2021). The lavas have affinity with ocean island basalts (OIB) rather than island arcs (Macpherson et al., 2010), and have modeled temperatures and pressures of melt equilibration at 1362 ± 39°C and 1.66 ± 0.24 GPa (55 ± 8 km; Roberts et al., 2018). These results are in line with the temperatures estimated from Equation 2 (1350-1370°C at 50-60 km depth) and indicate the replacement of the lithosphere with hotter than average asthenosphere and the presence of melt in the upper mantle causing the observed low Vsv. Importantly, our velocity model suggests this has happened in localized patches 50-100 km across, rather than across the whole of Borneo. We speculate that this process could have occurred in other areas of Borneo, since other recent volcanic eruptions also share OIB affinity (Cullen et al., 2013).

Conclusions
Collectively, our Vsv model, interpreted LAB depth and crustal thicknesses argue that uplift across northern Borneo cannot be caused by collision or wholescale replacement of the lithosphere by hot asthenosphere. Instead, thinning of the lithosphere indicates a primary role for NW-SE directed extension driven by back-arc spreading behind the NW subducting CS plate (Hall, 2013;Tsikouras et al., 2021). Where the crust is thin (∼25 km), we image a sedimentary basin formed in the early Miocene between the Crocker Range and East Sabah. The β-factor is estimated to be between 1.3 and 2, well-below that expected for crustal breakup and onset of seafloor spreading. Beneath the Crocker Range, crustal thicknesses reach 55 km, consistent with continental collision.
Removal of the lithosphere has occurred in discrete patches across Sabah, such as beneath the Semporna Peninsula. Our calculated temperatures reveal that, in this location, hot asthenosphere has replaced the lithosphere, which may have dripped off because of a convective instability. Across the rest of Sabah, the lithosphere is approximately 100 km thick except where extension has caused thinning of the lithosphere. Our results show the variety of crustal and lithospheric processes that occur when subduction ends and how careful interpretation of the available data is needed to unpick the observed data.

Data Availability Statement
The supplementary information included in this submission includes the observations made using the data set, station locations and further details on the synthetic recovery tests and model quality. These data, and the final model are assigned the doi: https://doi.org/10.5281/zenodo.5465032 and are available on Zenodo. Waveform data from the nBOSS array (https://doi.org/10.7914/SN/YC_2018) is under embargo until February 2023 where it will be available for download from the IRIS. Access to data from the Malaysian national network (https://www.fdsn. org/networks/detail/MY/) is restricted. Two stations from a temporary network deployed in Kalimantan (https:// doi.org/10.7914/SN/9G_2018) were also used in the study, requests for this data should be sent to TG or NR.