The Signature of Lithospheric Anisotropy at Post-Subduction Continental Margins: New Insight From XKS Splitting Analysis in Northern Borneo

The relative paucity of recent post- subduction environments globally has meant that, so far, little is known about tectonic processes that occur during and after subduction termination, as previously convergent tectonic plates adjust to the new stress regime. The region of Southeast Asia that now encompasses northern Borneo has been host to two sequential episodes of subduction—both now terminated—since the mid-Paleogene. It is expected that these processes will have left signatures in the fabric of the upper mantle, which are manifest in the form of seismic anisotropy. We investigate the evidence for, and alignment of, anisotropic fabrics by measuring the splitting of a family of teleseismic shear phases. These observations provide a measure of the orientation of the effective anisotropic elastic tensor, in the form of the orientation of the fast shear-wave polarization, ϕ

. Map of the configuration of the northern Borneo Orogeny Seismic Survey (nBOSS) network across northern Borneo. Blue and red triangles indicate nBOSS 30 and 60 s instruments, respectively. Purple squares indicate the permanent broadband stations operated by MetMalaysia. The white dashed line delineates the extent of the Western Cordillera, which encompasses the Crocker and Trusmadi Ranges. A selection of relevant units of the surface geology (derived from Hall (2013)) are delineated by shaded polygons: Ranau peridotites (blue), Kinabalu granite (red), Telupid ophiolite (purple), and crystalline basement (orange). A map of the surface geology and a fully labeled network map are available in Figures S1-S2 in Supporting Information S1, respectively.

10.1029/2022GC010564
3 of 17 emerged above sea level by around 5 Ma (Hall, 2013). The exact cause of this uplift has not been clearly identified, but it is likely to be different between western and eastern Sabah. Possible explanations include Celebes Sea slab detachment (Hall, 2013) or the development of a gravitational instability with subsequent detachment of a lithospheric drip  for eastern Sabah, and proto-South China Sea slab detachment or lithospheric delamination from the thickened region beneath the Crocker range for western Sabah (Hall, 2013). Additionally, a change to ocean-island volcanism has been recorded in Semporna Peninsula at around 5-2 Ma (Macpherson et al., 2010). The tectonic evolution of northern Borneo from the late Miocene-and exactly how each tectonic event since the Neogene is related-remains a puzzle, with the published models discussed above focusing on the crust due to limited constraints on mantle structure and dynamics. Consequently, there exist a number of possible scenarios for variations in the thickness of the lithosphere across northern Borneo and the dynamic state of the asthenosphere below. The region has seen little prior coverage of seismic instrumentation, resulting in a lack of seismic constraints on the structure of the crust and mantle.

Seismic Anisotropy
Observations of seismic anisotropy-the directional dependence of seismic wavespeeds-have long been linked to deformational processes within the Earth (Hess, 1964;Silver & Chan, 1988;Vinnik et al., 1984). Provided there exists a relationship between this deformation and the orientation of the induced anisotropic fabric, such observations can be used to make inferences on the dynamic state of the mantle (Vinnik et al., 1989), as well as instances of large-scale lithospheric deformation (Nicolas, 1993;Silver & Chan, 1991). Under finite strain, intrinsically anisotropic minerals, such as olivine (the primary constituent of the upper mantle), form a preferential alignment with respect to the flow geometry. This allows the intrinsic, crystal-level elastic anisotropy to manifest on a macroscopic scale, a phenomenon known as lattice preferred orientation (LPO) anisotropy (Nicolas & Christensen, 1987;Zhang & Karato, 1995). Under typical mantle conditions, olivine has typically been thought to form A-type LPO, wherein the a-axis of the olivine crystals are aligned parallel to the mantle flow direction. It has been shown, however, that when deformation proceeds by dislocation creep, the resultant LPO is also a function of the physical and chemical conditions (e.g., fluid content, pressure, and temperature) of deformation, which can complicate the connection between observations of seismic anisotropy and the inferred state of the mantle (e.g., Jung et al., 2006;Katayama et al., 2004).
Shear-wave splitting is a key and near-unambiguous indication of the presence of seismic anisotropy (Savage, 1999;Silver, 1996). When a linearly polarized shear wave impinges on an anisotropic medium, it is partitioned into two quasi-S waves, which propagate at different velocities. A time lag, δt, develops between these two waves as they travel through the anisotropic medium, with the final integrated value proportional to both the path length and strength of anisotropy. The polarisation of these two waves, commonly called "fast" (denoted ϕ hereafter) and "low", are controlled by the symmetry and orientation of the anisotropic elastic tensor. The orientation of the fast axis of anisotropy, ϕ, can be related to both more recent asthenospheric flow or historic lithospheric deformation (e.g., Löberich & Bokelmann, 2020;Silver, 1996), with the dominant mechanism being the LPO of intrinsically anisotropic minerals such as olivine (e.g., Karato, 1995;Zhang & Karato, 1995). Studies of azimuthal anisotropy in the upper mantle commonly make use of the family of core-refracted phases that involve a P-to-S conversion at the receiver-side core-mantle boundary (CMB), hereafter collectively referred to as XKS phases (Savage, 1999;Silver, 1996;Silver & Chan, 1988, 1991Vinnik et al., 1984). These converted phases are polarised in the radial plane and retain no information on source-side anisotropy-thus, energy observed on the transverse component is diagnostic of anisotropy or lateral heterogeneity beneath the receiver. Consequently, shear-wave splitting measurements provide excellent lateral resolution of seismic anisotropy, but lack vertical resolution.
Null measurements are measurements that suggest that there has been no splitting of the wave between the CMB and the receiver. It is possible that these observations indicate that there is no radial anisotropy along the raypath, either due to the Earth being isotropic or because the axis of anisotropy is oriented vertically, which can be the case if there is vertical flow beneath the station (e.g., Merry et al., 2021;West et al., 2009). In this instance, one would expect to observe nulls at any azimuth. However, null observations may also arise if the initial polarisation of the core-refracted phase is aligned with either the fast or slow axes of anisotropy. The latter case is hereafter referred to as a "geometric" null and is commonly indicated by null observations being limited to two azimuths, 90° apart. Geometric nulls can still be useful for constraining the orientation of anisotropy, though they (inherently) lack any information on the strength of anisotropy.
In general, shear-wave splitting measurements provide excellent lateral resolution of seismic anisotropy, but lack vertical resolution (e.g., Savage, 1999;Silver & Chan, 1991). The geometry of mantle flow in active subduction zones, for instance, has primarily been constrained by measurements of shear-wave splitting in subvertically propagating core-refracted phases, due to their superior lateral resolution compared to surface wave inversions. There may be contributions from different parts of the subduction system, including the overriding lithosphere, the mantle wedge, the slab itself, and the sub-slab mantle. Local S phases originating within the slab can provide useful constraints on the depth distribution of anisotropy (Abt et al., 2009;Bowman & Ando, 1987;Eakin et al., 2015;Fischer & Wiens, 1996;Long & van der Hilst, 2005). In the case of northern Borneo, however, any remnant lithospheric material in the underlying mantle appears to be entirely aseismic, which we speculate to be a result of subduction termination and slab breakoff. In general, the anisotropy beneath the mantle wedge in active subduction zones appears to be oriented parallel to the trench (e.g., Russo & Silver, 1994). In the mantle wedge it is often a more complex picture, with fast orientations tending to be trench-parallel close to the trench and rotating to trench-perpendicular in the back-arc of many systems, such as the Tonga subduction zone (Fischer et al., 2000;Karato, 1995).
In some regions undergoing lithospheric shortening, such as northern Tibet, the fast axis of seismic anisotropy tends to be oriented parallel the orogenic belts (e.g., Kaviani et al., 2021;McNamara et al., 1994;Nicolas, 1993;Silver & Chan, 1988, 1991, likely as a result of the formation of fabric during the collision that is subsequently frozen into the lithosphere. Observations of shear-wave splitting can therefore not only image deformation related to more recent mantle flow geometries, but also constrain the tectonic history of continental regions (Gilligan et al., 2016;Liddell et al., 2017). In the Alps, slab-rollback related toroidal flow in the asthenosphere has also been shown play a role in generating mountain-belt parallel seismic anisotropy (Barruol et al., 2011;Löberich & Bokelmann, 2022;Salimbeni et al., 2018). It is therefore important to consider the spatial and temporal scales of a subduction zone when interpreting observations of seismic anisotropy. But what happens after subduction stops? How does order established in the upper mantle evolve with the changing stress conditions as the once converging plates cease their relative motion and subduction terminates? Furthermore, what does this mean for seismic anisotropy in systems transitioning from the subduction of oceanic lithosphere to continent-continent collisions and then orogen collapse in a post-subduction environment? Answering these questions will provide fresh insight into the termination and post-subduction phases of the subduction cycle.
Here, we present the first broad-scale study of teleseismic shear-wave splitting using a network of instruments across northern Borneo, with the goal being to understand the present-day dynamics and long-term deformation of the region. The results are interpreted in the context of models that have been proposed for the tectonic evolution of northern Borneo and provide new constraints on the dynamic state of the mantle, particularly the lithospheric mantle, in a post-subduction setting.

Data and Methods
The seismic waveform data used in this study were recorded by two networks of seismic instruments ( Figure 1): the temporary nBOSS (northern Borneo Orogeny Seismic Survey) network (FDSN network code YC, https:// doi.org/10.7914/SN/YC_2018; see also Pilia et al. (2019)); and the permanent monitoring network operated by the Malaysian meteorological office, MetMalaysia (FDSN network code MY). The nBOSS network, which operated between March 2018 and January 2020, consisted of 46 seismic stations with a mean station separation of 38 km. Two instruments types were used: the Güralp 6TD, a three-component broadband instrument with a flat response between 30 s and 100 Hz; and the Güralp 3ESPCD, a three-component broadband instrument with a flat response between 60 s and 100 Hz. The MetMalaysia permanent network consists of Streckeisen STS-2/2.5 seismometers, which are three-component broadband instruments with a flat response between 120 s and 50 Hz, with their deployment focusing on the seismically active regions around Mount Kinabalu in the north-west and Darvel Bay in the south-east ( Figure 1). Restricted data recorded by the MetMalaysia network between January 2018 and January 2020 were made available to the nBOSS working group.
A catalog of events fulfilling the criteria M b ≥ 5.8 and epicentral distance ≥85° was produced, containing a total of 129 earthquakes. A list of the events used in this study and a map showing their location are available in the Supporting Information S1. Model phase arrival times for the XKS phases were calculated using a traveltime lookup table for each event/station pair. Due to the short deployment period, there are notable gaps in the back-azimuthal data coverage (see Figure S3 in Supporting Information S1), which has implications for how the shear-wave splitting measurements can be used for interpreting anisotropic structure. Earthquakes that contributed an observation of shear-wave splitting used in this study are shown as gray circles in the inset of Figure 4.
The MetMalaysia site KKM (see Figure 1c) has been in continuous operation since 2005, with continuous waveform data archived at the Incorporated Research Institutions for Seismology (IRIS) Data Management Center (DMC). A separate catalog of events, using the same criteria, was created in order to produce a long-term benchmark to help validate the observations at the temporary/more recent networks. The distribution of the 1,030 events that are suitable for analysis are shown in the Supporting Information S1 ( Figure S2 in Supporting Information S1) and those which contributed an observation of shear-wave splitting used in this study are shown as black squares in the inset of Figure 4. Consequently, the azimuthal coverage at this site greatly exceeds that of any other station in the network, enabling us to investigate the possibility of more complex anisotropic structures beneath northern Borneo, such as dipping and/or multiple layers, features that might be expected given what is known of the region's tectonic history.
The shear-wave splitting analysis was performed using the SplitRacer code package (Reiss & Rümpker, 2017). All waveform data were bandpass filtered between 3 and 25 s before a signal-to-noise ratio was calculated in a window around the predicted arrival. Only those phase arrivals with a signal-to-noise ratio exceeding 1.5 were retained for splitting analysis. The waveform data were then visually inspected and the automatically assigned analysis windows were assessed. Where necessary, these time windows were adjusted to best capture the phase arrival and exclude any non-XKS arrivals. Poor quality, noisy waveforms were also removed at this stage. The initial splitting analysis was performed using the single channel transverse energy minimization method of Silver andChan (1988, 1991), assuming a horizontal anisotropic layer. This technique consists of a grid search over the fast direction (ϕ) and delay time (δt), for the combination that best removes splitting. The uncertainties in the measurements were assessed using the statistical measure laid out in Silver and Chan (1991), with the corrections identified in Walsh et al. (2013). The resulting splitting measurements were then visually inspected and classified as "good" (clear and well-constrained splitting), "fair" (clear evidence of splitting, but less well-constrained), "null" (clear absence of splitting), and "poor" (indeterminable result). Measurements classified as "poor" were discarded in further analysis. A measurement was classified as "good" if it exhibited the following three characteristics: (a) distinct energy on the transverse component prior to correction; (b) elliptical particle motion for the horizontal components before correction, which became rectilinear after correction; and (c) the 95% confidence contour was narrower than 60° along the ϕ axis and 0.25 s along the δt axis, with "fair" observations exhibiting at least criteria (a) and (b). Additional estimates of ϕ and δt were calculated using the joint analysis method of Wolfe and Silver (1998) and the splitting intensity technique of Chevrot (2000), which incorporate the full suite of information available across multiple observations to maximize the signal-to-noise ratio of the (ϕ, δt) search grids and thus minimize the uncertainties. For the joint analysis, the (ϕ, δt) search grids for the good and null observations at each site were stacked and the global minimum extracted. This technique assumes a single layer of anisotropy with a horizontal axis of symmetry, making it prone to inaccuracy in the case of more complex anisotropy, as well as being biased toward event clusters that contribute a large number of good observations. To balance this, the multi-channel splitting intensity method (Chevrot, 2000) is also used as a means of cross-validating the measured splitting parameters. This technique solves for ϕ and δt at a single station using multiple records from different azimuths simultaneously. For values of δt ≪ the dominant period of the signal, the equations for the radial and transverse components of a split shear wave arrival (set out in Silver and Chan (1991)) reduce to show that the transverse component can be approximated as the time derivative of the radial component multiplied by a sinusoidal factor whose amplitude and phase depend on δt and ϕ, respectively. With a good azimuthal distribution of sources, it is thus possible to extract the shear-wave splitting parameters by fitting a sinusoid to the splitting vector. It is independent of-and thus, a good complement to-the single-channel methods for measurement of shear-wave splitting parameters.

Reference Shear-Wave Splitting Measurements at KKM: 2006-2020
Variations in shear-wave splitting parameters as a function of back azimuth are often considered to be a key indicator of complex, multi-layered and/or dipping anisotropic fabrics. Identification of such patterns, however, is often hampered by the temporary nature of many passive seismic experiments and the azimuthally biased distribution of earthquakes in general. KKM-a permanent station operated by MetMalaysia since 2005, situated on the north-west coast of Sabah-offers an opportunity to explore potential complex anisotropy. Figure 3a shows the misfit grid resulting from the joint splitting analysis performed at KKM using all single-channel measurements that were previously graded as good or null, which returned a best-fitting (ϕ, δt) pair of (N068°E ± 4°, 0.92 ± 0.10 s). There appears to be little to no variation in ϕ or δt as a function of back azimuth (see Figure 3b), which may indicate there is a single dominant horizontal layer of anisotropic material responsible for the shear-wave splitting observed at KKM (see Figure S6 in Supporting Information S1). In addition to this, the splitting intensity technique was applied to independently measure the splitting parameters, the results of which are shown in Figure 3c. The method returned a strong fit to the data and a best-fitting (ϕ, δt) pair of (N066°E ± 3°, 0.71 ± 0.08 s), which closely agrees with the result of the joint splitting analysis. Together with the limited evidence of back-azimuthal variation, we hereafter treat the anisotropy beneath KKM as simple and perhaps representative of the shear-wave splitting observed at most of the stations along the north-west coast of northern Borneo.

Shear-Wave Splitting Observations: 2018-2020
A total of 1,437 phase arrivals were analyzed, resulting in 687 splitting measurements ranked as "good" (151), "fair" (328), or "null" (208). Examples of good and null observations are shown in Figures S4 and S5 in Supporting Information S1, respectively. From these 687 observations, 271 are from SKS phase arrivals, 353 from SKKS, 37 from PKS, 8 from SKIKS, and 18 from PKIKS. The average delay time for the "good" and "fair" splitting measurements is 1.41 ± 0.76 s, which is greater than the global average of 1.0 s for continents (Silver, 1996). However, this average value may be biased by a small number of stations for which there are a large number of observations. Station misalignments were assessed by comparing the orientation of the first eigenvector of the XKS horizontal particle motion with the back azimuth to the earthquake (e.g., Eakin et al., 2018). The average station misalignment was calculated to be −2 ± 5°, which is below the uncertainties of the measured ϕ values for every station (see Table 1) and is in line with the assessment of instrument orientations taken in the field during deployment and retrieval. Stacked results were possible at 42 stations, shown in Figure 4 and listed in Table 1. Taking the arithmetic mean of the joint measurements of δt results in an average delay time of 0.90 ± 0.33 s, which is consistent with the average delay time expected for continental regions. Two trends in ϕ emerge from the stacked observations. Stations in the west and north-west of northern Borneo show an average fast orientation of N063°E ± 14°, sub-parallel to the Crocker range and the north-west Borneo trough. In contrast, stations in the east and south-east have an average fast polarisation of N112°E ± 19°, which is sub-parallel to both the Absolute Plate Motion (APM, N120°E; Argus et al., 2011) and direction of spreading in the Sulu Sea. A small number of stations, located primarily to the south-east of Mount Kinabalu (see Figure 1), exhibit ϕ values that are perturbed from both the dominant NE-SW and NW-SE trends. Elsewhere, the transition between the two trends appears to be sharp, occurring over ≲ 40 km.
There are notably fewer measurements in the east of northern Borneo, despite the uniform network coverage (Figure 1), which may reflect the complex sedimentary successions observed in this area (e.g., Tongkul, 1991Tongkul, , 1993, as well as an elevated level of noise degrading the phase arrivals. Indeed, the low-lying regions of eastern Sabah are significantly more densely populated than central and southern Sabah. A number of stations exhibit nulls (represented by crosses in Figure 4)-on inspection, it is likely that these are geometric in nature (the observations typically come from a limited back-azimuthal band), as opposed to them indicating that the Earth beneath the stations is isotropic. That aside, there is evidence from mapping of the lithosphere-asthenosphere boundary across northern Borneo that the lithosphere is thinned around Tawau in the south-east of Sabah (Greenfield et al., 2022;Pilia et al., 2021), which is coincident with a number of stations that exhibit nulls (see Figure 4) and intraplate volcanism on the Semporna Peninsula (Macpherson et al., 2010). Thinning of the lithosphere can induce vertical flow of the upper mantle, resulting in nulls and, at the very least, disrupting any fossil anisotropy that has been previously frozen in. The black contour represents the 95% confidence interval. The blue plus symbol represents the optimal (ϕ, δt) pair, which is (N068°E ± 4°, 0.92 ± 0.10 s). Panel b shows ϕ as a function of back azimuth. The dashed line shows the average ϕ determined from the joint analysis method of Wolfe and Silver (1998). Black squares and gray diamonds represent good and fair measurements, respectively, with corresponding uncertainty measurements. Panel c shows the best-fitting splitting vector (dashed line) to observations of splitting intensity. The phase and amplitude of the sinusoid give ϕ and δt, respectively. The black squares and gray diamonds represent the same measurements as in panel (b) The light gray circles represent null measurements.
There is some variation in δt across northern Borneo, most notably in the south-central region around the Maliau Basin, where there is a marked reduction compared to measurements at nearby stations. This correlates somewhat with the edge of the orogenic belt, but may also be a result of more complex anisotropy as a function of depth beneath this region; without sufficient back-azimuthal coverage, it is difficult to be conclusive. Overall, however, there appears to be little variation in ϕ as a function of back azimuth, suggesting that a single dominant horizontal layer of anisotropy is sufficient to explain the shear-wave splitting observations (Silver & Savage, 1994). This is consistent with the conclusions of Song et al. (2021), who looked at teleseismic shear-wave splitting at a handful of permanent stations across the Sundaland block.

Lithospheric Origin of Azimuthal Anisotropy
The observed patterns of shear-wave splitting are consistent with the signatures of some of the major tectonic events that have occurred in the last 20 Ma: termination of southward subduction of the PSCS and associated orogeny; and extension induced by Sulu Sea rifting and/or slab rollback due to northward subduction of the Celebes Sea. While the contribution to the observed delay time from the crust cannot be constrained, it has been shown that it typically amounts to around 0.1-0.5 s (Barruol & Mainprice, 1993;Silver, 1996), thus making a mantle contribution necessary to explain our observations. It is possible to estimate the thickness of the anisotropic layer, L, from δt, using the expression L ≈ (δt × V S )/dV S (where V S is the shear velocity and dV S is the average percent anisotropy). Using a V S of 4.48 km s −1 (ak135; Kennett et al., 1995) and a dV S of 4% (an upper limit of the degree of anisotropy in the upper 200 km; Savage, 1999), the mean delay time of 0.9 ± 0.33 s corresponds to a layer thickness of 100 ± 40 km. An estimate of the depth of the lithosphere-asthenosphere boundary beneath the nBOSS network was recently extracted from a shear-wave velocity model, calculated through an inversion for surface wave phase velocities at periods between 25 and 200 s (Greenfield et al., 2022;Pilia et al., 2021). This study used a two-plane-wave approach and converted from V S to temperature, before taking the lithosphere-asthenosphere boundary to be the 1333°C contour (see Figure 4). This model exhibits an average lithospheric thickness of ∼100 km, with notable thinning (∼50 km) around the Semporna Peninsula. Hence, the lithosphere beneath northern Borneo-with frozen-in anisotropy-is thick enough to produce the observed delay times for reasonable values (∼4%) of shear wave anisotropy.
Our results are supplemented by a small number of observations from three previous studies (Cao et al., 2021;Song et al., 2021;Xue et al., 2013), which measure ϕ at the permanent MetMalaysia stations KKM and LDM. The ϕ values measured in these studies are consistent with our value at KKM and with the stations in the vicinity of LDM (see Figure 4). Xue et al. (2013) only published a single result at KKM, with a δt significantly larger than that observed in this and other studies (Cao et al., 2021;Song et al., 2021).
The orogen-parallel trend observed across western and north-western Sabah indicates that the observed seismic anisotropy may reflect continent-continent collision between the extended Dangerous Grounds continental margin and the continental margin of northern Borneo, which occurred at roughly 21 Ma. This is consistent with the conceptual framework of transpressional vertically Note. See Figure S2 in Supporting Information S1 for station locations. Instrument types are available in the station file provided as part of the Supporting Information S1.  (Silver & Chan, 1988, 1991, in which continental plates deform coherently over their depth extent and can play a significant role in mantle anisotropy. This trend spans the entire area known as the Western Cordillera in northern Borneo (dashed white line in Figure 1), extends into the offshore fold-and-thrust belt, and also appears to continue off the northern tip of Borneo toward Palawan. Fast orientations of seismic anisotropy are commonly observed to trend parallel to orogenic belts in continental collision zones, both in the crust (e.g., Fry et al., 2010;Pilia et al., 2016) and in the lithospheric mantle (e.g., Bastow et al., 2007;Gilligan et al., 2016;Helffrich, 1995;Kaviani et al., 2021;Liddell et al., 2017;McNamara et al., 1994;Silver & Chan, 1988, 1991, with transpression being proposed as the source of strain (through which LPO anisotropy is induced) normal to the relative motion between the two plates (Nicolas, 1993). Any fabric within the upper mantle relating to subduction of the proto-South China Sea has likely been overprinted by continent-continent collision. This has implications for the lifespan of anisotropic fabrics within the ductile asthenosphere, suggesting they are much more transient than fabrics left in the lithosphere during large-scale tectonic deformation events.
The mean fast orientation in the east of Sabah is sub-parallel to the APM of the Eurasian plate, though the slow-moving nature of this plate might mean that the relative motion between the plate and the underlying asthenosphere is insufficient to organize the flow within the mantle. In this case, one could discount the dominance of simple asthenospheric flow in controlling the large-scale coherence and alignment of shear-wave splitting observations (Silver, 1996) in favor of an alternative mechanism. While the orientation inferred from our observations of shear-wave splitting is not incompatible with APM-induced LPO anisotropy, the plate velocity is only around 2 cm yr −1 (Argus et al., 2011), which is less than the empirical 4 cm yr −1 proposed by Debayle and Ricard (2013) as necessary for a plate to organize the flow in its underlying asthenosphere. Hence, it is worth also considering other potential sources for the observed trend, the most prominent among these being the extension and rifting of the nearby Sulu Sea. Until recently, it was believed that northern Borneo was, and still is, under a compressional tectonic regime. However, geochemical analysis of a suite of samples from an ophiolitic complex near Telupid, central Sabah, provided evidence for the continuation of rifting in the Sulu Sea inland, dated to around ∼9 Ma ( Figure S1 in Supporting Information S1; Tsikouras et al., 2021). This inference is supported by crustal thickness measurements made across northern Borneo using the nBOSS data set , which show a thinner crust (∼30 km), extending from the Sulu Sea toward central Sabah, coinciding with the exposed ophiolite complex (see Figures S7 and S8 in Supporting Information S1). From these observations, it is also possible to conclude that the mean fast orientation across eastern Sabah is sub-parallel to the direction of the extension that would be expected due to rollback of the subducting Celebes Sea slab and opening of the Sulu Sea. In such extensional and/or rifting environments with vertically coherent deformation, stretching lineation, and hence ϕ, is expected to be parallel to the extension direction (Silver, 1996).
The lack of evidence for back-azimuthal variations in fast polarizations at KKM provides an additional line of evidence for horizontally layered anisotropy that is vertically coherent. It remains possible that there are multiple layers with similar orientations, though without earthquakes within this layer the vertical integration of δt makes this indistinguishable from a single layer. Additionally, it only provides such a constraint for a small region around KKM, corresponding roughly with the width of the first Fresnel zone for the XKS phase arrivals. There may be multiple layers of anisotropy beyond this region that we are unable to identify due to the limited back-azimuthal coverage of our shear-wave splitting observations.
A small number of stations, located primarily to the south-east of Mount Kinabalu (see Figure 4), exhibit ϕ values that lie between the two dominant NE-SW and NW-SE trends. The surface geology in the region around Mount Kinabalu exhibits an exposed suite of peridotites, dated to ∼10 Ma-thus post-dating the formation of the Crocker Range. Geochemical analysis of zircons found within these peridotites points toward a subcontinental lithospheric mantle origin for these rocks, meaning they must have been uplifted and exposed within the last 10 Ma. A recent P-wave tomographic model for the mantle beneath Sabah has revealed a narrow, fast anomaly extending from the region of thinned lithosphere around the Semporna Peninsula (south-east Sabah), toward Mount Kinabalu, interpreted to represent a lithospheric drip that formed from the root of the Sulu Arc . Thermo-mechanical modeling of a Rayleigh-Taylor gravitational instability, seeded by a small density perturbation representing the root of the volcanic arc, has been found to provide a reasonable explanation for the observed extension, and lithospheric thinning, around Mount Kinabalu . Consequently, these observations may be an indication of a gradual overprinting of the orogen-parallel trend by localized extension related to the formation of a lithospheric instability from the root of the Sulu Arc. The relationship between events and processes in the tectonic history of the region and our observations is summarized in Figure 5.
Regional-scale tomographic models of Southeast Asia exhibit a high-velocity anomaly at between 200 and 300 km depth beneath northern Borneo (Hall & Spakman, 2015;Zenonos et al., 2019), which is attributed to cold, remnant subducted lithospheric material from either the Celebes Sea or the proto-South China Sea. Song et al. (2021) attribute the fast orientation measured at KKM to LPO anisotropy induced by the deflection of mantle flow around the fossil slab segment of the proto-South China Sea subduction (inferred from tomographic studies, e.g., Hall & Spakman, 2015) or a thick lithospheric keel. Pilia et al. (2021) have demonstrated the absence, however, of any strong anomaly in the mantle above 200 km that could potentially suggest that the imaged slab material has broken off and begun to sink into the lower mantle, the minimum depth of which is thought to be ∼250 km. Consequently, there is unlikely to be any remnant of 3-D mantle flow around the down-going slab or within the former mantle wedge, as is commonly seen in active subduction zones. For simple  (Argus et al., 2011). These results are shown on top of the lithosphere-asthenosphere boundary depth (i.e., the lithospheric thickness) as determined by Greenfield et al. (2022). Inset: Distribution of earthquake sources in a polar projection centered on northern Borneo, with squares and circles representing earthquakes in the periods 2006-2020 and 2018-2020, respectively. Figures S7-S9 in Supporting Information S1 show the same results on top of two crustal thickness models and radial anisotropy, derived from a recent FWI tomographic study of Southeast Asia (Wehner et al., 2022). asthenospheric flow, the common mechanism by which strain is localized in the asthenosphere, the memory of this flow direction is thought to be only a few million years (Silver et al., 1999) based on the amount of strain required to completely reorient olivine aggregates (Nicolas et al., 1973). This interpretation-that anisotropy is due only to simple asthenospheric flow-also neglects the possible contribution due to fossil anisotropic fabric induced during the collision between the South China Sea continental margin and northern Borneo, which has since been frozen into the lithosphere. There is a strong correlation between the lateral extent of this ∼NE-SW  trend and the region showing a thickened lithosphere in Greenfield et al. (2022), likely a result of the aforementioned continent-continent collision.
Without local, deep earthquakes beneath the study region, it is difficult to build any strong constraints on the depth distribution of seismic anisotropy. The short length scales (≲40 km) over which the fast orientation of anisotropy is observed to vary, however, does indicate that at least some of the observed shear-wave splitting is likely lithospheric in origin, so as to avoid a significant overlap in the first Fresnel zones of the phase arrivals at the network, which have a width of ∼150 km at the base of the lithosphere. Future work could make use of P-to-S conversions at significant boundaries, such as the Moho, or deep regional earthquakes from the Philippines and Sunda subduction zones (though the latter would primarily be useful for improving the azimuthal coverage, rather than the depth resolution).

Radial Anisotropy
Radial anisotropy within the Earth can be constrained by the extraction of laterally averaged and from Rayleigh and Love wave phase velocities, where and are the wave speeds for vertically and horizontally polarised shear waves, respectively. The ratio of these two averaged velocities can be used as a measure of the degree of radial anisotropy as a function of depth, which can in turn be related to the tectonic processes that generate seismic anisotropy. In contrast to shear-wave splitting measurements, this technique has good vertical resolution, but at the cost of poorer lateral resolution. For A-type olivine fabric (the dominant fabric formed under standard mantle conditions i.e., poorly hydrated mantle at a range of stresses (Karato et al., 2008)), ∕ > 1 (positive radial anisotropy) indicates horizontal flow or simple shear, and ∕ < 1 (negative radial anisotropy) indicates vertical flow or pure shear shortening. Deformation by pure shear shortening may be an important factor in the generation of anisotropy in regions undergoing lithospheric shortening and thickening, proposed as an important stage in the generation of stable continental roots (Priestley et al., 2021). A key feature one would expect to see if there were significant anisotropy generated by simple asthenospheric flow is positive radial anisotropy below the lithosphere (≳100 km).
We extracted 1-D depth profiles (Figure 6) of the radial anisotropy parameter, = ( ∕ ) 2 , from the SASSY21 model, a recent 3-D full-waveform inversion tomographic model of Southeast Asia (Wehner et al., 2022). This model included waveform data down to a period of 20 s (slices through the ξ model are shown for various depths in Figure S10 in Supporting Information S1). Consequently, the velocities (and hence anisotropic structure) of the upper 50-100 km of the Earth are smoothed somewhat over this depth range, though the observed trend and sign in ξ remain valid. The averages of and were calculated at each depth slice within a small region of the model centered on northern Borneo, noting that SASSY21 also exploited nBOSS waveform data, so there is good coverage in our study area. In addition, this region was subdivided into two sub-regions, encompassing the two trends observed in the teleseismic shear-wave splitting data set, and 1-D profiles of the radial anisotropy were calculated.
The radial anisotropic structure of the upper 50-75 km appears to be positive (ξ > 1), indicating an anisotropic fabric with a primarily horizontal orientation that is consistent with LPO induced by both shortening and transpression at convergent boundaries (the orogen-parallel trend in the north-west of Sabah) and extension of continental material (Sulu Sea spreading trend in the south-east of Sabah). There is no indication of positive radial anisotropy at asthenospheric depths, which supports the conclusion that the observed seismic anisotropy is principally attributable to the lithosphere and corresponding mechanisms of deformation. There is some difference in the amplitude of negative radial anisotropy between 125 and 250 km depth, which may indicate some changes in the dynamic state of the mantle or differences in the composition/fabric between these two regions, but this is beyond the scope of this study. The thickness of this observed layer is roughly consistent with the thickness of the lithosphere observed in Greenfield et al. (2022).

Conclusions
We have investigated seismic anisotropy across northern Borneo using shear-wave splitting analysis of XKS phases recorded between March 2018 and January 2020. This has been supplemented by shear-wave splitting measurements at KKM, a permanent station operated by MetMalaysia, between January 2006 and January 2020.
Across western Sabah, most notably beneath the Crocker Range, we find orogen parallel anisotropy striking at N063°E ± 14°, likely reflecting "fossil" anisotropy in the lithosphere formed during the collision of the continental Dangerous Grounds with the continental margin of northern Sabah. Indeed, the observed delay times of 0.6-1.8 s are consistent with a 60-180 km thick layer with 4% anisotropy (or, conversely, a 100 km thick layer with 2.4%-7.2% anisotropy). Our observations do not definitively constrain whether the anisotropy beneath northern Borneo is simple or complex (e.g., due to multiple or dipping layers), though we do show that a simple model is sufficient to explain our observations. This may indicate that the sub-lithospheric mantle is either isotropic, weakly anisotropic, or possesses the same anisotropic fabric as the lithosphere, but further analysis, particularly incorporating information derived from surface waves, would be beneficial. In eastern Sabah, the orientation of anisotropy is nearly orthogonal to the trend observed in the west and is sub-parallel to both the absolute plate motion and the direction of spreading of the Sulu Sea, striking at N112°E ± 19°. This trend extends inland to central Sabah, where it terminates, correlating well with the extent of onshore spreading inferred from recent U-Pb dating. The rapid transition between these two dominant trends (over ≲40 km) suggests that the anisotropic source is shallow. A reduction in the observed delay times seen in south-central Sabah, including a number of null stations, may indicate the presence of vertical mantle flow, such as that induced by a lithospheric drip. A number of stations in the vicinity of Mount Kinabalu exhibit fast orientations that lie between the two main trends, possibly indicating the gradual overprinting of the fossil fabric in response to localized extension and thinning of the lithosphere. These results constitute strong evidence for a system of mechanisms focused predominantly within the lithosphere as the primary controls on seismic anisotropy in this post-subduction continental setting, with little influence from present-day mantle flow.