Azimuthal Anisotropy Tomography of the Southeast Asia Subduction System

A detailed 3‐D model of P wave isotropic velocity and azimuthal anisotropy beneath SE Asia is obtained by jointly inverting local earthquake arrival times and teleseismic relative travel‐time residuals. Our results show that the high‐velocity (high‐V) subducting Australian slab has penetrated through the mantle transition zone (MTZ) and reached a depth of ∼1,200 km beneath Sumatra, whereas the high‐V slab has subducted toward the north and trapped within the MTZ beneath Java. The Hainan mantle plume is revealed clearly as a significant low‐velocity anomaly beneath the MTZ, which extends down to the bottom of our model (∼1,600‐km depth). The upwelling of the Hainan plume is resisted by the stagnant slab in the MTZ, resulting in divergent fast‐velocity directions of azimuthal anisotropy beneath the stagnant slab. The detailed results of seismic tomography and azimuthal anisotropy provide important new information on the complex mantle structure and dynamics of the SE Asian region, in particular, the subducting slabs and the Hainan plume, as well as their interactions.

Isotropic seismic tomography conducted so far provides a snapshot of the present-day structure of the Earth but does not directly constrain the ongoing dynamic processes (Zhao, 2015). The velocity of a seismic wave in the Earth usually changes with its traveling direction. This phenomenon is called seismic anisotropy (e.g., Long, 2013), and it widely exists in the Earth's interior, due to lattice-preferred orientation (LPO) and shaped-preferred orientation (SPO) of minerals (e.g., Long & Becker, 2010;Savage, 1999;Silver, 1996;Zhao et al., 2016). Mineral fabrics in the mantle can be preserved for a long time if no thermal and tectonic events affect the preexisting fabric (e.g., Ben-Ismail et al., 2001). Hence, seismic anisotropy is a powerful tool for probing mantle dynamics, because it directly links to mantle deformation (Long, 2013). Measuring shear wave splitting (SWS) has become a popular tool for detecting anisotropy in the crust and mantle. In the Sumatra subduction zone, trench-parallel fast-velocity direction (FVD) was observed using local SWS measurements in the fore-arc area, which may be associated with SPO of cracks/fractures in the overriding sediments (Collings et al., 2013;Hammond et al., 2010). The FVD of the SKS splitting is roughly parallel to the absolute plate motion (APM) direction of the subducting Indo-Australian Plate, suggesting that entrained flow may exist in the asthenosphere below the slab (Collings et al., 2013;Kong, Gao, Liu, Zhang, et al., 2020). Source-side SWS observations show that the APM-parallel FVD in Sumatra changes to roughly trench-parallel in Java (Lynner & Long, 2014). This transition occurs in the oldest part of the subducting slab and is consistent with the SKS splitting result, which may be caused by the different slab ages beneath Sumatra and Java (Hammond et al., 2010). The younger oceanic lithosphere beneath Sumatra may cause strong coupling between the subducting slab and the surrounding mantle, whereas the older oceanic lithosphere beneath Java may cause weak coupling or decoupling between them, resulting in differences in mantle flow pattern and anisotropy (Lynner & Long, 2014). APM-parallel FVDs also exist beneath northern Borneo, the Nansha Block, and the Indochina Peninsula, which may be associated with the partial coupling between the lithosphere and asthenosphere (Song et al., 2021;Yu et al., 2018). A disadvantage of the SWS technique is that it is a path-integrated measurement and so has low depth resolution even when multiple phases are used (Long, 2013). This problem can be resolved by 3-D P wave velocity (Vp) anisotropic tomography that includes parameters describing anisotropy in the velocity inversion .

10.1029/2021JB022854
3 of 21 In southern Sumatra, Liu et al. (2021) determined a local model of Vp anisotropic tomography down to ∼100-km depth, which shows trench-normal FVDs prevailing in the lower crust and trench-parallel FVDs in the fore-arc and back-arc mantle wedge. Using P wave arrival-time data of local earthquakes, Huang et al. (2015) determined a 3-D model of Vp anisotropic tomography down to 900-km depth beneath SE Asia, which shows significant low-V mantle flow beneath SE Tibet, Hainan, and Vietnam. However, the two previous Vp models have limited resolution because they did not use teleseismic data and adopted strong smoothing and damping constraints in their tomographic inversions. Due to the limited distribution of seismic stations and raypaths, the detailed 3-D anisotropic structure of the SE Asia region is still unclear, especially beneath the MTZ. In the present study, we investigate the detailed 3-D Vp anisotropic structure down to 1,600-km depth beneath the whole SE Asia subduction system by jointly inverting over 1 million arrival times of local earthquakes and relative travel-time residuals of teleseismic events. Our model of anisotropic tomography sheds new light on the complex tectonics and subduction dynamics beneath SE Asia, in particular, the slab-MTZ-plume interactions.

Data and Methods
We use abundant P wave arrival times of local earthquakes and relative travel-time residuals of teleseismic events (Figures 2 and 3). These events occurred during 1990-2019 and were recorded at 807 seismic stations (Figure 4). The P wave arrival times were generated by 19,525 local earthquakes ( Figure 2) and 5,688 teleseismic events (Figure 3), which are collected from the bulletins of the International Seismological Center (ISC). All the local and teleseismic events were relocated precisely using multiple P and S waves including depth phases based on the ISC location algorithm (Engdahl et al., 1998). Each of the selected local earthquakes was recorded at over . Tectonic settings of the study region. The colors show the seafloor age whose scale is shown at the lower-left corner. The age data are from Müller et al. (2008). The red triangles denote active volcanoes. HN, the Hainan volcano. TC, the Tengchong volcano.

of 21
30 seismic stations in the study region ( Figure 4). To make these events distribute as uniformly as possible in the study volume, we adopted an event selection scheme to divide the study volume into small cubic blocks (10 km × 10 km × 5 km). In each block, we select one earthquake that has the maximum number of arrival-time data. Travel-time residuals of the local earthquakes spread in a range of ± 8 s, and we take ± 4 s as the cutoff residual in the final tomographic inversion ( Figure S1 in Supporting Information S1). As a result, about 1.05 million P wave arrival times of 19,525 local earthquakes recorded at 807 stations ( Figure 4) are selected for the final tomographic inversion. Although most of these local events occurred at depths <100 km, many deep events with focal depths >400 km exist in the Java-Banda subduction zone ( Figure 2).
The epicentral distances of the 5,688 teleseismic events selected are limited in a range of 30°-90° ( Figure 3). Each of the selected teleseismic events was recorded at more than 30 seismic stations in the study region (Figure 4). To avoid the effect of hypocentral parameter uncertainties and structural heterogeneities that exist outside of the study volume, relative travel-time residuals are determined and used in the teleseismic tomographic inversion (Zhao, 2015;Zhao et al., 1994). The relative travel-time residuals of the teleseismic events spread in a range of ± 8 s, and we take ± 4 s as the cutoff residuals in the tomographic inversion ( Figure S1 in Supporting Information S1). As a result, more than 320,000 teleseismic data are selected to conduct the final tomographic An isotropic tomography method (Zhao, 2015;Zhao et al., 1994) is used to investigate the 3-D isotropic Vp structure of the study volume. A boundary-grid approach is adopted in the model parameterization. The velocity discontinuities (i.e., the Conrad, the Moho, the 410, and 660 km discontinuities) and the station elevations are taken into account in the 3-D ray tracing and the tomographic inversion. In the present study, lateral depth variations of the Conrad and Moho discontinuities are derived from the CRUST1.0 model (Laske et al., 2013). The starting 1-D velocity model used in the 3-D tomographic inversion was modified from the AK135 velocity model (Kennett & Engdahl, 1991). A 3-D grid with a lateral grid interval of 1.0° is set up in the study volume to express the 3-D Vp structure ( Figure S3 in Supporting Information S1). The grid nodes are set at depths of 15,50,100,150,200,250,330,430,530,630,750,900,1,100,1,300, and 1,600 km. Vp perturbations at the grid nodes relative to the 1-D starting velocity model are taken as unknowns to be inverted. An efficient 3-D ray tracing scheme (Zhao et al., 1994) is used to compute theoretical raypaths and travel times. The system of observation equations is usually very large but sparse. Here, we adopt the L-BFGS-B algorithm for bound constrained optimization to solve the system of observation equations (Morales & Nocedal, 2011;Zhu et al., 1997), instead of the LSQR algorithm (Paige & Saunders, 1982) that was usually adopted by previous tomographic studies (e.g., Wang & Zhao, 2008, 2013Zhao et al., 1994). This L-BFGS-B inversion scheme does not require damping regularization as the LSQR algorithm does, hence the inversion results are less influenced by human subjectivity . Smoothing is adopted in the tomographic inversion to avoid sharp variations of the 3-D velocity model. The optimal value of the smoothing parameter is determined after conducting many tomographic inversions with different values of the parameter. The optimal smoothing regularization can greatly suppress small and unstable Vp anomalies but keep important structural details that can be solved with good resolution. The root-mean-square (RMS) travel-time residual is 1.69 s before the tomographic inversion and reduced to 1.35 s after the isotropic inversion. The method of Zhao (2008, 2019) is used to conduct azimuthal anisotropy tomographic inversions. We assume that weak Vp azimuthal anisotropy with horizontal hexagonal symmetry exists in the modeling space (Becker et al., 2006;Browaeys & Chevrot, 2004;Liu & Zhao, 2016b). For the azimuthal anisotropy tomographic inversion, at each grid node there are two anisotropic parameters expressing the anisotropic structure, in addition to an isotropic Vp parameter. Thus, to obtain a robust model of 3-D azimuthal anisotropy structure, the ray coverage should be more uniform and denser than that for the 3-D isotropic Vp structure. A coarser 3-D grid with a horizontal grid interval of 2.0 • is set up in the study volume to express Vp azimuthal anisotropy ( Figure S3 in Supporting Information S1). Following Wang and Zhao (2019), we use the L-BFGS-B algorithm to conduct the inversion for Vp azimuthal anisotropy. Smoothing regularization is adopted in the tomographic inversion to avoid sharp changes in the anisotropic model. The P wave RMS travel-time residual is 1.29 s after the Vp anisotropic inversion, which is smaller than that (1.35 s) after the isotropic Vp inversion. To know whether the residual reduction is statistically significant or not, we conducted F test (Hua et al., 2017;Zervas & Crosson, 1986;Zhao et al., 1995). The test result shows that the F ratio for the isotropic and anisotropic Vp models is 8.78, which is much greater than the corresponding value (1.0) listed in the mathematical table, indicating that the RMS residual reduction from the isotropic inversion to the anisotropic inversion is indeed statistically significant.

Isotropic Tomography
The local and teleseismic raypaths crisscross quite well in our study volume ( Figure S2 in Supporting Information S1), enabling us to obtain a robust model of isotropic and anisotropic Vp tomography beneath SE Asia. The spatial resolution of our tomographic model is carefully evaluated. The checkerboard resolution test (CRT) is usually used to evaluate the model spatial resolution. We conducted two CRTs with different lateral grid intervals of 1.0° and 2.0°.
In the input velocity model, we assign ± 4% isotropic Vp perturbations at adjacent grid nodes. Synthetic travel times are calculated for the checkerboard model, and Gaussian noise (−0.25 to +0.25 s) with a standard deviation of 0.1 s is added to the synthetic data. Then we invert the synthetic travel times to see whether the input velocity model can be well recovered or not. The obtained CRT results show that most of the input Vp perturbations can be well recovered (Figures S4 and S5 in Supporting Information S1), except for the South China Sea (SCS) at shallow depths (<100 km) where the raypath coverage is not good enough. To evaluate whether and how much the noise amplitude affects the CRT results, we performed three more CRTs with larger standard deviations (0.2-0.4 s) of Gaussian noise. The test results ( Figure S6 in Supporting Information S1) show that the input checkerboard model can be well recovered even if the stronger noise is added to the synthetic data.
To confirm the main features of our tomographic model, we also performed restoring resolution tests (RRTs). The RRT procedure is similar to that of the CRT. The obtained isotropic tomography results are modified to be the input model of RRT. The RRT results show that main high and low Vp perturbations in the input velocity model are well recovered ( Figure S7 in Supporting Information S1). Another RRT is conducted with an input model in which the obtained velocity anomalies are reversed, and the results show that the main features of the input velocity model are well reconstructed ( Figure S8 in Supporting Information S1). To evaluate whether the Vp anomaly beneath the Hainan area is robust or not, we conducted a synthetic test, which shows that the Hainan plume can be well recovered, if it exists ( Figure S9 in Supporting Information S1).
Our data were recorded in a long period from 1990 to 2019. To investigate whether the different periods of the data would affect the tomography result, we divided our data set into two subdata sets ( Figure S10 in Supporting Information S1). The first subdata set contains the local and teleseismic events that occurred during 1990-2006, and the second subdata set contains the local and teleseismic events that occurred during 2007-2019. Then we conducted tomographic inversions of the two subdata sets. The obtained results ( Figure S11 in Supporting Information S1) show that main features of Vp anomalies derived from the two subdata sets are quite consistent with each other, though there are some minor differences between them, which are caused by the different distributions of earthquakes and seismic stations in the two data sets. These results indicate that the effect of the different periods of the data is very small.
We used the bootstrap method to estimate uncertainties of the 3-D isotropic Vp model. We constructed 50 subdata sets from the original data set. Each subdata set contains 30% of earthquakes that are selected randomly from the original data set. By conducting tomographic inversions of the 50 subdata sets, we obtained 50 isotropic Vp models. The variation range of Vp perturbation (dVp) at each grid node is estimated from the 50 models and is taken to represent uncertainties of the Vp model ( Figure 6). The dVp uncertainty is smaller than 0.4% in most parts of the study volume, except for the edge parts where the ray coverage is not good enough.
Our isotropic Vp tomography reveals significant low-V anomalies beneath the volcanic front of the Sumatra-Java and the Philippine subduction zones, which may represent hot upwellings induced by slab subduction and dehydration (Figures 7-9). This feature has been detected by previous tomographic studies made in SE Asia (e.g., Fan Huang et al., 2015;Koulakov et al., 2016;Liu et al., 2019;Pesicek et al., 2008) and other subduction zones (e.g., Gou et al., 2019;Huang et al., 2011;Zhao et al., 1994). Another important feature of our model is that significant high-V anomalies exist beneath the Sumatra-Java subduction zone, which are interpreted as the subducting Indo-Australian slab (Huang et   . Distribution of isotropic Vp perturbation uncertainties (in %). The layer depth is shown at the lower-left corner of each map. Blue and red colors denote low and high uncertainties, respectively, whose scale is shown at the bottom. Liu et al., 2019). Beneath Sumatra the high-V slab is subducting toward the northeast with a near-vertical dip angle (BBʹ and CCʹ profiles in Figure 8). This high-V slab has penetrated through the MTZ and reached a depth of ∼1,200 km. Beneath Java the high-V slab is subducting toward the north and trapped within the MTZ (DDʹ Figure 7. Map views of the isotropic P wave velocity (Vp) tomography. The layer depth is shown at the lower-left corner of each map. Red and blue colors denote low and high Vp perturbations, respectively, whose scale is shown at the bottom. Prominent high Vp anomalies exist in the MTZ (410-660-km depths), which represent the stagnant slab. Significant low Vp anomalies exist beneath the MTZ, reflecting the Hainan plume ascending from the lower mantle. and GGʹ profiles in Figure 8). Our isotropic tomographic model suggests that the slab geometry under Sumatra is different from that under Java. The slab beneath Sumatra has a near-vertical dip angle and penetrates across the MTZ, while the slab beneath Java has a relatively small dip angle and is trapped within the MTZ (Figure 8). This transition takes place between the CCʹ and DDʹ profiles, corresponding to the isochrons of ∼90-95 Ma. The differences in the slab geometry may be associated with the age difference of the subducting slab beneath Sumatra and Java (Figure 1). Beneath the MTZ, a low-V zone extends from 660-km to 1,600-km depth, which may reflect the Hainan plume (Huang, 2014;Lebedev & Nolet, 2003;Montelli et al., 2004;Xia et al., 2016;Zhao, 2004).

Azimuthal Anisotropy Tomography
As mentioned above, to determine a reliable 3-D model of Vp azimuthal anisotropy tomography, a uniform azimuthal coverage of raypaths is needed in the study volume. To further evaluate the ray coverage (e.g., ray density and azimuthal coverage) and the spatial resolution of our model, we performed CRT tests. In the input checkerboard model, besides the isotropic Vp perturbations of ± 4%, we also assign a perturbation of ± 3% to Figure 8. Vertical cross-sections of isotropic Vp tomography along the profiles as shown on the inset maps. The red and blue colors denote low and high Vp perturbations, respectively, whose scale is shown below the maps. The black circles denote background seismicity within a width of 50 km of each profile. The red triangles denote active volcanos. The two dashed lines denote the 410 and 660 km discontinuities. The slab beneath Sumatra has a near-vertical dip angle and penetrates across the MTZ, while the slab beneath Java has a relatively small dip angle and is trapped within the MTZ. the two azimuthal anisotropic parameters at adjacent grid nodes for Vp anisotropy. The ± 3% perturbations of the two anisotropic parameters result in azimuthal anisotropies with an amplitude of 4.25% and FVDs of 112.5 • or 22.5 • . The CRT procedure of anisotropic tomography is the same as that of isotropic tomography. The CRT results (Figures S12 and S13 in Supporting Information S1) show that the input isotropic velocity and azimuthal anisotropy model with a lateral grid interval of 2.0° can be well recovered, except for the SCS at shallow depths (<100 km), suggesting that the ray coverage in the study volume is good enough to obtain a robust azimuthal anisotropy model.
The trade-off between isotropic Vp and azimuthal anisotropy must be evaluated before interpreting the obtained anisotropic model. To do this, we assign perturbations of ±3% to the two azimuthal anisotropic parameters and zero perturbation to the isotropic Vp at adjacent grid nodes in the input model to compute the synthetic data. The test results (Figures S14 and S15 in Supporting Information S1) suggest that the trade-off between isotropic Vp and azimuthal anisotropy does exist but it is not significant. To further evaluate the robustness of our Vp anisotropic tomography, we performed extensive RRTs and synthetic tests with different input isotropic and anisotropic models. Results of the RRTs and synthetic tests show that the trade-off exists but is not obvious, and the obtained model of azimuthal anisotropy tomography is quite reliable (Figures S16-S21 in Supporting Information S1). Figure 10 shows map views of Vp azimuthal anisotropy tomography obtained by this study. The isotropic Vp images determined by the anisotropic tomography are very close to those obtained by the isotropic tomography (Figure 7). At depths <400 km, the FVDs are NE-SW in the Sumatra fore-arc area, which are nearly normal to the trench and parallel with the absolute plate motion direction. This feature is similar to the previous SKS splitting results (Hammond et al., 2010;Lynner & Long, 2014). The FVDs in the lower mantle are roughly trench-parallel. In the back-arc areas, the FVDs within the low-V anomalies are trench-normal at shallow depths (<100 km) and trench-parallel at greater depths (150-250 km). In the Java subduction zone, the FVDs are NW-SE at shallower depths (<250 km) and nearly E-W at greater depths (>400 km). The previous SKS results show that the FVDs in the Java subduction zone are generally trench-parallel and nearly E-W (Collings et al., 2013;Hammond et al., 2010;Lynner & Long, 2014). Wang and He (2020) showed that NE-SE FVDs exist in the Java back-arc area. The FVDs in the region from Tengchong to Hainan are mainly NW-SE, which may represent the extruded mantle materials from the Tibetan Plateau (Huang et al., 2015).

Aseismic Deep Slab Beneath Sumatra
As mentioned above, the Wadati-Benioff deep seismic zone is usually used to represent the general geometry of a subducting slab (e.g., Hasegawa et al., 1978;Hayes et al., 2018;Umino & Hasegawa, 1975). Beneath the Sumatra Island, earthquakes within the subducting slab generally occurred at depths <300 km (Figure 2). Using teleseismic tomography, Li et al. (2018) detected the subducting slab that reaches ∼500-km depth beneath Sumatra. A recent tomographic model obtained with local and teleseismic data shows that the subducted slab reaches ∼800-km depth in northern Sumatra (Liu et al., 2019). Regional and global tomographic models show that the high-V subducting slab has reached at least 1,000-km depth (Pesicek et al., 2008;Widiyantoro & van der Hilst, 1997). Although there is still no consensus on the slab geometry, an aseismic high-V slab subducting to a great depth (>500 km) exists indeed beneath Sumatra (Huang et al., 2015;Li et al., 2018;Liu et al., 2019;Pesicek et al., 2008;Widiyantoro & van der Hilst, 1997). In our model, a significant aseismic slab with a high-V anomaly is visible beneath the Sumatra Island and reaches ∼1,200-km depth (Figure 8). Our RRTs and synthetic tests suggest that the high-V anomaly beneath Sumatra is a quite robust feature (Figures S7-S9 in Supporting Information S1). Uplifting of the 410 km discontinuity and depression of the 660 km discontinuity detected by receiver functions suggest a colder MTZ beneath Sumatra (Kong, Gao, Liu, Ding, et al., 2020), which is consistent with our tomographic results.
The aseismic slab is an interesting feature and has been revealed in many subduction zones, such as Alaska (Gou et al., 2019;Qi et al., 2007), SW Japan (Asamori & Zhao, 2015;Huang et al., 2013), the Philippines (Fan & Zhao, 2019, and the Alps (Hua et al., 2017). The mechanism of intermediate-depth and deep earthquakes is controlled by temperature, pressure, and other factors (e.g., Jiang & Zhao, 2011;Jiang et al., 2015). Intermediate-depth and deep earthquakes take place within the subducting slabs because of low temperature and brittleness of the slabs. When the subducting slab is heated by the surrounding hot mantle, it will lose its brittleness and so intermediate-depth and deep earthquakes cannot be generated. However, the slab can continue to subduct downward in the mantle because of the gravity effect (Huang et al., 2013). The subducting slab is still colder than the surrounding hot mantle and manifests as a high-V anomaly in seismic tomography.
The slab geometry beneath Sumatra suggests that the subducting aseismic slab penetrates through the MTZ and slab roll back may occur in the lower mantle, whereas the high-V slab beneath Java has a relatively smaller dip angle and becomes stagnant in the MTZ. Goes et al. (2017) suggested that stagnant slabs usually exist in the western Pacific with old lithosphere age (90-145 m.y.), whereas some deeply penetrating slab sections are found along the eastern rim of the Pacific with relatively young lithosphere age (0-50 m.y.). Whether the subducting slab is stagnant in or penetrates through the MTZ is controlled by trench retreat, which is relevant with the strength and buoyancy of the subducting plate (van der Hilst & Seno, 1993). For a thermal aging trend of oceanic plates, both density and viscosity of the slab will increase with the slab age, retreating modes will be dominant, and the subducting slab is expected to stagnate in the MTZ. Numerical simulation models suggest that hot young slabs penetrate into the lower mantle more easily than cold old slabs (e.g., Agrusta et al., 2017). A plate reconstruction study (Müller et al., 2008) suggested that the slab beneath the Sumatra trench is younger (<95 Ma), whereas the slab beneath the Java trench is older (>100 Ma). The difference in the slab age may have caused the different slab behaviors (e.g., the slab dip angle and the maximum depth of slab seismicity) beneath Sumatra and Java.

Stagnant Slab in the Transition Zone
There are different styles of slab-MTZ interaction, one of which is the stagnant slab in the MTZ. The stagnant slab in the MTZ has been observed beneath subduction zones around the circum-Pacific and in Mediterranean (e.g., Fukao et al., 2009;Huang & Zhao, 2006;Lei & Zhao, 2005;Zhao, 2004). Another style of slab-MTZ interaction is that the subducted slabs pierce straight through the MTZ to the low mantle. Goes et al. (2017) suggested that there are two main factors that can affect whether slabs stagnate or penetrate: (a) the mantle resistance to flow through the MTZ and (b) the shape and strength of the slab as it starts interacting with the MTZ. Dynamic models suggested that the older lithosphere tends to lie down flat above the MTZ due to the high-viscosity lower mantle, whereas younger lithosphere can penetrate into the lower mantle more readily, buckles, and thickens. Slab thickening enhances the negative buoyancy, thus facilitating fast lower-mantle penetration (Goes et al., 2008). As mentioned above, the subducting older Australian slab is stagnant in the MTZ, whereas the younger Sumatra slab is penetrating through the MTZ to the lower mantle.
Our tomographic results show that significant high-V anomalies exist in the MTZ beneath SE Asia (Figure 7), which represent the stagnant slab in the MTZ. The subducted slab meets strong resistance at the base of the MTZ, leading to flattening, buckling, and thickening and further deformation . Wu and Suppe (2018) suggested that part of the high-V anomalies beneath the SCS represent the Proto-SCS slab. Significant MTZ thickening beneath the SCS southern margin has also provided evidence for the existence of the Proto-SCS slab . Considering the limited spatial resolution of our tomographic model, we cannot precisely estimate the real volumes of stagnant slab that came from the Australian plate collision. The period of a slab stagnation in the MTZ is unlikely to exceed 60 Ma . Assuming that the Australian plate keeps its rigidity after it subducted beneath the Sunda plate since 60 Ma with a velocity of 6 cm/year, the Australian slab has subducted into the mantle with a length of ∼3,600 km. In the H-Hʹ profile (Figure 8), the high-V anomalies in the MTZ may represent not only the stagnant Australian slab but also the subducting Philippines Sea slab, and part of the Proto-SCS plate (Wu & Suppe, 2018;Yu et al., 2021).
Our tomographic model shows that high-V anomalies exist in the MTZ beneath almost the entire SE Asia subduction zone. To show the real velocity structure, our 3-D Vp model keeps the background velocity at each depth, which is different from the previous tomographic models in which the mean Vp perturbation at each depth is removed to show the lateral Vp variations more clearly (e.g., Huang et al., 2015). In those models, high-V anomalies only exist in part of the MTZ beneath the SE Asia subduction zone and low-V anomalies under the Hainan volcano extend from the lower mantle through the MTZ to the upper mantle. Note that Huang et al. (2015) used a few more stations in Hainan than this work, which may cause some differences between the two models. In addition, removing the mean Vp perturbation at each depth or not may also generate some differences in the tomographic images.

The Hainan Plume
Surface wave tomography detected a significant low-V anomaly in the upper mantle beneath the Hainan Island (Lebedev & Nolet, 2003); however, its resolution is not high enough to reveal the detailed structure of the Hainan plume. Lei et al. (2009) conducted local tomography and revealed the Hainan plume that dips toward the southeast to 300-km depth with a diameter of ∼80 km. Regional body-wave tomography revealed an obvious low-V anomaly beneath the Leiqiong area down to 1,100-km depth, which reflects the Hainan mantle plume (Xia et al., 2016). This continuous low-V plume has a head in and around the MTZ beneath South China, which transforms into smaller low-V columns at shallow depths (Xia et al., 2016). Global tomography revealed a significant low-V anomaly under the Hainan Island down to the lower mantle (e.g., Montelli et al., 2004;Zhao, 2004;Zhao et al., 2021). Teleseismic receiver functions revealed a thinner MTZ beneath Leiqiong, suggesting that the Hainan plume has a lower-mantle origin (Wang & Huang, 2012;Wei & Chen, 2016). Although the previous tomographic studies revealed the Hainan plume with a significant low-V anomaly, its geometry and depth range detected by these studies are not consistent with each other, partly due to different data sets and methods used.
Local tomographic models usually have higher resolution in the upper mantle, but they are hard to resolve the 3-D structure of the MTZ and the lower mantle. The data set used in this study has fewer seismic stations in the Hainan Island, so our model has lower resolution at the shallow depths around the Hainan Island. Previous tomographic studies used arrival-time data recorded at seismic stations in and around the Hainan Island and showed 3-D structure of the Hainan plume at the shallow depths (Huang et al., 2015;Lei et al., 2009;Xia et al., 2016). However, these studies could not resolve the structure of the Hainan plume in the lower mantle. Using multiscale global tomography, Zhao et al. (2021) show that low-V anomalies exist in the whole mantle beneath SE Asia and these low-V anomalies are not derived from a single plume beneath Hainan. They suggested that a cluster of plumes may exist in the mantle beneath SE Asia and the Hainan plume may be the strongest one. Being determined by a joint inversion of abundant local and teleseismic data, our tomography model also reveals the 3-D velocity structure of the lower mantle and shows a significant low-V body beneath the stagnant slab in the MTZ, suggesting that the Hainan plume is ascending from the lower mantle. Beneath the Hainan Island and adjacent areas, the Vp perturbation in the MTZ is close to zero, which may represent weak areas or holes of the stagnant slab (Figures 7 and 8).
Pieces of geochemical evidence suggest that the Hainan plume-related magmatism was initiated at ∼22 Ma and has continued to the present-day, spreading over a wide area including Hainan Island, the Leiqiong Peninsula, the SCS, and Indochina (Xu et al., 2012;. It is considered that the lower-mantle-rooted plume (or plume cluster) was generated by deep subduction of slabs surrounding the SCS Zhao et al., 2021). The plume-ridge interaction may have occurred during the last spreading stage of the SCS. Geodynamic modeling suggested that the low-V zone (Hainan plume) in the upper mantle beneath the northern SCS may be interpreted as a result of subduction-induced mantle-return flows (Lin et al., 2019). Our tomographic model shows that a low-V anomaly appears beneath the MTZ and extends down to the bottom of our model, i.e., 1,600-km depth (Figures 7-9), which represents the lower-mantle-rooted Hainan plume. The low-V zone has a head spreading laterally beneath the MTZ. At depths of 100-430 km, an obvious low-V anomaly exists under the Leiqiong and Hainan areas, which may have caused the widely distributed Cenozoic basalts.
Most of the basalts in the SCS were generated in the Miocene, but the basalts in Hainan and Vietnam were mainly generated in the Quaternary (Xu et al., 2012). Such a temporal and spatial distribution of volcanism suggests that the Hainan plume was continuous from the lower mantle to the upper mantle in the early Miocene, which generated the relatively older basalts in the SCS. Then the subducted slab became stagnant in the MTZ as shown by our model (Figures 7-9), and the upwelling of the Hainan plume from the lower mantle is largely resisted by the stagnant slab. We think that the low-V anomalies at shallow depths (100-430 km) beneath Leiqiong and Hainan result from the hot upwelling from the lower mantle through the weak areas or holes of the stagnant slab within the MTZ. The widely distributed Cenozoic basalts in and around the SCS were related to the upwelling of the lower-mantle-rooted Hainan plume.

Seismic Anisotropy and Mantle Dynamics
Mantle flow induced by slab subduction could be detected by seismic anisotropy (Long, 2013;Zhao et al., 2016). Previous SWS studies made in SE Asia have provided important information on the subduction dynamics (Collings et al., 2013;Hammond et al., 2010;Lynner & Long, 2014;Song et al., 2021;Yu et al., 2018). The SKS splitting results show APM-parallel FVDs under the Sumatra Island (Collings et al., 2013), which may reflect the entrained flow below the subducting slab. The APM-parallel FVDs in Sumatra change to trench-parallel in Java (Hammond et al., 2010;Lynner & Long, 2014). However, these SKS studies cannot determine the source depth of anisotropy due to the limitation of the SWS method. Huang et al. (2015) revealed APM-parallel FVDs at shallow depths (<50 km) but trench-parallel FVDs at greater depths (50-600 km). However, their model has a limited resolution because they did not use teleseismic data and adopted strong smoothing and damping constraints in their tomographic inversion.
Our results show significant low-V anomalies in the upper mantle wedge above the subducting slab (Figures 8  and 9), which have been detected in most subduction zones all over the world, e.g., the Japan Islands (Liu & Zhao, 2016a, 2016bNiu et al., 2016;Zhao et al., 1994), Kamchatka (Jiang et al., 2009), Alaska (Gou et al., 2019;Tian & Zhao, 2012), Cascadia (Chen et al., 2015;, and the Philippines (Fan & Zhao, 2019;Zhao et al., 2021). These low-V anomalies are usually interpreted as subduction-driven hot and wet upwelling flows in the mantle wedge associated with the slab dehydration, which cause arc magmatism and back-arc spreading (e.g., Iwamori & Zhao, 2000;Liu & Zhao, 2016a;Zhao et al., 1994). Although the obtained Vp azimuthal anisotropy is complex in the low-V mantle wedge, there is an obvious feature that the FVDs within these low-V zones are trench-normal at shallow depths (<100 km) and trench-parallel in deeper areas ( Figure 10). The trench-normal FVDs in the low-V mantle wedge at shallow depths (<100 km) may represent the mantle wedge corner flow and its anisotropy generated by the slab subduction. This feature may indicate that the axis of A-type olivine is parallel to the direction of mantle flow generated by the slab subduction, which may also explain the results of seismic anisotropy in other subduction zones (e.g., Eberhart-Phillips & Reyners, 2009;Fan & Zhao, 2019;Gou et al., 2019;Huang et al., 2011;Liu & Zhao, 2016a;Wang & Zhao, 2008).
We also compare our 3-D Vp anisotropy model with the SKS splitting measurements in the study region. Our 3-D anisotropy model ( Figure 10) shows that the FVDs in the Sumatra-Java fore-arc are generally parallel with the APM direction at shallow depths (≤200 km), while the FVDs change to trench-parallel at greater depths (>430 km). This feature is different from the SKS splitting results, which show trench-normal FVDs beneath Sumatra and trench-parallel FVDs beneath Java (Hammond et al., 2010;Lynner & Long, 2014). The previous researchers suggested that this transition was caused by the different slab ages beneath Sumatra and Java. The younger lithosphere beneath Sumatra leads to strong coupling between the subducting slab and surrounding mantle, whereas the older lithosphere beneath Java may lead to weak coupling or decoupling between them, generating different mantle flow patterns and anisotropies (Lynner & Long, 2014). However, our 3-D Vp anisotropic model shows that most of the APM-parallel FVDs exist within the high-V subducting slabs, which cannot be interpreted by the subduction-induced mantle flow.
The APM-parallel FVDs may reflect fossil anisotropy in the subducting lithosphere, which was acquired at midocean ridges from the alignment of olivine crystals (crystallographic a axis) in the spreading direction (Blackman et al., 2002;Christensen, 1984;Hess, 1964;Shinohara et al., 2008;Tono et al., 2009). When the oceanic lithosphere is formed, it creates a large anisotropic signature and the anisotropy can be changed and reproduced during plate motion and subduction (Hammond et al., 2010;Wang & Zhao, 2012). The subducting slabs beneath Java and Sumatra may have reproduced their anisotropies during the plate motion, generating the APM-parallel FVDs.
SKS splitting measurements can constrain the mantle anisotropy with high lateral resolution, but their depth resolution is very poor. Our Vp anisotropic model has high vertical resolution, so it provides critical information on the depth variation of anisotropy in the upper mantle. Following Wei et al. (2019), we estimate synthetic SKS splitting parameters in SE Asia from our 3-D Vp anisotropic model ( Figure 11). We assume that 3-D anisotropic features of S waves are consistent with those of P waves. For weak anisotropy with a horizontal symmetry axis, the MATLAB seismic anisotropy toolbox (Walker & Wookey, 2012) and the method of Silver and Savage (1994) are used to synthesize SKS splitting parameters. The result (Figure 11) shows that the FVDs of the synthetic and observed SKS results are roughly consistent along the Sumatra and Java trenches, except for the Indochina area.
Another interesting feature of the 3-D Vp azimuthal anisotropy is that divergent shape FVDs occur in the significant low-V zones beneath the MTZ (750-km depth, Figure 10). A lower-mantle-rooted plume is revealed beneath the Hainan-Leiqiong-Indochina area (Figures 7-9). The upwelling of the hot plume may be resisted by the slab stagnating in the MTZ, resulting in the divergent shape FVDs beneath the stagnant slab. A previous study revealed a broad low-V zone beneath the MTZ under the Hawaiian Islands, suggesting that the Hawaiian plume materials were held back beneath the MTZ before ascending to the upper mantle (Cao et al., 2011). Our present results show that the structure of the Hainan plume is more complex because of the widely distributed deep slabs around the plume, e.g., the slabs beneath Sumatra, Java, and the Philippines. The lower-mantle-rooted Hainan plume may be divided into several smaller thermal upwellings in the upper mantle because of the existence of the stagnant slab in the MTZ, generating the widely distributed Cenozoic basalts in and around the SCS Xia et al., 2016;Zhao et al., 2021). The coexistence and interactions of the mantle plume and the stagnant slab in the MTZ may have caused the complex tectonic history and activities in the SE Asian region, including the many arc and back-arc volcanoes, active intraplate volcanism, and the widely distributed Cenozoic basalts (Figure 1).

Conclusions
A high-resolution 3-D tomographic model of P wave isotropic velocity and azimuthal anisotropy beneath SE Asia is determined by jointly inverting abundant arrival times of local earthquakes and relative travel-time residuals of teleseismic events, which provides important new information on the 3-D mantle structure and subduction dynamics beneath the SE Asia subduction system. The main findings of this study are summarized as follows: 1. Beneath the Sumatra Island, the subducting aseismic slab has penetrated through the mantle transition zone (MTZ) and slab rollback may take place in the lower mantle. In contrast, the high-V slab beneath the Java trench has a relatively lower dip angle and the slab is stagnant in the MTZ. The different slab geometries beneath Sumatra and Java may be associated with the lithosphere age of the subducting slab 2. A low-V anomaly is revealed beneath the MTZ and extends down to the bottom of our model (1,600-km depth), which represents the lower-mantle-rooted Hainan plume. The plume ascends to the upper mantle through weak zones or holes of the stagnant slab in the MTZ, generating the widely distributed Cenozoic basalts in the Leiqiong, Indochina, and SCS areas 3. APM-parallel FVDs are observed within the subducting high-V slab beneath Sumatra at shallow depths (<200 km), which reflect fossil anisotropy produced at the mid-ocean ridge and modified during the plate motion and subduction Figure 11. A comparison of the synthetic (red bars) and observed (blue bars) SKS wave splitting parameters. The orientation and length of the bars denote the fast polarization direction and delay time, respectively. The observed SKS splitting parameters are taken from http://splitting.gm.univ-montp2.fr/ (Barruol et al., 2009). The fast directions of the synthetic and observed SKS results are generally consistent along the Sumatra and Java trenches.
4. The upwelling of the Hainan plume is resisted by the stagnant slab in the MTZ, generating divergent FVDs that are centered in the significant low-V zones below the MTZ. The coexistence and interactions of the mantle plume and the stagnant slab may have caused the complex tectonic history and activities in SE Asia, including many arc and back-arc volcanoes, active intraplate volcanism, and the widely distributed Cenozoic basalts

Data Availability Statement
We used P wave arrival-time data that were downloaded from the International Seismological Center database (http://www.isc.ac.uk/). The 3-D tomographic models of isotropic Vp and azimuthal anisotropy obtained by this study are available at the website , Zenodo, https://zenodo.org/record/5115260). All the figures were plotted by using the free software GMT (Wessel et al., 2013).