Fate of Forearc Lithosphere at Arc‐Continent Collision Zones: Evidence From Local Earthquake Tomography of the Sunda‐Banda Arc Transition, Indonesia

A new 3‐D seismic P wave velocity model of the Sunda‐Banda Arc Transition from local earthquake tomography reveals (i) northward subduction of oceanic lithosphere, associated with the convergence of Australia and Sundaland, as a high‐velocity zone extending down to ~200 km depth; (ii) two distinct low‐velocity zones, one immediately above the slab, which is likely a zone of partial melt, and one in the 0–40 km depth range, which is probably a magma chamber associated with active volcanoes above; and (iii) a northerly dipping high‐velocity zone that bisects the two low‐velocity anomalies, which we interpret as an underthrust forearc sliver of continental origin. Based on He3/He4 measurements from volcanic outflows, it is possible that the magma supply is contaminated by interaction with this forearc sliver as it ascends from the melt region below.


Introduction
Indonesia lies in one of the most active tectonic regions of Southeast Asia where convergence between the Eurasian, Australian, and Pacific plates has occurred since the late Oligocene (Hall, 2002;Hamilton, 1979). Accommodation of this convergence is not limited to the Sunda Arc, which wraps its way around the southern limits of the Indonesian archipelago but extends over 1,000 km into the back-arc region (Simons et al., 2007). A range of features can be found throughout the influence zone, including transfer of material at subduction zones, ingestion of microcontinents, lateral escape, island arcs, and abundant magmatism (Hall, 2011). In this study, we seek to understand how rapid northward motion of the Australian plate has impacted on the unique tectonic setting of eastern Indonesia, with a particular focus on the fate of continental lithosphere, which it appears must have been subject to significant deformation, subduction, underthusting, and/or obduction in order to accommodate convergence (Fichtner et al., 2010).
The Sunda-Banda Arc has a tectonic setting that transitions from subduction along the Sunda Arc to arc-continental collision along the Banda Arc (Shulgin et al., 2009). This setting represents a unique location in which to investigate how plate margins evolve at the interface between subduction and collision. The current phase of subduction beneath the Sundaland core of SE Asia began in the Eocene, coincident with the rapid northward migration of the Australian plate (Curray, 1989;Hall, 2011). This northward movement is currently associated with the subduction of oceanic lithosphere of the Indo-Australian plate beneath Java and Sumatra and with the collision of Australian continental lithosphere with the western Banda arc along the islands of Flores, Sumba, and Timor that began in the Pliocene (Hall, 1997(Hall, , 2002Keep & Haig, 2010;. We refer to this region as the Sunda-Banda Arc Transition (SBT), which is marked by a change from Australian to Eurasian tectonic units, as seen in the Savu and Sumba Islands . The existence of back-arc thrusting to the north of Sumbawa and Flores Islands is suggested by thrust and fold structures detected in marine seismic surveys (Silver et al., 1983) and by the distribution of earthquake activity (McCaffrey & Nábělek, 1984). This back-arc thrust system, which we refer to as the Flores Back-arc Thrust, seems likely to have an important bearing on the evolution of the convergent margin.
Global Positioning System measurements of surface deformation in this region demonstrate how the convergence between the Australian Plate and Sundaland transitions from 95% accommodation at the subduction megathrust off eastern Java to having most of the accommodation taking place along the Flores Back-arc Thrust (40%); this is in addition to strike-slip motion along the Semau Fault, which stretches from the western tip of Timor to the island of Alor east of Flores (Koulali et al., 2016;Nugroho et al., 2009). This variation in the accommodation of convergence in the SBT is associated with the lateral offset of oceanic subduction and continental collision . Tate et al. (2015) show that between 215 and 229 km of Australian continental crust has been consumed during arc-continent collision.
The northerly dip of the Indo-Australian slab is 10-30°from the trench to the arc, after which it dips more steeply down to the mantle transition zone (Widiyantoro et al., 2011). The average slab dip in the Sunda Arc, where oceanic lithosphere is being subducted, is~55-60°between a depth of 200 and 400 km, while the slab dip along the Banda arc, where the remaining oceanic lithosphere has effectively been consumed by ongoing subduction, is~75-80° (Royden & Husson, 2009). The initiation of this subduction is likely to have involved reorganization of the microplates in the SBT, for which there are several competing models for explaining the current configuration of slabs and troughs (Audley-Charles, 2004;Linthout et al., 1997). Finally, we note that studies of volcanic products along the Sunda-Banda Arc show that there are pronounced along-arc variations in radiogenic isotope data, with Sr, Nd, and Pb (Vroon et al., 1993); Pb and He (Elburg et al., 2004); and O (Vroon et al., 2001) isotopes all indicating that volcanoes in the SBT contain strong signatures of continental material, which contrasts with volcanoes elsewhere along the Sunda-Banda Arc.
A National Science Foundation funded project  in collaboration with the Bandung Institute of Technology and BMKG installed 30 temporary seismic stations to form the YS network (http://www.fdsn.org/networks/detail/YS_2014/) in the Nusa Tenggara to Timor region (the green inverted triangles in Figure 1). The goal of this study is to use data collected by the YS network to image the P wave velocity structure (Vp) beneath the SBT, in order to better understand the interaction between oceanic and continental plates. Our new 3-D velocity model of the crust and upper mantle structure beneath the SBT is based on arrival time data from local earthquakes recorded by the YS network. The application of SIMULPS (Evans et al., 1994;Thurber, 1983) facilitated the recovery of a detailed Vp model, which not only agrees with earlier findings  but also exhibits a number of previously unseen features which help to explain the structure and dynamics of this unique tectonic setting.

Data and Methods
The temporary YS seismic network consisted of 30 broadband seismic stations deployed in the SBT between March 2014 and March 2016 ( Figure 1). Of the 30 stations that were deployed, 8 were Nanometrics Trillium 120PA sensors with Taurus digitizers and 22 were Streckeisen STS-2 sensors with Kinemetrics Q330 dataloggers. All stations recorded continuously on three components at 50 samples per second . In total, 567 earthquakes across a depth range of 6-455 km were recorded by the network over a 24 month period of operation. The arrival times of direct P and S waves were detected using manual picking ( Figure S1 in the supporting information) with Seisgram2K (Lomax & Michelini, 2009). Initial locations were determined using the Hypoellipse code (Lahr, 1979) with the 1-D velocity model AK135 (Kennett et al., 1995).
We then applied the double-difference method (Waldhauser & Ellsworth, 2000) using the HypoDD program (Waldhauser, 2001) to determine accurate relative hypocenter locations based on the picked arrival times of P and S waves. The method achieves its accuracy by taking advantage of the approximation that if the hypocentral distance between two earthquakes is small relative to the event station distance and scale length of velocity heterogeneity, then the raypaths from the two events can be considered identical beyond the source region. Following application of HypoDD to our data set, the mean epicentral correction is 8.7 km, and the mean depth correction is 11.3 km.
We compare the HypoDD relocated events with our initial locations derived using Hypoellipse in Figure S2. The relocated hypocenters plotted in the vertical cross sections ( Figure S3) indicate that most of the earthquakes are located within the Wadati-Benioff zone, which is consistent with the subduction setting. Shallow seismicity (6-25 km depth) appears to cluster in the north of the islands that lie between Sumbawa and Flores and is probably related to the Flores Back-arc Thrust that is thought to be present in this region (Widiyantoro & Fauzi, 2005). We employ the criterion that each earthquake used for tomographic inversion is recorded by at least six YS seismic stations, which reduces the initial set of 567 detected earthquakes down to 415. Figure S4 shows the distribution of epicenters, seismometers, and the inversion grid (which is used to define variations in velocity structure) employed in this study. To recover the 3-D velocity structure and relocate the hypocenters more precisely, we applied the well-known and commonly used 3-D inversion code SIMULPS (Evans et al., 1994;Thurber, 1983), which implements a pseudo-bending method of ray tracing and a damped least squares inversion approach. The code allows for the simultaneous inversion of P and S arrival times for Vp and Vp/Vs; this has been shown to be more effective than separately

Geophysical Research Letters
inverting for Vp and Vs and then computing Vp/Vs by division, due to the much smaller number of S picks available compared to the number of P picks and the application of ad hoc regularization (Eberhart-Phillips, 1986).
Due to the underdetermined nature of the inverse problem, regularization is necessary in order to produce a stable and robust solution. The regularization approach we use is damping, which favors models that are close to the starting model. In this study, the damping factor was empirically determined by using the trade-off between data and model variance (Eberhart-Phillips, 1986). The appropriate damping value was chosen to be where the data variance is minimized with as small a model variance as possible. Based on the results of this approach and several inversion tests, the appropriate damping factor is 60 for the Vp inversion ( Figure S5). Vp anomalies are computed with reference to the 1-D velocity model AK135 (Kennett et al., 1995), which is the starting model used in the iterative nonlinear inversion process. To recover the seismic velocity structure beneath the SBT zone and its surroundings, we used a total of 2,743 P wave arrival times ( Figure S6). In this study, we do not draw inferences from the Vp/Vs results since the number of good S picks was only around 50% of the number of good P picks, which resulted in a Vp/Vs model of much lower resolution than Vp; consequently, it was unable to usefully constrain the primary structures that we interpret from the Vp model.

Resolution Tests
To investigate the robustness of the recovered anomalies, we carried out checkerboard recovery tests. Positive and negative perturbations of ±10% relative to the 1-D reference velocity model were used as inputs in the tests (with the same grid spacing as used in the inversion of observations), and the associated synthetic traveltime data set was generated by the forward solver in SIMULPS. If the inversion results from the synthetic data-which use the same damping as the model produced from the observations-are able to recover the corresponding input anomaly, even with a reduction in amplitude due to the implementation of damping, we consider that anomaly to be resolved by the data. This provides a means of deciding which regions of the model can be interpreted with confidence, while keeping in mind the limitations of checkerboard tests (Rawlinson & Spakman, 2016). The results of checkerboard tests for horizontal slices and cross sections are presented in Figures S7 and S8. It is clear from these results that structural recovery beneath much of the islands of Flores, Sumba, and Timor is sufficient to permit interpretation across the depth range of the model, while some smearing can be observed between Sumba and Timor Islands.

Vp Model
The arrival time residuals associated with the initial and final P wave models are plotted as histograms in Figure S9; the RMS residual associated with the initial model is 2.8 s, while for the final model the RMS residual is 0.8 s. We visualize the final 3-D model as a series of key depth sections in Figure 2 and vertical cross sections in Figure 3. A significant feature at 20 km depth is the low-velocity anomaly which underlies Timor and may be a signature of passive continental margin crust being underthrust as a result of continent-arc collision; the location of this anomaly correlates well with the northern limit of the Australian plate according to Harris (1991). The high velocities at 50 km depth correspond to lithospheric mantle of the upper plate, with some contribution from the lower plate, although this can be observed more clearly in the cross sections discussed below. Below 70 km depth, high velocities tend to be associated with the northerly subduction of oceanic lithosphere.
Section A-A′ of Figure 3 depicts a low-velocity anomaly below Sumba and Flores islands at depths of 20 and 20-50 km, respectively. A clear and continuous high-velocity anomaly that dips north was detected between 50 and 200 km depths and is likely related to the subducted slab beneath this region. Cross-section B-B′ depicts a low-velocity anomaly detected beneath volcanoes on the island of Flores. As in section A-A′, the high-velocity anomaly beneath the Savu Sea, which is continuous to 200 km depth, is inferred to be the subducted slab. Section C-C′ crosses the islands of Rote and Flores; this section shows a high-velocity feature in the upper crust under the Savu Sea and a significant low-velocity zone in the lower crust or uppermost mantle below the island of Flores. Below this low-velocity anomaly is a high-velocity anomaly at a depth of 70 km, which is underlain by another low-velocity anomaly that is continuous to 120 km depth, which we interpret as the signature of partial melting. The depth sections in Figure 2  observations, although they also show that the low-velocity anomaly below Flores extends eastward beneath Alor, and there is a relatively low velocity anomaly to the NE of Sumba at~50 km depth. We also detect a low-velocity anomaly below Timor at depth of~10-20 km and a high-velocity anomaly at depth of 50-120 km; Porritt et al. (2016) discovered a low shear wave velocity anomaly below 15-25 km depth in this region, but at depths greater than 35 km their model exhibits a high shear wave velocity anomaly.
Based on cross-section B-B′, C-C′, and D-D′ of Figure 3, the sandwiched high-velocity layer at~70 km depth is interpreted as likely being underthrust hyperextended forearc continental mantle lithosphere. Similar underthrusting of a sliver of the Luzon forearc beneath Taiwan has been imaged in a tomographic study of the Taiwan arc-continent collision (Cheng, 2009). Our interpretation is corroborated by (i) the thickness of the high-velocity sliver, which at 10-20 km is thinner than one might expect from underthrusting of the Australian plate and (ii) the progressive disappearance of the forearc from west to east (Harris, 2006(Harris, , 2011. The fate of this forearc, which changes in length from 300 km in the west (Line A of Figure 1) to~30 km near Dili, has long been a mystery, but our new tomographic images appear to resolve this longstanding question.
Both Fichtner et al. (2010), using earthquake waveform tomography, and Porritt et al. (2016), using ambient noise tomography, inferred the presence of subducting continental lithosphere in the SBT; in the latter case this was to a maximum depth of~40 km, which represented the resolution limit of the surface wave data they invert. McBride and Karig (1987) find a mass deficit, implied by a regional gravity anomaly, in the upper mantle or crust beneath the Savu Sea that they interpret to be caused by a combination of subducted and lighter continental crust, downward flexing of the upper crustal plate in the forearc area, and possibly anomalous low-density upper mantle below the Savu Sea. This is consistent with previous geochemical studies (Hilton & Craig, 1989;Hilton et al., 1992) of volcanoes in this region which show that they contain anomalously low He 3 /He 4 ratios, which indicate continental contamination, that is, the presence of melted continental material in igneous samples. Based on the results from our study, it may be that this contamination comes partly from the underthrust forearc (Figure 3b), which Harris (2011) associates with the Banda Terrane, originally part of a proximal forearc basin of a continental arc.
We also detect a southward dipping, high-velocity anomaly in the lower crust and uppermost mantle north of the island of Flores (Figure 2, section C-C′), which coincides with the position and orientation of the EW

Geophysical Research Letters
striking Flores Back-arc Thrust (Figure 1). This could represent incipient southward subduction of Flores lithosphere, in a polarity reversal of subduction as suggested originally by Hamilton (1979).
Based on the Vp model and the results of prior work in this area (Widiyantoro et al., 2011), we present a schematic interpretation for the SBT zone in Figure 4 in terms of a western section (based on Figure 3a) and an eastern section (based on Figure 3c). The western section (Figure 4a) is more consistent with what is traditionally observed at a subduction zone, with a high-velocity slab descending beneath the forearc, although evidence from a number of studies (Audley-Charles, 1985;Fortuin et al., 1994;Rigg & Hall, 2011, 2012) have demonstrated that Australian crust of the Scott Plateau has been thrust beneath Sumba.
By contrast, the eastern section ( Figure 4b) represents a more advanced collision zone, spanning the region from east of Sumba to Timor Island. In this region, we interpret that a forearc sliver has become entrained in the subduction process and resists significant plate bending due to increased buoyancy, resulting in shallow underthrusting beneath the Savu Sea, central and eastern Flores, and western Timor. This produces the high-velocity anomaly above the main subducting oceanic slab. The "thinness" of the underthurst mantle forearc sliver is likely due to the crust being scraped off (obduction) and the fact that the forearc comprises hyperextended lithosphere, and hence a thinner than normal high-velocity zone.
The presence of the forearc sliver between the zone of partial melting at depth and the shallower magma chamber (Figure 4b) may well contribute to the contamination of the ascending melt with continental material, which helps explain the anomalous isotope ratios observed in eruptive materials from the overlying volcanoes describe above. While it is possible that a significant contribution also comes from the basement terrane of the upper plate, it is interesting to note that the volcanism in western Flores-where there is no underthrust forearc sliver-is not anomalous despite also having a basement of continental provenance.

Concluding Remarks
Local earthquake data recorded by the temporary YS seismograph network have permitted us to image several important structural features beneath the SBT, especially at depth above 120 km. We have detected the

Geophysical Research Letters
presence of a high-velocity anomaly north of Flores Island that could represent southward underthrusting of Flores Sea lithosphere along the Flores Back-arc Thrust, consistent with the hypothesis of subduction polarity reversal initiation originally proposed by Hamilton (1979). We suggest that low-velocity anomalies in the crust directly beneath the volcanoes in the inner arc islands of Flores to Alor, and in the uppermost mantle just above the north dipping slab, could indicate the presence of partial melt and/or fluids that drive the active arc volcanism. Finally, we found that between the low-velocity zones at the top of the subducting slab and in the lower crust, there is a north dipping high-velocity anomaly, which we interpret to be evidence of an underthrust sliver of forearc lithosphere. This implies that fluids expelled by the subducting slab, and the consequent partial melting, would have to traverse the entrained fragment of continental lithosphere before they reach the lower crust, where they drive volcanism in the SBT. This could help explain measurements of isotope ratios in igneous rocks of the SBT that exhibit clear continental signatures (Elburg et al., 2004;Vroon et al., 1993;Vroon et al., 2001). The new tomographic results provide one of the clearest examples to date of what happens to forearc lithosphere at arc-continent collision zones and demonstrates its importance in influencing dynamic process at plate margins.