Natural Hazards and Earth System Sciences Multi-parametric investigation of the volcano-hydrothermal system at Tatun Volcano Group , Northern Taiwan

Abstract. The Tatun Volcano Group (TVG) is located in northern Taiwan near the capital Taipei. In this study we selected and analyzed almost four years (2004–2007) of its seismic activity. The seismic network established around TVG initially consisted of eight three-component seismic stations with this number increasing to twelve by 2007. Local seismicity mainly involved high frequency (HF) earthquakes occurring as isolated events or as part of spasmodic bursts. Mixed and low frequency (LF) events were observed during the same period but more rarely. During the analysis we estimated duration magnitudes for the HF earthquakes and used a probabilistic non-linear method to accurately locate all these events. The complex frequencies of LF events were also analyzed with the Sompi method indicating fluid compositions consistent with a misty or dusty gas. We juxtaposed these results with geochemical/temperature anomalies extracted from fumarole gas and rainfall levels covering a similar period. This comparison is interpreted in the context of a model proposed earlier for the volcano-hydrothermal system of TVG where fluids and magmatic gases ascend from a magma body that lies at around 7–8 km depth. Most HF earthquakes occur as a response to stresses induced by fluid circulation within a dense network of cracks pervading the upper crust at TVG. The largest (ML ~ 3.1) HF event that occurred on 24 April 2006 at a depth of 5–6 km had source characteristics compatible with that of a tensile crack. It was followed by an enrichment in magmatic components of the fumarole gases as well as a fumarole temperature increase, and provides evidence for ascending fluids from a magma body into the shallow hydrothermal system. This detailed analysis and previous physical volcanology observations at TVG suggest that the region is volcanically active and that measures to mitigate potential hazards have to be considered by the local authorities.


Introduction
The reliability of any volcano monitoring effort is dependent in part on an accurate understanding of the nature of background activity and eruptive history of the particular volcano. Seismic monitoring is the most commonly used tool for assessing the potential of future eruptive activity (Mc-Nutt, 1996). Geochemical monitoring of fumarolic gases and spring water could also be expected to provide important information regarding changes in the magmatic-hydrothermal system. In many cases the combination of these methodologies has proved very effective in constraining the source of fluids (magmatic versus hydrothermal) and their path through the crust (Moran et al., 2000;Tassi et al., 2004;Madonia et al., 2008;Martini et al., 2010 among others). Additionally, an understanding of the relation between background seismicity, fumarolic chemistry and previous volcanic eruptions can help to establish a baseline against which any future changes in the volcano-magmatic system can be evaluated.
The Tatun Volcano Group (TVG) is about 15 km north of Taipei, the largest city and capital of Taiwan, which has more than seven million inhabitants (Fig. 1). Taipei is a modern metropolis and has a well developed urban infrastructure, such as a high-speed railway system as well as two nuclear power plants located only few kilometers northeast of TVG. In addition to that, TVG itself includes the Yangmingshan National Park which attracts tens of thousands of visitors every year. During the last decade a debate has been initiated regarding the status of TVG as a dormant or extinct volcano (Yang et al., 1999;Song et al., 2000;Lin et al., 2005;Konstantinou et al., 2007;Belousov et al., 2010). The aforementioned studies provide several lines of evidence suggesting that TVG can be considered as a potentially active (thus the broader Taipei area and the study area of TVG which is highlighted with a square. The thick lines indicate the Kanchiao (KF) and Chinshan Fault (CF) respectively while the stars represent the locations of two nuclear power plants, (c) the seismic network distribution used in the study is also presented (YM01-YM12). With C we indicate the location of the Chihsinshan edifice and with T the location of the Tayiokeng area. The yellow lines indicate the fault distribution in the area obtained from LIDAR observations (from Konstantinou et al., 2009). presently dormant) volcano. This evidence includes mantlederived Helium in gases/spring waters, volcano-seismic signals similar to those observed in active volcanoes and radiocarbon dating of the latest eruption at about 6 ka.
A detailed study of the volcanoclastic deposits and the morphology of lava flows and domes in TVG has been conducted by Belousov et al. (2010) in an effort to infer the characteristics of previous eruptions. Volcanism in TVG was mostly andesitic, exhibiting both extrusive and explosive activity, with the former being more common in the form of lava domes and lava flows with lengths up to 5.6 km and volumes up to 0.6 km 3 . Explosive volcanism on the other hand, has produced fallout deposits of plinian eruptions, lithic ashfall deposits of vulcanian eruptions and explosive breccias of phreatomagmatic eruptions. The highest (∼ 1120 m) and most well-preserved volcanic cone in TVG is Mt Chihsinshan which along with the nearby Tayiokeng area exhibit the strongest fumarolic activity and highest seismicity in TVG (cf. Fig. 1). Additionally, Mt Chihsinshan was the focus of several cycles of effusive and explosive activity followed by lahars and gravitational collapses during the last 20 ka. Consequently, monitoring TVG becomes a crucial task since an eruption of any type and scale may have serious consequences for the surrounding region.
In this paper we present a detailed analysis of the spatial and temporal variation of seismicity and its relation to rainfall, gas geochemistry and fumarole temperature variations. This work builds upon the earlier studies of seismicity published by Lin et al. (2005) and Konstantinou et al. (2007Konstantinou et al. ( , 2009 as well as the geochemical study of Lee et al. (2008). Initially, seismic signal classification is performed followed by determination of event magnitude and accurate absolute locations. Among the signals that we have recorded are also events caused by oscillations of fluid-filled cracks and we analyzed them using the Sompi method in order to infer possible fluid compositions. These results are then compared with temporal variations of several gas geochemistry indicators and fumarole temperatures published previously. Finally, we discuss possible mechanisms for the occurrence of seismicity and geochemical anomalies in TVG that can be used in order to understand future periods of unrest.

Signal classification and temporal distribution of seismicity
In a previous study, Konstantinou et al. (2007) Fig. 1 inset). Stations YM01, YM02, YM03 and YM04 were equipped with short-period (Lenartz 3D5s) and broad-band (Guralp CMG-40T) sensors while all other stations had only short-period ones. The first step of the data processing was to manually isolate local earthquakes and other volcano-seismic signals from regional or teleseismic events. The selected signals were classified in different groups by using a simple scheme based on the observed frequency content and waveform appearance (McNutt, 1996). This classification approach was also advisable in order to be consistent with the previous work from 2003-2005 and hence allow a direct comparison. Around volcanic areas propagation effects are usually quite severe and can lead to biased classification of signals if there is a large distance between the source-receiver or the signal is recorded to only a few stations. To avoid this drawback in our analysis we used only events that were recorded by more than five stations. From the classification four groups were created: (a) high frequency earthquakes (HF) with clear Pand S-phases and frequencies 1-20 Hz (Fig. 2a), (b) mixed frequency earthquakes with clear high-frequency (> 20 Hz) P-but no S-phase followed by a low-frequency (∼ 5 Hz) coda ( Fig. 2b), (c) low frequency earthquakes (LF) with emergent P-phase and no S-phase and frequencies up to 5 Hz (Fig. 2c), (d) spasmodic bursts, which are variable duration sequences (total durations between 1 min to 1 h) of HF events placed closely in time (Fig. 3); they most probably represent rupture through a fracture mesh driven by transient variations of fluid pressure (Hill and Prejean, 2005).
In order to keep the same quality level of observations with previous studies, the HF events that were selected and manually picked were recorded by at least five stations yielding a total of 486 events. Also, recorded to five or more stations were 3 LF events, 2 mixed frequency events and 10 spasmodic bursts. Figure 4a shows the temporal distribution of seismicity observed at TVG from October 2003 until end of December 2007. At this point it should be noted that the seismic network experienced several technical problems throughout this period, resulting in about three stations out of eight being down for substantial periods of time. The most serious of these periods was between July to November 2005 where only two to three stations were operational; therefore, it was not possible to have a clear picture of the seismicity in the area. Even though the average number of HF events was about 12 events per month, in several instances this number increased to 40 or even 50 events per month. This increase was either caused by the occurrence of many HF events through a whole month, or the occurrence of one long-duration spasmodic burst within a few minutes.
Rainfall levels can influence seismicity rates around volcanoes through saturating the highly fractured upper crust with fluids, leading to increased pore pressure thus promoting rock failure (Richter et al., 2004). This may be especially true for regions like Taiwan that are subject to a very humid climate throughout the year. We investigated this possibility by checking the monthly rainfall levels at ANBU which is the closest (∼ 2 km) station to TVG and the temporal variation in mm of rain can be seen in Fig. 4b. For the observed period, the average level of rainfall was about 460 mm. During all instances when the seismicity peaked (October 2003, February 2004, April 2006, February and April 2007) the rainfall was close or even below this average. Even though there does not seem to be any obvious link between seismicity (HF events, LF events or spasmodic bursts) and heavy rain in TVG, we should note that the quantity of available meteoric water plays an important role in diluting exsolved volatiles and in cooling the fumarole fluids (this is referred to as the "buffering effect" in the literature, e.g. Martini et al., 2010). We will analyze this point further when we discuss the relation between rainfall and the temporal variation of fumarole temperature/geochemical parameters.  Table 1).

Duration magnitudes
For the purpose of allowing a uniform magnitude determination in time and avoid problems arising from different instrument types with different responses, we estimated duration magnitudes for all HF events. According to observations performed by Aki (1969) and Aki and Chouet (1975), the duration of the seismograms for local earthquakes is independent of the nature of the seismic source, epicentral distance or regional geology, but is strongly dependent on the strength of the seismic source.
From the Konstantinou et al. (2007) study, we measured the seismogram duration (wherever the noise level allowed it) for all the HF events which had their local magnitude (M L ) estimated. To these events we also added all the events that were located in the vicinity of TVG by Taiwan's Central Weather Bureau (CWB) and had a calculated M L . We plotted the magnitude-duration for all these 150 events ( Fig. 5) and applied the maximum likelihood method in order to obtain the linear regression relationship where τ is the seismogram duration in seconds. According to Havskov et al. (2003), who also tested the same method in a volcanic environment, the wide dispersion of the data can be attributed to both the errors in the determination of the signals duration and in the different noise levels present in the data. We estimate the goodness-of-fit of this relationship to the data by comparing the resulting chi-squared value to the degrees of freedom. These two quantities are very close for our case (degrees of freedom 150 − 2 = 148, chi-squared value is 140) therefore the fit can be considered acceptable (Press et al., 1992). For the events from March 2005 to December 2007, we measured the coda duration in all the stations that an event was recorded and obtained the average value. This value was then used to estimate the duration magnitude. Hereafter, when we refer to magnitude for events after the year 2004 we mean M d . The smallest of these magnitudes was equal to 0.1 while the largest was 3.1 and corresponded to a HF event that took place in 24 April 2006 at 03:00:57 (GMT). Figure 6 shows the magnitude distribution for each of the 4 full years spanned by our dataset. For the year 2004, the distribution peaks at magnitude 1.2, while the other three years the peak is lowered to a value of 1.0 as a result of the deployment of more stations.

Absolute locations
All the selected HF and the spasmodic bursts events were located using a probabilistic non-linear method implemented in the software package NONLINLOC (Lomax et al., 2000). When comparing NONLINLOC with the traditional linear location methods several differences are apparent. First, the nonlinear location algorithm is less sensitive to the problem of local minima in the solution space. Second, for events occurring outside the seismic network, it can provide a much better constraint on the hypocentral depth and its corresponding vertical error estimate. Third, it has been found that NONLINLOC can provide a more realistic assessment of horizontal/vertical error estimates than linear algorithms which tend to underestimate such errors (e.g. Lippitsch et al., 2005). A probabilistic solution can be expressed as a posteriori probability density function (PDF) provided that the errors in the observations (phase picking) and in the forward problem (travel time calculations) are Gaussian. The sampling algorithm selected was the Oct-Tree which uses recursive subdivision and sampling of cells in 3-D to generate a cascade of sampled cells (Lomax and Curtis, 2001). NONLINLOC can be used with either a 1-D or a 3-D velocity model and in this study we use the minimum 1-D velocity model estimated by Konstantinou et al. (2007), following the methodology introduced by Kissling et al. (1994). This model is valid only for the upper 5 km of the crust and for this reason we constrain the velocities at larger depths from the tomographic study of Taiwan published by Kim et al. (2005).
By default NONLINLOC provides two kinds of location estimates for each event: (a) one corresponding to the maximum likelihood point of the PDF, and (b) the Gaussian estimate based on the covariance matrix that is obtained from samples drawn from the PDF. The two solutions are in good agreement only if the complete, nonlinear PDF has a single maximum and ellipsoidal form (Lomax et al., 2000). To assess the errors involved in the location estimates, we estimate the differences between the Gaussian hypocenters and the resulting maximum likelihood. In this way we identify 102 events with distorted PDF for which the confidence ellipsoid estimates are biased; this means that from the total 486 HF events the 384 have well-constrained locations. From 2005 to 2007, we find that the absolute error in the epicenter location is between 0.5 km and 0.4 km and the errors in the focal depth between 1.3 km and 1.02 km. The observed improvement in the error estimates for 2007 is due The spasmodic burst (with their durations), mixed frequency and LF events are also indicated. The light gray bars indicate the events analyzed by Konstantinou et al. (2007). The almost zero number of events from July to December 2004 is due to network problems that did not allow a continuous recording of data at all stations, (b) rainfall levels for station ANBU, the nearest one to TVG. The black dashed line indicates the average level of rainfall (460 mm) for the observed period.  Konstantinou et al. (2007) dataset are plotted as blue circles and the 16 CWB events as red circles (linear correlation co-efficient equal to 0.7). The least-squares linear regression result is shown with the black solid line (see text for more details). to a denser network, since the number of stations operating increased to 11 in that year. Based on the confidence ellipsoid dimensions of the well-constrained events from 2005-2007 we estimated mean absolute error in the epicenter location 0.39 km and in focal depth 1.13 km. Konstantinou et al. (2007) compared the absolute earthquake locations as derived from NONLINLOC with the double-difference relative locations. They found that the two sets of locations are compatible which suggests that the results from NONLINLOC are precise and well-constrained. Figure 7 shows the locations of well-constrained HF events and spasmodic burst events for the period 2003-2007 (thus including 144 HF and 59 spasmodic burst events from Konstantinou et al., 2007). In particular, the largest HF event that occurred during the period of this study on 24 April 2006 was located very near Chihsinshan (121.5538 • E, 25.1616 • N) at a depth of 5.16 km. Events belonging to spasmodic burst episodes are all located to similar locations and depths (0-2 km), as in our previous study of TVG seismicity, close to Chihsinshan and at Tayiokeng area. The locations of the 3 LF events are not as accurate as those of HF events, since their P-phase is usually emergent and there is no clear S-phase. Table 1 shows for each of these events both the maximum likelihood and Gaussian locations which differ especially with respect to hypocentral depth. However, for one of the events (Event #2, Table 1) a deeper source is favored by both solutions. HF earthquake activity forms several clusters around the Chihsinshan edifice and Tayiokeng, while the activity in between these two areas is much more diffuse. Once surface fractures determined using LIDAR data are superimposed on these locations, it can be seen that some of the clusters coincide with fracture branches.

Source properties of the 24 April earthquake
Volcanic and geothermal areas are prone to exhibit seismic sources with anomalous ("non-double couple") characteristics as a result of the interaction between fluids and solid rock (Chouet, 1996). The 24 April 2006 earthquake wascompared to the rest of the seismicity observed at TVG -unusual because of its large magnitude and hypocentral depth. In order to investigate whether this earthquake had an anomalous nature, we perform an inversion for determining its moment tensor. We use only waveform data from three stations (YM01, YM03, YM04, unfortunately station YM02 was not in operation at the time of the earthquake) that have wellcalibrated broadband sensors and are azimuthally distributed around the epicenter. A linear, time-domain inversion method described by Herrmann and Ammon (2002) was utilized for calculating: (a) the pure double-couple (DC) solution, (b) the deviatoric moment tensor (compensated linear vector dipole or CLVD and DC components), (c) the full unconstrained moment tensor (isotropic, CLVD and DC components).
The preparation of the waveform data for the inversion involved the deconvolution of the instrument response and rotation of the horizontal components into radial and transverse with respect to the calculated absolute location of the earthquake. After this the data were filtered in the passband 1-4 Hz using a two-pole Butterworth filter. The choice of this passband is dictated by the small signal-to-noise ratio we encountered when trying to perform the inversions at lower frequencies (0.08-0.1 Hz), resulting in inverting noise rather than the earthquake signal. As the frequency range we chose is quite high, we focus on fitting just the first few seconds of each waveform in a similar way to that done by Madonia et al. (2008) for a small event at Vesuvius. Green's functions were calculated for a depth range 1-9 km with a step of 0.5 km using the method of wavenumber integration and utilizing the minimum 1-D velocity model to approximate the upper crustal structure. The traces were then filtered in exactly the same way and aligned along with the data according to their arrival times.
The overall results for the three different inversions show that the highest fit occurs at depths of 5-6 km which agree well with the hypocentral depth of the absolute location (Fig. 8). The fit is the smallest when a DC constraint is applied to the inversion and it is much higher when deviatoric or full moment tensor is assumed. It is generally expected that a model with more parameters can fit better the observations than a model with fewer ones, as in the case of the DC source and deviatoric/full moment tensor. However, it should be determined whether the model with the extra parameters can fit significantly better in a statistical sense the data. For this purpose we employed an F-test to define which of the three models best represents the seismic source (see for example Dreger et al., 2000). The F-test compares the variances of the misfits for a pair of models and accepts or rejects the null hypothesis that they come from the same distribution at a given confidence level. We found that for a confidence level of 95 %, the deviatoric moment tensor fits the data better than the DC solution, but the fit is not significantly better if a full moment tensor is assumed. The preferred solution exhibits a large CLVD (∼ 72 %) and a smaller DC (∼ 28 %) component, while the focal mechanism geometry indicates thrust with a strike-slip component.

Complex frequencies of LF events
It is generally thought that low-frequency seismicity reflects pressure fluctuations resulting from unsteady mass transport (flow) of the subsurface magmatic and or/hydrothermal fluids, hence provide a glimpse of the internal dynamics of the volcano-hydrothermal system (Chouet, 1996). The characteristic properties of the resonator system at the source of the LF event may be determined from the complex frequencies of decaying harmonic oscillations in the tail of the seismogram (Kumagai and Chouet, 1999). This type of seismicity has been linked with eruptive activity and fluid transport processes in many volcanoes (e.g. Nakano et al., 2003;Matoza and Chouet, 2010). Kumazawa et al. (1990) developed a spectral analysis method named Sompi that can quantify the spectral properties of harmonic signals. The method performs a spectral analysis based on a homogeneous autoregressive (AR) model, where a signal is represented as a linear combination of coherent oscillations which are called wave elements and exhibit decaying amplitudes plus some white noise. This kind of analysis can then be used to determine the complex frequencies (frequency and quality factor Q) of the decaying oscillations of the signal under study. The results are presented in a diagram of the frequency (f ) versus the growth rate (g), defined as g = −f/2Q, where the complex frequencies for different trial AR orders are plotted. For more details about the method, the reader is referred to the original papers by Chouet (1999, 2000).
This methodology was applied to the decaying part of the waveforms of three LF events that were identified and located previously. The complex frequencies of the wave elements for AR orders ranging 4-60 were estimated. As the Sompi method also calculates amplitudes of the wave elements, we have added this information to the frequency-growth diagram by using a color scale. This additional parameter can help constrain better the values of Q since high-amplitude clustered wave elements are considered dominant spectral components of the LF signal, while low-amplitude scattered wave elements represent most likely noise (Alparone et al., 2010). Figure 9 shows two of the three LF event waveforms that were analyzed, their Fourier amplitude spectrum and the results of the Sompi analysis. All the LF events analyzed exhibit high-amplitude, clustered wave elements having Q values ranging from 160 to 400. According to Kumagai and Chouet (1999) oscillation damping decreases for ash-laden or water-droplet gas mixture with decreasing weight percent gas. Therefore, the presence of few weight percent gas in the Fig. 9. Examples of results obtained after the application of the Sompi method to LF events #1 and #2 (see Table 1). Top panel: velocity waveforms where the gray shaded part is the one used for the calculation of the complex frequencies. Middle panel: Fourier amplitude spectra of the time windows shown above. Lower panel: Growth rate -frequency diagrams where circled areas indicate high-amplitude wave elements according to the color scale shown at the right (see text for more details). Solid black lines represent constant Q lines for a range of values (20-1000). fluid mixture can generate oscillations with larger Q like the ones observed here.

Gas geochemistry and seismicity
The fumarolic gases of TVG have been analyzed for their chemical composition and isotopic components by various researchers previously (Yang et al., 1999;Lee et al., 2005;Ohba et al., 2010). Recently, Lee et al. (2008) performed regular sampling in 11 sites from the TVG fumaroles using the Giggenbach bottle technique. They published the temporal variations of several gas geochemistry indicators starting from October 2003 up to and including December 2006 (Fig. 10). The meaning of these indicators can be understood if we take into account that SO 2 and H 2 S are the dominant sulfur species in volcanic gases and that they usually have magmatic origin. Additionally, the authors also performed for the same period measurements of the temperature of fumarolic gases even though these were not as regular, thus some gaps exist in the resulting time series.
In order to be able to compare the temporal variations of gas geochemistry and seismicity, we calculate the normalized seismic energy E * of all HF events which is given by (De where d is a constant, M is the magnitude of each earthquake and M min is the completeness magnitude. E * represents the energy release normalized to the energy of the event with magnitude M min so that E * = E/E min = 10 c+dM /10 c+dM min . A value of 1.0 is selected as the completeness magnitude threshold based on the magnitude distributions shown earlier. We also use two values for the constant d, one proposed by Gutenberg and Richter (1956) equal to 1.5 and the other d = 1.96 suggested by Kanamori et al. (1993). Figure 10 11. Possible configuration of the volcano-hydrothermal system beneath TVG suggested by Konstantinou et al. (2007). Main elements of this model are the magma body at depths 7-8 km and the highly cracked, fluid-saturated rock where HF seismicity occurs. Hsiaoyiokeng is an area adjacent to the Chihsinshan edifice which also exhibits strong fumarole activity fed by a system of fractures. be seen to extend from January to August 2005. The geochemical measurements seem to vary very little during the years 2004-2005 and do not correlate well with the normalized seismic energy distribution. On the other hand, the peak in April 2006 when also the largest HF earthquake occurred, coincides with increases in four geochemical indicators and an increase above 100 • C in fumarole temperature. The ratio SO 2 /H 2 S and HCl show a steady increase during late March and early April signifying that the source of the fluids became deeper and their nature was more magmatic rather than hydrothermal. Finally, it should be noted that all geochemical indicators as well as temperature readings after May 2006 exhibit a period of fluctuating values shifting from high to low values and vice versa. This is accompanied by a similar trend in normalized seismic energy mostly due to the occurrence of long-duration spasmodic bursts.

Dynamics of the volcano-hydrothermal system
A model for the volcano-hydrothermal system in TVG has been previously proposed by Konstantinou et al. (2007) based on available geological and seismicity observations (Fig. 11). This model has two main elements, namely the existence of a magma body at depth and on top of it a zone of brittle rock where all HF/spasmodic burst events occur. The depth of the magma body is inferred from geochemical and thermal observations to be 7-8 km, while the brittle zone extends from the surface down to about 5 km where after this depth seismogenesis ceases. The frequent occurrence of shallow spasmodic bursts indicates that the top 2 km in the TVG area is not only brittle, but also highly cracked and fluid-saturated. Such a suggestion is supported by both a high V p /V s ratio (∼ 1.8) and an average shear wave anisotropy of 2.9 % for these depths (Konstantinou et al., 2009). In general, LF events appear only occasionally at TVG and are of shallow origin (< 2 km) even though their hypocentral locations are not as well constrained as for the HF events. Their source process likely involves fluid pressurization/resonance processes inside cracks, while the fluid composition indicated by the Sompi method points consistently to a dusty or misty gas (Kumagai and Chouet, 1999). The lack of any observed sign of volcanic activity points however, to the possibility that the fluid corresponds to water droplet-gas mixture. The quality factors inferred for the LF events analyzed here, are somewhat smaller than those estimated in our previous study of TVG seismicity (Q ∼ 200-1000(Q ∼ 200- , Lin et al., 2005. This change could reflect an increase of the weight percent gas (CO 2 or SO 2 ) in the mixture and give an indication of a variable degassing source. Unfortunately, all the analyzed LF events occurred during February 2007 when there are no additional geochemical observations that could be used in order to corroborate this suggestion.
Even though the source process of HF events is shear failure of rock, the underlying cause for this may be explained by multiple mechanisms, such as regional stresses, stresses due to gravitational loading or due to magma movement. Konstantinou et al. (2009) have performed detailed stress inversions in the TVG area using fault plane solutions of HF events assuming a double-couple source. Their results revealed that NW-SE extension can explain only the deeper (> 3 km) seismicity beneath Tayiokeng and is consistent with the regional stress field in northern Taiwan. On the contrary, the seismicity around Chihsinshan exhibits a localized stress field with subhorizontal σ 1 − /σ 3 -axis which is a configuration that allows opening of cracks and fluid ascent in the upper crust. It has been suggested that HF events can be labeled as "volcano-tectonic" only if they are generated by a local rather than a regional stress field (Moran, 2003). In this context the earthquakes around Chihsinshan represent such events that are caused by stresses induced by movements of volcanic/hydrothermal fluids. In particular for the 24 April 2006 event, it is interesting to note that the azimuth and plunge of the P-axis (37 • , 10 • ) of its DC mechanism is very close to the orientation of σ 1 -axis (azimuth: 34 • , plunge: 7 • ) estimated for the area of Chihsinshan.
The gas geochemistry anomalies observed at TVG may be attributed to three generic mechanisms (Madonia et al., 2008 and references therein): (a) variations of a degassing magma source, (b) changes in rock permeability due to the opening of new cracks, and (c) hydrological modifications due to fluctuations of the quantity of meteoric water in the system. Even though the geochemical measurements during [2004][2005] show no apparent correlation with the normalized seismic energy, it is interesting to note the temperature increase during 2005. The gas temperature reaches about 120 • C in March, April, May, up until June when it drops slightly and increases again in July and August. This variation seems to follow the rainfall levels that are extremely low in April (∼ 100 mm), much higher than average in May (∼ 900 mm) and again lower in June (∼ 300 mm) (cf. Fig. 4b). It is likely, therefore, that this period is characterized by enhanced thermal output from the magma body which provides extra heat to the shallow hydrothermal system resulting in higher fumarole activity. This extra heat may get diffused when cooler meteoric water is of sufficient quantity and mixes with the hydrothermal fluids. In this case the HF events accompanying this thermal anomaly can be interpreted either as the result of the generation of new cracks and/or of transient variations of fluid pressure inside existing cracks.
The temperature of the gas is again increasing during February and March, reaching 140 • C which is the highest value for the whole study period and suggests another episode of enhanced thermal output. Despite the gaps in the measurements, it can be seen that temperature falls to 120 • C in May and this is consistent with rainfall levels that are relatively low in February-March (∼ 300-350 mm) but increase in April and May (∼ 460-520 mm). From late April it is evident that the fumarole fluids become more magmatic in composition (increase of both SO 2 / H 2 S and HCl) even though rainfall levels have also increased and should theoretically favor the dilution of any magmatic components. This apparent contradiction can be explained if the results of the source properties of the 24 April earthquake are taken into account. The combination of a large CLVD component coupled with a smaller DC one can be explained by a tensile crack coupled with an oblate spheroid magma chamber (e.g. Nakamichi et al., 2004). In this scenario, when the crack opens and the chamber shrinks simultaneously, the total moment tensor is equivalent to a large CLVD and smaller DC component with a net volume change equal to zero. Both the travel time location and moment tensor inversion for this event favor a source at a depth between 5-6 km which is close to the top of the inferred magma chamber beneath TVG. We therefore interpret the occurrence of this earthquake as evidence of fluid input to the hydrothermal system of TVG from a deep magma body. Such input is probably responsible for the shift of the fumarolic gas to more magmatic compositions and the subsequent fluctuations of the geochemical parameters between July and November 2006.

Associated hazards
The observations presented in this study suggest that the volcano-hydrothermal system in TVG is active and at times shows strong temporal variations. The HF events can reach magnitudes of 3 with intensities that can be felt in the suburbs of Taipei. Events with magnitudes above 2 can take place very near the surface and be felt in the area around Tayiokeng. From recent physical volcanology observations (Belousov et al., 2010) valuable information can be extracted in relation to the TVG eruptive history. This detailed investigation found that the lava flows at TVG are exceptionally thick and long, with average speeds around 6 m h −1 and the eruptive periods lasting 4-5 yr. Explosive eruptions were in most cases mild, but in Chihsinshan a phreatic eruption occurred very recently (6 ka) followed by gravitational collapse. More than five scars and debris avalanche deposits formed by gravitational collapse were found in TVG and these events occurred long after the eruption had ceased. Lahars of various types and scales were very common due to the very wet climate of Taiwan.
The detailed seismic analysis presented here together with information extracted from other studies should remind people that the TVG region is volcanically active and that measures to mitigate the risks have to be considered by the local authorities. In the hazard analysis we should also take into account that TVG itself includes the popular Yangmingshan national park visited by tens of thousands of people every year and that two nuclear power plants are located very close to its vicinity. For instance, the Kuosheng nuclear power station is located at a distance only 2 km from the front part of one of the largest lava flows (Belousov et al., 2010). Authorities should consider conducting a multidisciplinary volcano monitoring program that will only focus on the TVG area. The design of a public program on volcano hazard education is crucial, even though we note that such a program is something new for Taiwan. This should focus on the issue of living near an active volcano environment and how to take protective measures in case of an eruption. In addition, territorial infrastructure planning and emergency/crisis management plans should be developed for a metropolitan city as densely populated as Taipei.