P and S Wave Anisotropic Tomography of the Banda Subduction Zone

In subduction zones with slab‐slab interactions, the pattern of mantle convection is very complex and still unclear. In this study, we jointly invert a large number of P and S wave arrival time data of local earthquakes for 3‐D isotropic and anisotropic velocity structures of the Banda subduction zone. Along the curved Banda arc, the subducting Indo‐Australian slab is detected clearly as a high‐velocity zone, and its azimuthal anisotropy changes along the arc strike, representing fossil anisotropy within the slab and modified anisotropy by the subduction processes. Around the northern edge of the Banda slab, a semi‐toroidal pattern of anisotropy appears in low‐velocity anomalies, representing mantle flow extruded from the Banda arc and escaped from a gap of the Banda‐Molucca slab toward the northeast. Our 3‐D anisotropic tomography uncovers the mantle convection pattern induced by the slab‐slab interactions, shedding new light on the complex dynamical processes in this curved subduction zone.

• The first P and S wave anisotropic tomography of the Banda subduction zone is obtained • Trench-parallel anisotropy along the Banda arc reflects lateral mantle flow caused by the highly curved Banda slab • Semi-toroidal anisotropy around the northern edge of the Banda arc reflects extruded mantle flow from the arc

Supporting Information:
Supporting Information may be found in the online version of this article.
The Philippines-Banda area is located at the junction of the Philippine Sea (PHS) plate, the Eurasian plate, and the Indo-Australian plate (Figure 1; Cardwell & Isacks, 1978;Hall, 2002;Hall & Spakman, 2015).The PHS plate is subducting beneath the Eurasian plate along the southern Philippine trench.The Celebes oceanic lithosphere is subducting beneath the Mindanao and Sulawesi islands along the Cotabato and Sulawesi trenches, which are located in the west of the Philippine trench.The Molucca Sea plate has subducted beneath the Sulawesi-Sangihe and Halmahera islands in the southern part of the Philippine islands along the Sangihe and Halmahera trenches, resulting in an arc-arc collision and double-side subduction system.In southern Indonesia, the oceanic lithosphere of the Indo-Australian plate is subducting beneath the Eurasian plate, generating a long volcanic arc that extends from Sumatra to the Flores island (Miller et al., 2016;Titu-Eki & Hall, 2020).The Indo-Australian plate transforms into the continental lithosphere of Australia as it collides with the Banda arc, known for its 180° curvature, near the east of the Flores Island (Miller et al., 2016).The shape of the volcanic arc and the deformed slab has been attributed to the interactions between the PHS plate, the Eurasian plate, and the Indo-Australian plate.The subduction and collision of the three plates have produced a unique tectonic system with frequent earthquakes, arc volcanoes, islands, and microcontinental blocks.The tectonics in this region is rendered particularly complex due to the subduction and interactions of multiple slabs in the mantle.Hence, the Philippines-Banda area is an ideal natural laboratory to study the mantle flow pattern associated with the slabslab interactions.
Seismic anisotropy is a phenomenon that seismic wave velocity changes with its polarization or propagation direction, due to lattice-preferred orientation (LPO) and shaped-preferred orientation (SPO) of minerals (e.g., Long & Becker, 2010;Savage, 1999;D. Zhao et al., 2023).In the crust, orientations of local geological structures and tectonics can cause seismic anisotropy, such as active faults, folds, fracture zones, and mountain ranges.In the upper mantle, seismic anisotropy is mainly caused by the LPO of olivine crystals induced by mantle flow (Savage et al., 2016).Hence, investigating seismic anisotropy can provide critical information on mantle flow.Previous studies have investigated seismic anisotropy beneath the Philippine-Banda area by making shear wave splitting measurements (Cao et al., 2021;Di Leo et al., 2012;Hammond et al., 2010;L. Wang & He, 2020).However, the shear wave splitting method is hard to determine the source depth of anisotropy, especially in the presence of complex slab geometry and mantle flow patterns induced by the slab-slab interaction beneath a study region.
In this work we investigate the 3-D azimuthal anisotropy structure of the Banda area by jointly inverting a large amount of P and S wave arrival times of local earthquakes.The obtained high-resolution 3-D images of isotropic velocity and azimuthal anisotropy shed new light on mantle flows induced by the multiple plate subduction and slab-slab interactions.

Data and Methods
In this work we use P and S wave arrival time data of local earthquakes, which are collected from the International Seismological Center (ISC)-EHB Bulletins during 1999-2018 (Engdahl et al., 2020, http://www.isc.ac.uk/isc-ehb).Our data set includes 8,844 local earthquakes recorded at 175 seismic stations in the study region (Figure 1).All the selected events are located in the study region, and each event was recorded at over 10 seismic stations.P and S wave travel time residuals spread in an extensive range, and we take ±3 s as the cutoff residual for the P wave data and ±4 s as the cutoff residual for the S wave data in the final tomographic inversion (Figure S1 in Supporting Information S1).As a result, our final data set includes ∼200,000 P wave arrivals and ∼30,000 S wave arrivals from the 8,844 local earthquakes (Figure 1).These events are relocated during the tomographic inversion.However, for those events with large relocation errors, we adopt the well-determined hypocentral parameters of the ISC-EHB catalog so as to obtain a more robust tomographic model.
Most previous studies of anisotropic tomography used only P wave data to constrain the 3-D anisotropy structure (see a recent review by D. Zhao et al. (2023)).On the basis of the tomographic method of D. Zhao et al. (1992) and J. Wang and Zhao (2008), Liu and Zhao (2016b) developed an updated technique to constrain the 3-D anisotropic structure by jointly inverting P and S wave travel time data.Orthorhombic anisotropy with a horizontal symmetry plane is assumed to exist in the modeling space.Here we apply the technique of Liu and Zhao (2016b) to determine isotropic Vp and Vs tomography and 3-D azimuthal anisotropy beneath the study region.We set up two 3-D grids in the study volume: one fine grid for expressing the 3-D isotropic velocity structure, the other coarse grid for expressing the azimuthal anisotropy (Figure S2 in Supporting Information S1).Following the approach of Liu and Zhao (2016b), before conducting the anisotropic tomography, we determined the included angle between the polarization directions of a general S wave and an SV polarized wave with the same propagation direction for the S wave data set.The RMS (root-mean-square) S wave travel time residual becomes the minimum when the included angle is zero.Hence, in the final joint inversion of P and S wave data, the included angle is set to be zero.

Results and Resolution Analysis
Our tomographic results show that significant high-Vp and high-Vs anomalies exist along the Java-Banda arc at depths >50 km (Figure 2), which represent the subducting Indo-Australian slab.Intermediate-depth and deep earthquakes occurred actively within the high-velocity (high-V) slab in the Banda arc with a 180° curvature.Significant low-velocity (low-V) anomalies exist in the bay area of the Banda arc, which represent hot upwelling flow induced by the slab subduction.The fast velocity directions (FVDs) of azimuthal anisotropy within these low-V zones seem to be parallel with the strike of the Banda arc.To the north of the Banda arc, the subducting Molucca Sea slab is imaged as a significant high-V anomaly at depths >50 km.Trench-normal FVDs appear within the high-V Molucca Sea slab.At depths <75 km, significant low-V anomalies exist beneath the arc volcanoes (Figure 2).These features are similar to the tomographic results obtained by inverting only P or S wave data (Figure S4 in Supporting Information S1).The FVDs to the east of Sulawesi island are nearly N-S, whereas Wei et al. (2022) showed NE-SW FVDs there.This difference could be due to several factors, such as differences in the data set, range of the study area, model parameters, and damping and smoothing parameters, etc.  of the study volume where the rays crisscross very well.To perform the RRTs, the obtained azimuthal anisotropic results (Figure 2) are adopted as the input model.The RRT results (Figures S27-S34 in Supporting Information S1) show that the main features of the tomographic images are well recovered, suggesting that the trade-off is insignificant and the results of azimuthal anisotropy tomography are quite reliable.We also investigated the influence of the P and S wave data weighting on the inversion results.The results (Figure S35 in Supporting Information S1) show that increasing the weight of S-wave data leads to slightly larger amplitude of anisotropy, but the FVDs have almost no change, indicating that the different weights of the P and S wave data have little effect on the anisotropic result.

Slab Geometries
The geometry of subducting slabs is essential for understanding slab subduction and its dynamic processes.Since almost all the intermediate-depth and deep earthquakes (at depths of 60-670 km) occur within the subducting slab, the deep seismic zone can somewhat represent the subducting slab geometry (e.g., Hayes et al., 2018).Due to the unique tectonic setting of the study area, a large number of seismic events, including deep earthquakes, have been recorded in this region.The distribution of these earthquakes provides a preliminary insight into the complex geometry of the slabs, such as the arc shape of the subducting Indo-Australian slab and the double-side subduction of the Molucca Sea plate (e.g., Cardwell & Isacks, 1978;Hutchings & Mooney, 2021).
Our results show clearly arc-shaped high-V anomalies at depths of ∼50-400 km below the Banda arc (Figures 2  and 3), representing the subducting Indo-Australian slab, that is, the Banda slab.Taking into account previous large-scale tomographic images (e.g., Hua et al., 2022;Huang et al., 2015;Spakman & Hall, 2010;Wei et al., 2022;Widiyantoro & Hilst, 1997), we suggest that the subducted slab is trapped within the mantle transition zone (MTZ).Note that only local earthquake arrival-time data are used in this study, so our model depth is limited by the distribution of local deep earthquakes.Our model is generally consistent with the results of full waveform tomography (Fichtner et al., 2010), which show a high-V anomaly beneath the Banda arc.Full waveform tomographic models suggest that the western edge of the subducted continental lithosphere is located near central Flores (Fichtner et al., 2010;Shulgin et al., 2009), as shown also by the isotope ratios and abundances of He, Sr, and Nd in the intrusive and extrusive igneous samples (Elburg et al., 2004(Elburg et al., , 2005)).An interesting feature is the reversal of the polarity of slab subduction at the eastern end of the Banda arc.The slab subducts northward beneath the Timor trough and southward beneath the Seram trench.The two slab-related high-V anomalies contact at a depth of ∼450 km and appear in the MTZ, forming a unique structure with a spoon shape (Wehner et al., 2022;Wei et al., 2022;Widiyantoro & Hilst, 1997), which is also visible in the deep seismicity (Cardwell & Isacks, 1978;Hutchings & Mooney, 2021).Wei et al. (2022) used regional earthquake data to reveal the stagnant slab in the MTZ, but our model has low resolution in the MTZ and so could not detect the stagnant slab due to the limited ray coverage at depths greater than 400 km.
Another interesting feature is the double-side subduction of the Molucca Sea plate.Despite the limited resolution, our model shows an apparent high-V anomaly beneath the Molucca Sea, which dips westward beneath Sangihe and eastward beneath Halmahera (Figure 3).This feature is consistent with previous tomographic results (Fan & Zhao, 2018;Wei et al., 2022) and the distribution of seismicity (Cardwell & Isacks, 1978;Hutchings & Mooney, 2021).The high-V anomaly below Sangihe continues to the MTZ, while the high-V anomaly below Halmahera seems to extend only to ∼250-300 km depths (Figure 3).In other words, this double-side subduction system and the arc-arc collision system are strongly asymmetric.This difference may be caused by the different subduction histories and subduction rates beneath Sangihe and Halmahera (Király et al., 2018).The Sangihe trench has been active since ∼20-30 Ma, producing at least 1,000 km length of plate convergence and the deepest slab in the Philippine mobile belt (Lallemand et al., 1998).In contrast, the subduction beneath the Halmahera trench started ∼10 Myr ago and produced ∼300 km long plate convergence with an average convergence rate of ∼3 cm/yr (Lallemand et al., 1998).Similar asymmetric double-side subduction has been also found in Europe, where the Adria Sea slab dips toward the west beneath the Apennine trench and dips toward the east beneath the Dinaride trench (Hua et al., 2017;Király et al., 2018).Numerical simulations suggested that different plate tensions between the two sides of subduction caused the plate to move to the side with a faster convergence rate, generating the asymmetry of the double-side subduction (Király et al., 2018).

Anisotropy Within the Subducting Slab
Our tomographic model shows an apparent high-V anomaly along the Banda arc (Figure 2), representing the subducted Banda slab with abundant deep seismicity.In the northern part of the Banda arc, the subducted high-V slab dips to the south.The subducted high-V slab dips to the north in the southern part of the Banda arc.The two parts of the subducting slab are interconnected and stagnant in the MTZ, as shown by the deep seismicity.
Beneath the Sumba island, the oceanic lithosphere of the Indo-Australian plate is subducting beneath the Eurasian plate (Harris et al., 2020).The FVDs within the subducted slab are NW-SE, which are perpendicular to the isochrons of the plate.Similar features have also been detected in the subducting Pacific slab beneath Alaska (Gou et al., 2019;Tian & Zhao, 2012) and the subducting Cocos slab beneath Costa Rita (Rabbel et al., 2011).We think that the observed NW-SE FVDs beneath the Java and Sumba islands represent the fossil anisotropy acquired by the plate at the mid-ocean ridge.
The slab boundary in the south of the Sumba Island represents a change in the tectonic regime from the Indo-Australian oceanic lithospheric subduction to the continental-island arc collision.To the east of the Sumba Island, the FVDs in the high-V anomalies from Timor to below the Babar Island are predominantly trench-parallel, which may represent the fossil anisotropy of the continental lithosphere.Below the Weber and Seram Troughs, the subducting plate bends with a 180° curvature (Figure 2), which is consistent with the background seismicity.
Focal mechanism solutions depict that the slab beneath the contorted eastern Banda arc has undergone lateral compression, probably due to the slab bending (Cardwell & Isacks, 1978).The FVDs within the high-V slab beneath the Banda Sea at 250 km depth are perpendicular to the trench (Figure 2).A previous study has suggested that the formation of the 180° curved Banda volcanic arc was caused by a slab being extruded and then bent by an external force (Bowin et al., 1980).External forces or thermal events may change the LPO of minerals and modify the fossil anisotropy to cause such trench-normal FVDs.Our observed FVDs of anisotropy could be produced by this process.

Mantle Flow in the Banda Subduction Zone
Due to the coupling between the subducting slab and its surrounding mantle, different styles of mantle convection may exist during the slab subduction.There are two simplified mantle flow endmember models (Figure 4a; Long, 2013).The first model is about two-dimensional corner flow generated by the coupling between the subducting slab and its surrounding mantle.This flow pattern has been detected by seismic anisotropy tomography in Japan (Liu & Zhao, 2016b;J. Wang & Zhao, 2008), Alaska (Gou et al., 2019;Tian & Zhao, 2012), Cascadia (Eakin et al., 2010;D. Zhao & Hua, 2021), and Chile (Huang et al., 2019).The second model is about 3-D mantle convection caused by trench migration and slab rollback, which generate particular circulation at the slab edge or slab window (D.Zhao & Hua, 2021).Both models are highly simplified, and mantle convection in an actual subduction zone is usually explained by a combination of the two models.In particular, the Banda-Philippines region is tectonically complex due to the existence of multiple subducting plates and the spoon-shaped slab, so the two end-member models cannot explicitly explain its mantle flow.The curvature of the slab and the interaction between the adjacent slabs must have affected the mantle flow pattern in the subduction zone.
The most striking tectonic feature of the Banda subduction zone is the 180° curved volcanic arc.Previous studies have suggested that this volcanic arc was caused by the opposite subduction of two plates or by bending of one plate.Either way, the particular curvature of this subducting plate would produce a unique pattern of mantle convection.Di Leo et al. (2012) showed that FVDs in the Banda arc were trench-parallel, and they suggested that the oblique supraslab shear layers caused the FVDs in the Banda mantle wedge.Our 3-D azimuthal anisotropy tomography shows trench-parallel FVDs in the significant low-V anomalies at depths >100 km beneath the Banda Sea (Figure 2).The FVDs are usually parallel to the mantle flow direction in the A-type olivine fabric.Due to the subduction of the curved Banda slab or the opposite subduction of two slabs, asthenospheric materials may be squeezed out because of the limited space in the mantle wedge (L.Wang & He, 2020).The trench-parallel FVDs in the north of the Banda arc at depths >100 km may represent the lateral mantle flow caused by the subduction of the highly curved Banda slab.These characteristics are different from those in other subduction zones, such as Japan (Liu & Zhao, 2016a, 2016b), Cascadia (D. Zhao & Hua, 2021), and Alaska (Gou et al., 2019;Tian & Zhao, 2012), which show trench-normal FVDs in the mantle wedge under back-arc areas.This difference may be caused by the spoon shape of the subducting Indo-Australian slab.The big curvature of the Indo-Australian slab leads to a unique mantle flow pattern, which is quite different from the two end-member models (Figure 4).
In the north of the Banda slab and the east of the Molucca Sea slab, nearly NE-SW FVDs occur in the low-V anomaly (Figure 2).Around the northern edge of the Banda slab, our anisotropic tomography shows a semi-toroidal pattern of FVDs.In addition, the Vp and Vs images show a certain gap between the Banda slab and the Molucca Sea slab (Figure 2).Considering the mantle flow pattern of the Banda arc as mentioned above, we think that the semi-toroidal pattern of FVDs may represent mantle flow extruded from the Banda arc and escaped from the middle gap of the Banda-Molucca slab toward the northeast (Figure 4).The complexity of the slab shape could lead to toroidal mantle flow around the slab edge, as suggested by numerical simulations (Husson et al., 2022).This kind of toroidal mantle flow around the slab edge has also been revealed in other subduction zones, such as Cascadia (Eakin et al., 2010;D. Zhao & Hua, 2021) and Kamchatka (Jiang et al., 2009;L. Zhao et al., 2021).Our 3-D anisotropic tomography actually reveals the mantle flow pattern induced by the highly curved slab and slab-slab interactions (Figure 4), providing new insight into mantle convection in subduction zones.

Conclusions
The 3-D anisotropic structure of the Banda subduction zone is investigated by jointly inverting a large number of P and S wave arrival time data of local earthquakes.Our results shed new light on the mantle flow pattern affected by the curvature of the subducting slab and slab-slab interactions.The main findings of this study are summarized as follows.
1.The first 3-D model of the P and S wave azimuthal anisotropy down to 500 km depth beneath the Banda region is obtained.2. The subducting Indo-Australian slab is clearly revealed along the curved Banda arc.The fast-velocity direction (FVD) of azimuthal anisotropy in the slab changes along the arc strike, representing fossil anisotropy within the slab modified by the subduction processes.3. Low-velocity anomalies with trench-parallel FVDs along the Banda arc are revealed, which reflect lateral mantle flow caused by subduction of the highly curved Banda slab.4. A semi-toroidal pattern of FVDs appears around the northern edge of the Banda arc, which may represent mantle flow extruded from the Banda arc.We appreciate helpful discussions with Drs.Fan Yang and Yang Yu.This work was supported by research grants from the National Natural Science Foundation of China (41890812, 42106066), Key Special Project for Introduced Talents Team of Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou) (GML2019ZD0202), Tuguangchi Award for Excellent Young Scholar (SKLaBIG-TJ-22-01), Japan Society for the Promotion of Science (19H01996), the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan under its The Second Earthquake and Volcano Hazards Observation and Research Program (Earthquake and Volcano Hazard Reduction Research, THK-05), and the "14th Five-Year Plan" Independent Project of GIGCAS.Prof. Lucy Flesch (the Editor) and two anonymous referees provided thoughtful review comments and suggestions, which have improved this paper.

Figure 1 .
Figure 1.(a) Tectonic settings of the study region.The red triangles denote active volcanoes.The two red lines denote locations of vertical cross-sections in Figure 3.The slab age data are from Müller et al. (2008).(b) Epicentral distribution of the 8,844 local earthquakes (color dots) used in this study.The dot colors denote the focal depths whose scale is shown at the upper-right corner.(c) Distribution of the 175 seismic stations (blue squares) used in this study.The red box denotes the target area of this study.

Figure 2 .
Figure 2. Map views of the tomographic results.The top row shows isotropic P-wave velocity (Vp) tomography at depths of 25, 150, 250, and 430 km.The second row shows isotropic S-wave velocity (Vs) tomography.The third and fourth rows show Vp and Vs azimuthal anisotropy tomography obtained by a joint inversion of P and S wave arrival times.The blue and red colors denote high and low velocity perturbations, respectively.The orientation and length of the white bars represent the fast-velocity directions (FVDs) and the amplitude of azimuthal anisotropy, respectively.The scales of the isotropic velocity perturbation and the anisotropic amplitude are shown at the bottom.

Figure 3 .
Figure 3. Vertical cross-sections of isotropic Vp and Vs tomography along the two profiles as shown in Figure 1a.The red and blue colors denote low and high velocity perturbations, respectively, whose scales are shown on the right.The open circles denote background seismicity within a width of 30 km of each profile.The red triangles denote active volcanos.The blue patches above each panel denote the surface topography.The black dashed line denotes the 410 km discontinuity.

Figure 4 .
Figure 4. 3-D views of the subducting slab and mantle flows.(a) Two endmember models of mantle flow in a normal subduction zone.(b) Trench-parallel flow in the mantle wedge of the curved Banda subduction zone.