Unveiling Tatun volcanic plumbing structure induced by post-collisional extension of Taiwan mountain belt

The Tatun Volcanic Group (TVG) is proximal to the metropolis of Taipei City (population of ca. 7 million) and has long been a major concern due to the potential risks from volcanic activity to the population and critical infrastructure. While the TVG has been previously considered a dormant or extinct volcano, recent evidence suggests a much younger age of the last eruption event (~ 6000 years) and possible existence of a magma reservoir beneath the TVG. However, the location, dimension, and detailed geometry of the magma reservoir and plumbing system remains largely unknown. To examine the TVG volcanic plumbing structure in detail, the local P-wave travel time data and the teleseismic waveform data from a new island-wide Formosa Array Project are combined for a 3D tomographic joint inversion. The new model reveals a magma reservoir with a notable P-wave velocity reduction of 19% (ca. ~ 19% melt fraction) at 8–20 km beneath eastern TVG and with possible northward extension to a shallower depth near where active submarine volcanoes that have been detected. Enhanced tomographic images also reveal sporadic magmatic intrusion/underplating in the lower crust of Husehshan Range and northern Taiwan. These findings suggest an active volcanic plumbing system induced by post-collisional extension associated with the collapse of the orogen.

To illuminate the volcanic plumbing structure of the TVG, a deployment of dense broadband seismic array, named Formosa Array 23 , was launched in 2017 and has installed 120 out of 140 planned stations by October 2019 (Fig. 1b, blue squares). The station spacing is about 5 km uniformly across plain and mountain areas, providing a new opportunity to improve tomographic imaging of northern Taiwan [24][25][26][27] . In this study, we combine the teleseismic waveform data from the Formosa Array and the local P-wave picking data from an integrated seismic network to conduct a 3D joint inversion (Fig. 1c,d). The new model unveils a clear magma reservoir beneath the east of the TVG and possible magmatic intrusion/underplating bodies related to post-collisional extension in northern Taiwan.

Results
Model results and verification. The model slices at different depths are shown in Fig. 2, where poorly resolved areas are masked and defined as a resolvability index R < 0.6 calculated from checkerboard tests (Supplementary information and Supplementary Fig. S5). Benefiting from the uniform and dense station distribution of the Formosa Array, the joint inversion integrating teleseismic data greatly improves the resolution to the north of latitude 25.3° and at depths of 5-50 km when compared to the inversion results without teleseismic data (Supplementary Figs. S1 and S6). With the enhanced images, a pronounced low velocity anomaly (L1) stands out beneath the Tatun volcano group (TVG) at the depth slice of 8 and 16 km (Fig. 2b,c). This slow anomaly is elongated sub-vertically at a depth range of 8-20 km as shown in cross-sections AA′ and BB′ (Fig. 3b,d), with most of the seismicity lying above the anomaly. According to the projected locations of craters (green triangles in Fig. 3), this slow anomaly is not located at the center of TVG (i.e. Chihsingshan, CHS) but slightly to the east beneath the Huangzuishan (HGS) and Dayoukeng fumarole (DYK). In the cross-section AA′, it can be seen that the reservoir extends upward to a shallower depth of ~ 5 km and northward to the offshore areas near where active submarine volcanoes (SV and Ksv in Fig. 1b) have been documented 9 . Deep velocity structures beneath northern Taiwan are shown in regional cross-sections in Fig. 4. In crosssection CC′, the subducting PSP is characterized by the seismicity with high velocity anomaly, H1 (Fig. 4b,c). The low velocity anomaly L2 at a deep depth around 70-80 km has been identified and suggested as the subductioninduced partial melting in the previous studies 23, 25 . This slow anomaly L2 seems not to be connected to the Tatun slow anomaly (L1) as a feeding source. Other two pronounced low velocity anomalies, L3 and L4, at shallower depths have been observed previously and inferred as the serpentinized mantle wedge 25,27 and sediments of Ilan plain 25 , respectively. Cross-section DD′ provides another angle of view for these structures (Fig. 4d,e). Further to the west in cross-section EE' , it is noteworthy that an abrupt transition in crustal structure is observed from www.nature.com/scientificreports/ thickened crustal root in the south to a few high velocity anomalies, denoted as H2, in the north (Fig. 4f,g). This transition occurs roughly along the boundary between the Central range and Hsuehshan range (Fig. 2d). Characteristic-model tests are conducted to evaluate the robustness of the geometry of the TVG L1 slow anomaly ( Supplementary Fig. S7). We test vertically-elongated and concentrated synthetic slow anomalies at depth ranges of 8-20 and 12-16 km and find that the vertically-elongated shape of L1 is resolvable and not a consequence of smearing effect. Other tests including data noise and different initial models are also conducted: adding random noise with 0.12-s standard deviation (based on the residual RMS of actual inversion) to the synthetic data and using a different initial model (e.g. 1-D velocity model) do not change the tomographic imaging that much (Supplementary Figs. S8 and S9). Detailed descriptions of the tests are referred to the Supplementary Information.

Discussion
Magma reservoir of Tatun volcano. The new model derived from joint inversion of local earthquake and teleseismic data unveils a pronounced low velocity anomaly (L1) beneath the TVG, interpreted as the magma reservoir (Figs. 2, 3). It is vertically elongated in shape and is about 10 km in diameter at depths of 8-20 km, which is shallower than the previously proposed depth 18,19 . This vertically extended geometry of the imaged reservoir seems to be robust through characteristic-model tests ( Supplementary Fig. S7-S9). The maximum P-wave velocity (Vp) reduction of the reservoir is about 12%. Since the regularization of tomographic inversions often damps the magnitude of velocity perturbations, we conduct synthetic tests with different input Vp reductions for the reservoir to assess the actual velocity anomaly magnitude (Supplementary Information and Supplementary  Fig. S10). The results find that a greater Vp reduction of 19% is needed to fit the inversion results (Supplementary The high value of 19% is difficult to explain by temperature and composition alone and implies the presence of melt. We estimate the melt fraction of the magma reservoir by using Gassmann's relations 30,31 . The calculation and used parameters are described in "Data and methods" and Supplementary Table S2 in detail. Assuming andesitic melt and granite frame rock, 19% Vp reduction yields a high melt fraction of ~ 19% (Fig. 5). Alternatively, replacing andesitic melt by water or CO 2 for a hydrothermal reservoir results in a pore fraction of ~ 15% (Supplementary  Table S2). While it is unlikely to have hydrothermal systems down to as deep as 20 km due to high pressure, a certain mixture such like partial melt intruded at depths and hydrothermal gas/fluid accumulating at the top is possible. Although the uncertainty of this simplified estimation is not trivial, the high melt fraction indicates that the TVG plumbing system could remain active with sufficient heat supply.
Pu et al. 32 recently relocated the micro-earthquakes in the TVG area and revealed a clear conduit-like structure beneath the Dayoukeng fumarole, which requires volcanic gas and fluid source ascending from below. The consistent location of the magma reservoir we resolved can potentially be the gas/fluid source to either replenish a shallower hydrothermal reservoir or feed directly to the fumaroles through fractured pathways (Fig. 3b). www.nature.com/scientificreports/ Origin of Tatun volcano magma reservoir. In a regional volcanism context, the next key piece of information is to know the origin of magma in the reservoir. The TVG is located above the edge of the northward subducting Philippine Sea Plate (Fig. 1a). But from Fig. 4, it is clear that the magma reservoir (L1) is not fed by current subduction-induced partial melting (L2) as no pathway is identifiable in between. The magmatism of the TVG, as well as the entire northern Taiwan volcanic zone (NTVZ, Fig. 1b), has been suggested to be part of the Ryukyu volcanic front that was initiated by the westward advance of Philippine Sea Plate before 2 Ma and ceased by the following slab rollback and opening of the Okinawa Trough 3,4 . Wang et al., in contrast, suggested a mechanism of post-collisional delamination to form the NTVZ since geochemical characteristics of NTVZ magmas show significant components of asthenosphere and metasomatized subcontinental mantle 4,6 . If extensional magmatism occurred, we envisage that the crust should have been subjected to extensive magmatic intrusion or underplating which is often characterized by high velocity bodies in the lower crust or beneath Moho 33 . In Figs. 2d and 4f, we do observe such plausible high velocity bodies (denoted as H2) sporadically present in the lower crust of the Hsuehshan range and its north according to the Moho depth around 30 km in this region [34][35][36] .
To the south, the crust is thickened and relatively intact along the Central Range (Fig. 4f,g). Distribution of these faster anomalies also largely overlays with the extent of NTVZ volcanism including the Tsaolingshan lava site (TLS) at the southern end of NTVZ 7 (Fig. 2d, purple diamond), as strong evidence of the magmatic intrusion/ underplating during the post-collisional extension of northern Taiwan mountain belt which is still ongoing to the present day 5,37 . A clearer 3-D perspective view of northern Taiwan velocity structure is shown in Fig. 6. In addition to features, such as subduction-induced partial melting, serpentinized mantle wedge, and thickened crustal root, the new model unveils the location and geometry of the Tatun magma reservoir and magmatic intrusion/underplating bodies in the lower crust of northern Taiwan, which sheds light on a clearer picture of post-collisional magmatism and tectonics of Taiwan orogeny. In fact, during the slab roll back of Teng's model 3,4 , extensional magmatism could also occur 38 . By current tomographic imaging solely, it could be difficult to differentiate the slab roll back 3,4 and the delamination model 4,6 . A higher-resolution imaging of slab interaction between the PSP and EP at deep depths beneath the northern Taiwan may help resolve the debate. However, limited by the relatively narrow aperture of Formosa array, the teleseismic data do not provide additional resolution deeper than 60-70 km (Supplementary Fig. S5).

Concluding remarks
The 120 uniformly-distributed broadband stations of the Formosa Array provide a unique opportunity to illuminate the Tatun volcanic plumbing system and deep subsurface velocity structure beneath northern Taiwan (Fig. 1b). The wealth of this new dataset will facilitate further studies and findings to come. The unveiled Tatun magma reservoir shows a shallower location (8-20 km) than previously thought and high melt content (~ 19%). Since no clear mantle feeding source is identified, whether the melt is residual from a recent eruption or is new from the mantle source is unknown. While the high melt does not imply an imminent eruption, it could probably sustain an active hydrothermal system of the Tatun volcano in the long run. Monitoring all forms of activity in the TVG (e.g. hydrothermal, microseismicity, gass/fluid chemistry, CO 2 /SiO 2 flux, ground deformation, etc.) is therefore of great importance for volcanic hazard forecasting and mitigation.

Data and methods
Seismic networks and data processing. We use both local and teleseismic earthquake data for a joint inversion in this study. The local earthquake data are compiled from an integrated regional seismic network that combines (1) (5) stations from the TAIGER project that include inland broadband stations and western offshore OBS arrays 43 . A total of 671 stations are used in the study area (Fig. 1b, purple squares). Based on the data criteria that (1) each event was recorded by at least 4 stations, (2) gap angle of recorded stations less than 180°, (3) focal depth shallower than 120 km, and (4) picking quality less than level 2 (according to the CWB picking criteria), 60,081 events are sorted out from 1991 to 2017. To avoid the dominance from compactly-distributed earthquake clusters for better conditioning the inversion, a 3-D event grouping method 25 with 2 km radius is applied to homogenize the event distribution and reduce the event number to 3587 with 37,373 readings in the end (Fig. 1d).
For the teleseismic data, we collect waveform data from 371 teleseismic events (M L ≥ 5.0) in a distance range of 30°-90° recorded by the newly-deployed Formosa Array (Academia Sinica, Institute of Earth Sciences, 2017) from April 2018 to October 2019 (Fig. 1c, gray dots). The vertical component of waveform data is bandpassfiltered in 0.08-0.15 Hz to measure the P-wave relative travel times between stations using the adaptive stacking method 43 . Through the iterative measuring process, only those with signal-to-noise ratio larger than 5 and crosscorrelation coefficients (CC) greater than 0.9 are retained. After final visual inspection and avoiding dominance in the inversion by the events in Tonga subduction zone from the southeast, we sort and even the events in four quadrants and obtain a total of 70 quality events with 4891 readings (Fig. 1c, red dots). Supplementary Fig. S2 shows the measured relative P-wave travel times for four example events coming from different quadrants. The travel times for the stations around the TVG are consistently delayed regardless of the event azimuth, implying the existence of low velocity anomalies underneath. 3-D distribution of ray paths from the local and teleseismic events are shown in Supplementary Fig. S1.

Tomographic joint inversion and resolution tests.
A recently developed tomographic code for multidataset joint inversion is employed in this study 25,26,31 , in which the absolute travel-time residuals from local earthquake data and the relative travel-time residuals from teleseismic data are simultaneously minimized (Supplementary Information). The model grids are parameterized with 4 km in longitude and latitude and increasing intervals from − 5 to 120 km in depth as shown in Supplementary Fig. S1a and Supplementary Table S1. The 3-D www.nature.com/scientificreports/ P-wave velocity (V P ) model of Huang et al. 25 is interpolated to the current model grids and used as the initial model for joint inversion. The damping and smoothing factors are chosen to be 20 and 10 after a series of tradeoff tests (Supplementary Fig. S3). We iterate the inversion until the root-mean-square of travel-time residuals (RMS) reduces insignificantly. After 3 iterations, the RMS is reduced from 0.29 to 0.18 s (~ 40% reduction). The relatively less RMS reduction is because of the 3-D initial model used. The residual distribution for local earthquake and teleseismic data is plotted separately in Supplementary Fig. S4. Both the residuals are effectively reduced and concentrated toward zero after the inversion.
The checkerboard test is conducted for assessing general resolution power of the obtained velocity model (Supplementary Information). In Supplementary Fig. S5, the results of the checkerboard tests show good recovery in most inland region and to the east of longitude 122.1° E as well as above 30 km depth. Below 30 km, the well-recovered region is gradually shifted eastward as the depth increases. At 90 km depth, the well-recovered region is primarily limited beneath the Ilan plan and its east offshore. The recovered checkerboard model is then used to calculate a resolvability index, R, at each model node (Supplementary Equation S5 and Supplementary  Fig. S5). R ranges from 0 to 1 where the value 1, 0.5, and 0 indicate the node is 100%, 0%, and − 100% recovered, respectively. R = 0.6 is used as a lower bound for the resolvable nodes in tomographic images 31 .
Estimation on melt fraction and volume of magma reservoir. We calculate effective elastic moduli of fluid-saturated porous materials using Gassmann's relations 44 as where K o , K d and K f represent bulk modulus for the material of the frame rock, the drained matrix, and the saturating fluid, respectively. µ d is shear modulus of the drained matrix. ϕ is porosity. For porosity below the critical value ( ϕ c ) that defines the transition from a medium which is supported by the frame to a medium where solid materials are suspended in the fluid 45 Table S2 collects the physical parameters of rocks and fluid for calculation, where we assume the frame rock of granite and saturating fluid of andesitic melt for imaged magma reservoir (L1, Fig. 3). The calculation then yields a relationship of the V p reduction and porosity (i.e. melt fraction) as shown in Fig. 5. So, given the V p reduction of 19% for the magma reservoir from characteristic-model tests in Supplementary Fig. S10, the melt fraction is estimated to be ~ 19%. (1)