1 Introduction

The Tatun volcano group (TVG) is composed of more than 20 volcanoes (Fig. 1a), which mainly erupted 0.4–0.8 Ma ago (Chen and Wu 1971; Juang and Chen 1989; Taso 1994; Teng 2007). The volcanic features, such as hot springs and fumaroles, are still significant on the surface (i.e., Chio and Lin 2017), and thus the major part of the TVG has been preserved as the Yangminshan (YMS) national park since 1985. The TVG was often considered to be an extinct volcano group in that there were no recorded volcanic eruptions in human history. However, a variety of recent studies from geochemical, seismic, geodetic, and geological data indicate that the TVG is still active. First, geochemical analyses for Helium isotope from samples collected at hot springs and volcanic gases show a significant mantle signature (Yang et al. 1999). Secondly, many shallow micro-earthquakes have clustered at the TVG (Lin et al. 2005a, b). Among those clustering earthquakes, some typical volcanic earthquakes such as tornillos, spasmodic events, and very-long-period events have been detected (Lin et al. 2005a, b; Konstantinou et al. 2007; Lin and Pu 2016). Recently, it has been interesting to see particular volcano-seismic events such as phreatic eruptions and heartbeat-like seismicity (Lin 2017a, b). Thirdly, some pressure sources in the uppermost crust have been observed by the high accuracy leveling survey of the TVG (Murase et al. 2014), which is a technique of calculating the height of one level on the surface relative to another. Fourthly, the results from dating volcanic ashes and magnetite shows the last eruption might be younger than 10,000 years (Belousov et al. 2010; Zellmer et al. 2014). Such recent volcanic activity fits well with the empirical criteria for defining an active volcano (Szakacs 1994). Finally, a deep magma reservoir beneath the TVG has been identified from both S-wave shadow and P-wave delay (Lin 2016). This result agrees with the phenomenological definition of an active volcano. Thus, The TVG is unambiguously considered an active volcano according to both empirical and phenomenological criteria.

Fig. 1
figure 1

Locations of the dense geophone array, background geology and tectonics in and around northern Taiwan. a Plots of a dense geophone array (triangle), Shanchiao fault (dashed line), nuclear power plants (hexagon) and Taipei city (square) on a topographic map in and around the Tatun volcano group. b Distribution of the array and three felt earthquakes (circles) in the northern Taiwan area (see Table 2). The insert map on the lower-right panel shows the general tectonics in and around the Taiwan area

Since the TVG is active, some interesting issues have recently been identified from different considerations of volcanic hazard and geothermal exploration. It is obvious that the potential for a volcanic eruption might not be totally excluded in the future. Since the TVG is close to the Taipei metropolis, with more than 7 million residents living in two large cities (Taipei City and New Taipei City), possible volcanic disasters such as lava flow, pyroclastic flow, and volcanic ash would have a significant social and economic impact. In fact, the distance between Mt. Chihsing (the youngest and highest volcano among the TVG) and the Taipei 101 building in downtown Taipei is less than 15 km. Additionally, two nuclear power plants are located along the northern coast of Taiwan, just north of the TVG (Fig. 1a). As a result, the potential volcanic threat cannot be ignored in the future. On the other hand, the TVG may provide some geothermal resources for replacing energy produced by traditional power plants from nuclear or fossil sources. For the purpose of both volcanic impact and geothermal exploration, therefore, it is very important to identify more detailed subsurface structures beneath the TVG.

In this study, we deploy a dense linear geophone array across the TVG to measure subsurface structures with better resolution. Due to the eco-preservation within YMS national park, it is almost impossible to install many ordinary seismometers in the TVG. Therefore, a geophone array is used instead of seismometers. The array is comprised of 50 geophones with an average spacing of ~ 200 meters, covering a length of more than 10 km. The array is designed for striking along the NW–SE direction to record the abundant earthquakes in the Ilan area of Taiwan, which has both deep earthquakes within the subduction zone and also many shallow earthquakes along the back-arc opening in the Okinawa trough. The crosscheck seismic ray-paths of both shallow and deep earthquakes from the SE direction provide a good opportunity to examine the subsurface structures beneath the TVG.

2 Tectonic and Geological Background

The TVG is located at the northern tip of Taiwan (Fig. 1b), where it is just on the top of the western end of the subducted Philippine Sea plate beneath the Eurasian plate (Tsai et al. 1977; Lin 2016). Thus, the possible volcanic genesis of the TVG is often directly linked to the volcanic arc (Teng 1996). However, this idea has been strongly argued against by other observations. First, the depth of the subducted slab beneath the TVG is significantly deeper than the typical depth of ~ 100 km beneath the volcanic arc in general (Lin 2016). Second, the isotope analyses indicate that the TVG does not have an arc signature, but it is more like the rifting magma in the continental margin (Wang et al. 1999, 2004). This is another possible reason for the volcanism at the TVG, considering the volcanoes lie along a subduction-transform edge propagator fault (Govers and Wortel 2005; Hale et al. 2010; Nijholt and Govers 2015). This volcanism is produced by the lateral leaking of mantle fluids at the end of the western boundary of the subducted Philippine Sea plate. Regardless of the exact genesis of the TVG, recent observations show it is still active.

It is interesting that all active volcanisms in the TVG are likely associated with the Shanchiao fault, which is the only active fault in and around the Taipei metropolis (Lin et al. 2008). The distribution of volcanic features such as hot springs and fumaroles at the TVG is roughly parallel to the Shanchiao fault (Fig. 1a). The Shanchiao fault is a normal fault striking in the NE-SW direction, and it is re-active along the Chinshan fault, which was a thrust fault active in the Miocene (Chen et al. 2014). It is also worth noting that all of the active volcanisms, including hot-springs, fumaroles, and volcano-earthquakes, are limited at the hanging-wall part of the Shanchiao fault. In other words, there is no volcanic activity at the footwall part of the Shanchiao fault. Therefore, it seems that the Shanchiao fault might be the western boundary of active volcanism among the TVG.

3 Dense Linear Geophone Array

In total, 50 vertical component high-frequency geophones are deployed with an average spacing of ~ 200 m across the TVG along the NW–SE direction (Fig. 1a and Table 1). The geophone sensors are flat to velocity at frequencies above 4.5 Hz. The array is roughly perpendicular to the strike of the main geological structures (Ho 1988). This dense geophone array is mainly designed for recording the many earthquakes in the NE Taiwan area in which the Philippine Sea plate is subducting beneath the Eurasian plate. The seismic data are recorded with a RefTek (Texan) recorder. The same instruments had been installed for studying non-volcanic tremors in the southern Taiwan area (Sun et al. 2015). The deployed geophones and equipment are shown in Fig. 2. The signal conditioning is provided by the Data Acquisition System (DAS), which consists of preamplifiers, analog-to-digital converters, and on-site data storage recorded to a flash memory card. All geophones were successfully operated for more than 6 months, from August 1, 2015 to January 6, 2016. Seismic data were recorded at a sampling rate of 100 Hz.

Table 1 Station coordinates and measured travel-times for Events 78 and 84
Fig. 2
figure 2

Photos for showing a geophone and other equipment deployed in the field

Since the major part of the TVG is preserved as the YMS national park, deployment of dense geophone arrays has significant advantages for examining the subsurface structures beneath the volcanoes. First, it is very easy to deploy geophones in the YMS national park. Since geophones are small and light, it is extremely convenient to carry them to any place, even locations that are very difficult to reach by cars and other vehicles. Secondly, the small battery provides enough power for recording data during the experimental period without concern for power supply from AC or solar panels. Thirdly, geophones are simple to install without any complicated construction that would be typical in deploying seismic stations. Finally, the individual noises recorded at a particular site might be removed by the nearby geophones, even though they are often recording short-period signals with poor signal-to-noise ratios.

4 P-Wave Delays

We checked all felt earthquakes (M > 4) that occurred during the period following array deployment in 2015, and discovered two events E78 and E84 were approximately located along the SE extended direction of the linear array; these were ideal for examining the possible subsurface structures along the NW–SE profile (Fig. 1b). One deep earthquake is located within the subduction zone; the other shallow earthquake occurs at the end of the Okinawa trough nearby Taiwan (Table 2). Careful examination of these two earthquakes shows some interesting variations of direct P-waves recorded at the dense linear geophone array because some significant subsurface structures might be discovered (Table 1).

Table 2 Earthquake parameters provided by Central Weather Bureau in Taiwan

Initially, the seismograms generated by the deep earthquake at a depth of 95.6 km on Nov. 11, 2015 (Event 84 at Table 2) show that the arrival of most direct P-waves are generally lineated along the blue dashed-line plotted by an apparent velocity of ~ 10 km/s (Fig. 3; Table 1). It is expected to see the P-wave arrival times with such an apparent velocity because the earthquake is located at the mantle depth. Although there are some uncertainties (less than 0.05 s), it is very interesting to see that several P-waves delays (marked by the red dashed-line in Fig. 3) are recorded at stations between T027 and T041. For those delays, the dramatic increase of the apparent velocity from ~ 10 to ~ 40 km/s indicates that the incidence angles of P-waves near the surface probably shift to the nearly vertical direction. This could be the result of low-velocity structures due to the ray refraction.

Fig. 3
figure 3

P-wave arrivals generated by Event 84 and recorded at the linear seismic array in the Tatun volcano group of Taiwan. Most of the first arrivals are lineated by blue dashed lines with an apparent velocity (Vap) of ~ 10 km/s, except some significant delays marked by a red-dashed line

Similar results are also found from the seismic data generated by another shallow earthquake on Nov. 11, 2015 (Event 78 in Table 2). Some P-wave delays are again found between stations T027 and T041 (Table 1). For most P-waves, an apparent velocity of ~ 7.4 km/s is expected because the earthquake source is located around the mid-crust. For those P-wave delays, the increase of the apparent velocity from 7.4 to 12.6 km/s suggests a large change in incident angles of the P-waves (Fig. 4). Again, the change of incident angles might be associated with the ray refraction from the high-velocity zone into the low-velocity zone. This phenomenon will be shown later as the forward modeling is applied.

Fig. 4
figure 4

P-wave arrivals generated by Event 78 and recorded at the linear seismic array in the Tatun volcano group of Taiwan. Most of the first arrivals are lineated by blue dashed lines with an apparent velocity (Vap) of ~ 7.4 km/s, except some significant delays marked by a red-dashed line

Although all P-wave arrivals are not corrected by station altitudes, the similar pattern of the P-wave delays in Figs. 3 and 4 reflects some particular subsurface structures across the TVG. The consistent pattern of P-wave delays gradually increasing from T027 to T041 is not the result of altitude variation of geophone sites because there is no significant correlation (see Fig. 5). Also, the significant 0.3 s delay increase from T041 to T042 cannot be largely attributed to the ~ 250 m altitude difference between both stations. Given reasonable velocities ranging from 3.0 to 4.0 km/s, the altitude difference of 250 m only causes P-wave delays between 0.06 and 0.08 s (Fig. 5). Those calculated differences caused by the station altitudes are significantly less than the observation delay of ~ 0.3 s between the two stations (T041 and T042).

Fig. 5
figure 5

Elevations (circles) and travel-time delays (dots) due to elevation variations of 50 seismic stations (T001–T050) shown in Fig. 1. Among those stations, the red ones mark the geophones with some travel-time delays generated by Events 84 and 78. Dots are the estimated travel-time delay due to elevation variations given by the P-wave velocities of 3 and 4 km/s, respectively

In order to further evaluate the possible effect of altitude difference between sites T041 and T042, we check P-wave arrivals generated from another earthquake located in a SW direction to the array (Event 98 at Table 2 and Fig. 1b). As can be seen in Fig. 6, it is clear that no significant delay is recorded at geophones between stations T041 and T042. This observation confirms that those significant P-wave delays from Events 84 and 78 are probably the result of some low-velocity zones beneath the array and are not primarily caused by altitude variation across the geophone sites.

Fig. 6
figure 6

P-wave arrivals generated by Event 98 and recorded at the linear seismic array in the Tatun volcano group of Taiwan. The seismograms recorded at geophones at sites of T041 and T042 are plotted by thick lines. The dashed line shows the arrivals with an apparent velocity (Vap) of ~ 6.0 km/s

5 Forward Modeling

In order to improve the understanding of possible subsurface structures causing P-wave delays recorded at parts of the geophone array, we employ a 2-D ray-tracing algorithm (Luetgert 1992) to perform the forward modeling. A simplified 1-D velocity structure is assumed for calculating two groups of ray-paths and arrival times of the P-waves generated by two earthquakes (Events 84 and 78 in Table 2). The velocity model is simply divided into 3 layers with depths divided at 2 km and 35 km. The P-wave velocities gradually increase from 4.0 to 5.0 km/s at the top layer (0–2 km) and from 5.0 to 6.6 km/s at the 2nd layer (2–35 km). The constant P-wave velocity is 7.8 km/s below 35 km.

For fitting the P-wave delays recorded at sites of geophones T027–T041 generated from two earthquakes at different focal depths and epicentral distances, we assumed a low velocity zone of extremely shallow depths (around 0.5–2.5 km) exists. The geometry of the low-velocity zone is selected to match the observations of a gradual increase in P-wave delays from T027 to T041 from both earthquakes shown in Figs. 3 and 4. The maximum delay due to the low-velocity zone (LVZ) is approximately equal to the observation of ~ 0.3 s (Figs. 7, 8). To account for both the features of the gradual increase of the P-wave delays from T027 to T041 and the significant jump between T041 and T042, an asymmetrical body with averaged velocities ~ 60% lower than the surrounding rocks is obtained (Fig. 6). For Event 84, the general pattern of the predicted arrivals are basically fitting well with the observations (Fig. 7). For Event 78, the general pattern between the observations and predictions might be still acceptable even though the differences between them are getting larger (Fig. 8). Those differences might be largely resulting from the difficulty of arrivals picked on the seismograms with larger uncertainties (Fig. 4). In general, it is well recognized that the result obtained here is not a unique solution and it has some geometrical uncertainties, but such a preliminary result still provides not only a reasonable explanation for the significant P-wave delays recorded at sites between T027 and T041 but also is related with results from some of our previous studies (i.e., Lin et al. 2005a, b; Murase et al. 2014; Lin 2016; Huang et al. 2017).

Fig. 7
figure 7

Ray-tracing (the lower panel) of Event 84 for modeling the travel-time delays (the upper panel) due to a low-velocity zone (red region) beneath the dense linear seismic array

Fig. 8
figure 8

Ray-tracing (the lower panel) of Event 78 for modeling the travel-time delays (the upper panel) due to a low-velocity zone (red region) beneath the dense linear seismic array

6 Discussion

The location of the low-velocity zone (LVZ) obtained here with some geometrical uncertainties is generally consistent with previous results from inversion of seismic ambient noise (Huang et al. 2017). Although the tomographic images from the ambient noise show several low-velocity regions at extremely shallow depths (Layers 1 and 2) by ambient noise inversion, only one continuously extends into deeper regions according to longer period images (Fig. 9b). The location of this major low-velocity zone below 1 km is roughly consistent with the location of the low-velocity zone obtained here. In other words, the low-velocity zone is consistently imaged from both the surface waves by ambient noise inversion (as described in the previous study) and the body wave delays (in this study).

Fig. 9
figure 9

modified from Huang et al. 2017)

Comparison between the low-velocity zone obtained here and one of the previous studies. a Schematic model showing the possible connection among the fault zone, hydrothermal reservoir in the extremely shallow crust, and the suspected magma reservoir in the upper crust, and b depth distributions of low-velocity zones (the dotted circles) imaged from the ambient noise inversion (

The low-velocity zone obtained here is likely associated with the major hydrothermal reservoir on the hanging-wall side of the Shanchiao fault (Fig. 9a). Identification of such a low-velocity zone in the shallower crust is associated with several previous observations (Fig. 10). First, the low-velocity zone is located closely beneath high fumarolic activity areas, such as Dayoukeng, Hsaioyoukeng, and Huangshan (Fig. 1), in which the volcanic degassing process is often strong. Secondly, a clustering of micro-earthquakes has been detected in and around the LVZ (Lin et al. 2005a, b) (Fig. 10b). In addition to the volcano-tectonic earthquakes, we have occasionally detected typical volcanic earthquakes and tremors such as tornillo and spasmodic events (Lin et al. 2005a, b), heartbeat-like seismicity (Lin 2017b), VLP events (Lin and Pu 2016), and seismic events generated by the phreatic eruption (Lin 2017a). Thirdly, a shallow pressure source beneath the Dayoukeng, determined by leveling survey (Murase et al. 2014), has been detected in and around the low velocity area (Fig. 10). The consistent location of both the low-velocity zone and the pressure source suggests the presence of a material rich in fluids, which would decrease the seismic velocity. Those fluids and volcanic heat might be released from the suspected magma reservoir in the crust and ascend along the fault zone to the hydrothermal reservoir in the uppermost crust beneath the TVG. Finally, a low resistivity zone beneath the Dayoukeng area and above the low-velocity zone was detected by the magnetotellurics (Komori et al. 2014). Such a low resistivity zone strongly indicate a hydrothermal activity beneath the Dayoukeng area. In summary, all of different observations consistently show a major hydrothermal reservoir exists in and around the Dayoukeng area even although its detailed geometry is still not well constrained yet.

Fig. 10
figure 10

Integrated locations of different observations in and around the Dayoukeng area. a Locations of Shanchiao fault, Mt. Chihsin, Dayoukeng, pressure source (PS, double rings), low resistivity zone (LR, thick-dotted) and a low velocity zone (thin a doted) by Huang et al. (2017). b Distribution of seismicity (pluses), PS, LR (green area), low velocity zone 1 (LVZ1, pink area) in this study and low velocity zone 2 (LVZ2) from Huang et al. (2017) along the NW–SE profile

The discovery of the major hydrothermal reservoir by the LVZ beneath the strong degassing fumaroles indicates the potential for phreatic eruptions taking place in the future. In fact, some small-scale phreatic eruptions have been detected recently (Lin 2017a). Strong acoustic waves released from the fumarole sites such as Hsiaoyoukeng, Huangshan, and Chungshenlo have been recorded by both the infrasonic and seismic sensors in the TVG. Also, several fracture zones, which compose a group of small craters or explosive sites, can be identified through a high-resolution topographic map such as produced by LiDAR or even in Google Earth. Each small crater or explosive site could have resulted from a phreatic eruption in the past. As a result, the phreatic eruptions not only occurred in the past but are also likely occur in the future because the major hydrothermal reservoir provides a significant amount of heat and fluids beneath the fumarole sites at the TVG.

Although the hydrothermal reservoir could have a volcanic impact in the future, it could also provide a sustainable geothermal resource in Taiwan. The major energy supplies in Taiwan are largely from traditional power plants run with nuclear or fossil resources. However, the huge amount of CO2 released by the fossil power plants causes global warming. The safety of nuclear power plants has been seriously challenged after the complex disasters of the 2011 earthquake in Japan. Thus, the exploration for geothermal resources could provide opportunities to develop sustainable energy to replace power produced by traditional nuclear and fossil sources in the future.

The LVZ obtained here with some geometrical uncertainties is one of the preliminary results for future study of both volcanic threats and geothermal exploration. Detailed images of its size and location will be obtained soon when several dense seismic arrays with more than 600 Zl and geophones are deployed in the TVG over the next 3 years (Lin 2016). In the interim, a dense broadband seismic array with a station-spacing of ~ 5 km will cover the northern part of Taiwan for imaging the magma chambers beneath the TVG, as well as Turtle Island (Kueishantao), which is offshore the Ilan plain of eastern Taiwan.

7 Conclusion

A low-velocity zone with some geometrical uncertainties is identified by the deployment of a dense linear geophone array across the Tatun volcano group of Taiwan. A detailed examination of several felt earthquakes shows that consistent P-wave delays are recorded at some particular stations of the array. Further forward modeling indicates that there is a low-velocity zone at depths of between ~ 0.5 and ~ 2.5 km beneath the major fumarole sites. Although such a preliminary result is not a unique solution and it has some uncertainties, the low-velocity zone is related with results from some of our previous studies. Therefore, we strongly suggest that the LVZ with some uncertainties in sharp might be associated with the major hydrothermal reservoir at the TVG. The identification of the hydrothermal reservoir implies that the potential volcanic threat cannot be totally excluded and there exists the potential for sustainable geothermal resources to be developed in the future. Detailed images of the LVZ and other volcanic structures will be obtained later after dense geophone arrays with more than 600 geophones are deployed in the next 3 years (2020–2022).

8 Data and Resources

Seismograms used in this study were collected by the Taiwan Volcano Observatory at Tatun (TVO) using the instruments provided by the instrument center of the Taiwan Earthquake Research Center (TEC). The data will be released by TEC. Some plots were made using the Generic Mapping Tools (version 4.5.2; URL: gmt.soest.hawaii.edu) and the SAC software (version 101.5; URL: ds.iris.edu).