In order to exclude the possibility of S-wave shadow caused by the source radiation pattern effect, the focal mechanism of the first earthquake was determined by the 1st motion P-waves recorded at seismic stations in the TVG as well as the Taiwan area (Fig. S4). Although the focal mechanism was not well constrained due to poor station coverage, the projection of the 1st motions of the P-waves recorded at seismic stations in the TVG were limited within a small zone of the dilatational quarter in the focal mechanism. Also the 1st motion of P-waves observed at those seismic stations were not far from one of the possible fault planes. Thus, the S-wave shadow observed at a few of neighboring seismic stations in the TVG couldn’t be caused by the source radiation pattern.

S-wave shadows (attenuation), like the early discovery of the liquid outer core of the earth22, provide evidence for the existence of a partially molten magma reservoir beneath northern Taiwan. The S-wave shadows recorded at stations YC03-04 and YL05-07 within an approximate area of 3 km × 8 km indicate molten magma exists somewhere along the ray-paths generated from the first earthquake (Fig. 6a). The S-wave shadows are hardly caused by any strong scattering from the Moho-discontinuity or complex structures in that the S-wave shadows are different from the coda waves, which often show seismic amplitudes decrease with time. The possibility of molten magma in the upper crust might be ruled out from the fact that no S-wave shadow was recorded from any seismic stations during the second earthquake. On the other hand, the velocity anomalies might not be located in the upper mantle because an extremely low velocity is required for causing P-wave delay of ~0.4 s within a small size of velocity anomalies limited by the area of the seismic stations that have P-wave delays. The extremely low velocity might not be reasonable in the partial melting of the upper mantle. Besides, P-wave delay of ~0.4 s recorded at seismic stations within a short distance (<10 km) in the Tatun volcano group of Taiwan is difficultly caused by the complex Moho-depth variation unless the seismic stations are just across the major plate boundary. Although there is currently not enough information to identify the exact depth of the magma reservoir, a reasonable estimation of its depth might be at the lower crust, such as the deeper reservoir beneath Yellowstone National Park in the USA23. Thus, the projection of the S-wave shadows down to the Moho depth, based on the calculated ray-paths, suggests the possible location of magma reservoir shifted approximately 10 km eastward at the lower crust, roughly beneath the Chinshan area (Fig. 6a).

Figure 6 (a) Location projection of seismic stations (triangles), S-wave shadows (broken lines in red), and P-wave delays (dotted lines) on the surface down to the Moho depth based on the ray-paths from the deep earthquakes beneath northern Taiwan. Schematic plots for a magma reservoir composed of either (b) a thin molten overlay (red) or (c) many molten sills (red) within thick partially molten rocks (pink). All figures were plotted with GMT v4.5.2 (gmt.soest.hawaii.edu). Full size image

In addition to the S-wave shadows, the existence of a magma reservoir can be strongly confirmed by the P-wave velocity delay at particular seismic stations at the TVG. In fact, the P-wave delay was detected not only at the seismic stations where S-waves attenuated, but also at surrounding stations in the broader area of 15 km × 6 km (Fig. 6a). It is well known that molten magma, as well as partially molten rock, will slow down the propagation speed of P-waves. Also, the delay from seismic wave propagating through the molten magma is significantly larger than that from wave propagation through partially molten rock. Thus, the lack of a significant difference in P-wave delay between seismic stations with and without S-wave shadows may simply suggest that a thin molten layer may be overlying thick partially molten rocks at the lower crust beneath the TVG of northern Taiwan (Fig. 6b). In other words, the travel-time delay through the thin molten layer was significantly smaller than that caused by the thick partially molten rocks. A similar magma reservoir model of a molten magma layer overlying a partially molten low velocity zone was observed at the East Pacific Rise24. Although the tectonic setting between the East Pacific Rise and the TVG are different, the generation of magma reservoir might be similar because both of them are under the extensional regime17.

Alternatively, the presence of many molten sills as inferred in the lower crust beneath the Toba caldera25 might be another plausible model explaining the lack of a significant difference in P-wave delay between seismic stations with and without S-wave shadows (Fig. 6c). Again, comparing with thick partially molten rocks, the travel-time delay caused by sills might not be significant. Based on the experimental data and numerical modeling, Annen et al.26 proposed a model of “deep crustal hot zones” for the generation of intermediate and silicic igneous rocks. In the model, the mantle derived basalt emplaced into the lower crust as a succession of sills. Then partial crystallization of basalt sills not only can generate residual H 2 O-rich melts, but also provide heat and H 2 O for partial melting of pre-existing rocks. For examples, at the higher pressure (1.2 GPa), melt fraction would be ~0.5% if one considered H 2 O is 2.5 w% at temperature of 1100 °C (Fig. 5 in Annen et al. 2006). Under the same condition, the melt fraction would increase to ~0.56% at the lower pressure (1.0 GPa). However, there is not enough information to distinguish which model might be the most acceptable at this moment. Regardless of the exact model, both S-wave shadows and P-wave delays suggest that a magma reservoir may exist beneath the Taipei metropolis due to the ascending melt from the asthenospheric mantle17.

The delay of approximately 0.4 sec in P-wave arrivals provides further information for the estimation of the possible thickness of a low-velocity zone either in the lower crust or around the Moho depth. Assuming a P-wave velocity of 6.5 km/s in the lower crust, the thickness of the low velocity zone would be about 3.9 km, 6.1 km, and 10.4 km for 60%, 70%, and 80% of the velocity of the lower crust, respectively. Based on the model calculation of seismic velocity from petrofabrics and average shape of the melt phase24, the low-velocity zones of 60%, 70%, and 80% velocity of the lower crust are roughly associated with melt fractions of approximately 34%, 23%, and 14% in gabbro, which almost has the same geological composition as basalt in the lower crust. In other words, the thickness of the low-velocity zone caused by 14% partially molten rocks is approximately 10.4 km, but the thickness of the low-velocity caused by 34% partially molten rocks will be only approximately 3.9 km. With a larger fraction of partial melting, the layer will be thinner. Considering 34% partial melting within the magma reservoir, its size could be approximately 350 km3 (3.9 km × 15 km × 6 km). If 14% partial melting was estimated, the reservoir’s size might be up to 936 km3 (10.4 km × 15 km × 6 km).

In summary, a deep magma reservoir has been identified beneath the Taipei metropolis of Taiwan based on both S-wave shadows and P-wave delays. The magma reservoir is probably composed of either a thin magma layer overlaying a thick partially molten zone or many molten sills within the partially molten rocks. The size estimate of the magma reservoir is dependent on the fraction of partial melting within it. Considering 34% partial melting within the magma reservoir, its size could be approximately 350 km3. Alternatively, the size of magma reservoir may be up to 936 km3 if the partial melting is only 14%. Although the exact depth is not yet well constrained, a detailed image of the magma reservoir will be obtained in the next four years (Fig. S4). A new project is planned to deploy 120 broadband seismic stations with a spacing of ~5 km to create a grid covering northern Taiwan, 45 ocean-bottom-seismometers with a spacing of ~10 km offshore from northern Taiwan, and 500 geophones (3-components) with a spacing of ~500 m at the TVG in 2017–2020.