Resistivity mapping in the Tatun Volcano Group , Northern Taiwan , revealed by VLF-MT surveys

In a volcanic region, unstable conditions of the surface geology result in natural hazards such as landslides. The present study aims to examine the properties of the formation near the surface in the Tatun Volcano Group (TVG), Northern Taiwan. Dense Very Low Frequency magnetotelluric (VLF-MT) surveys were conducted to obtain the near-surface resistivity. A geostatistical approach and stochastic simulation were used to examine the continuity of the data, and to render the spatial distribution of the resistivity. The results show that the estimated range is consistent with the length scale of the geological features such as the hydrothermal phenomena and volcanic edifices. The estimated distribution of the resistivity successfully identified a fresh lava (≥ 1000 Ωm), a lineament-bearing, weathered lava (≤ a few hundred Ωm), and hydrothermal areas (≤ 30 Ωm). On the basis of the present study, we expect further detailed geophysical investigations will be performed from the viewpoint of hazard mitigation in the TVG. Article history: Received 26 September 2016 Revised 10 January 2017 Accepted 20 February 2017


INTRoducTIoN
Geophysical exploration has been widely used for delineation of underground structure all over the world, for instance, to assess potential natural hazards (Jongmans and Garambois 2007), to monitor the groundwater conditions (Fitterman and Stewart 1986), and to exploit mineral and energy resources (Waters 1987).Particularly in Taiwan, there are increasing numbers of geophysics-based explorations programs for disaster investigation and mitigation (Chang et al. 2012).Volcanic activity is one of the potential sources of disasters in Northern Taiwan, where there is an active volcanic region, the TVG.
Resistivity is one of the geophysical parameters that can effectively visualize the shallow structure in a volcanic region, because it has a great sensitivity to cracks filled with fluids, rock alteration, and hydrothermal activity, which are characteristics of the region (Stanley et al. 1976;Benderitter and Gérard 1984;Ogawa et al. 1998;Komori et al. 2013a).In the present study, an electromagnetic method was applied to reveal the resistivity distribution in the TVG.The distribution was compared with the geological features in order to examine the distribution of hydrothermal phenomena and the properties of the surface lithology.

sITe descRIpTIoN -TVG
The TVG is located on the northern tip of Taiwan Island, covering an area of ~ 300 km 2 and overlying Tertiary basement, as shown in Fig. 1 (Chen and Wu 1971;Wang and Chen 1990).Significant magmatism occurred in response to the complex tectonic activity caused by the postcollisional processes involving the Philippine Sea Plate and the Eurasian Continent, and the opening of the Okinawa Trough (Wang et al. 1999;Chang et al. 2003;Shyu et al. 2005;Chen et al. 2010).Two major NE-SW trending faults, known as the Chinshan and Kanchiao Faults, were also Terr. Atmos. Ocean. Sci., Vol. 28, No. 6, 833-842, December 2017 developed under the above tectonic regime.The present extensional stress field in the area has resulted in a half-graben structure developing between the faults (Shyu et al. 2005).
The TVG is composed of two volcanic trends; an older NE-SW volcanic trend and a younger E-W trend (Belousov et al. 2010).In the NE-SW trend, the edifice has steep cliffs and deep valleys caused by prolonged weathering.In contrast, the shape of the domes and cones, and their collapses in the E-W trend are well preserved, as represented by Mt.Chishin, Mt.Shamao, Mt.Datun, and Mt.Huangzuei.The domes and cones are composed of weakly to moderately vesiculated, andisitic lava (Belousov et al. 2010).At the foot of these mountains, there are active hydrothermal vents such as Beitou, Siao-you-keng, Matsao, and Da-you-keng, indicating that the TVG has potentially active heat and/or magmatic sources below the surface (Yang et al. 1999(Yang et al. , 2003;;Lee et al. 2005;Ohba et al. 2010;Liu et al. 2011;Ohsawa et al. 2013).In addition, hot springs are also present in the limited area of the coastal strip of the Chinshan District.
Previous electromagnetic surveys have been carried out to investigate the underground resistivity structure in the area (MRSO 1969(MRSO , 1970(MRSO , 1971(MRSO , 1973;;Chen et al. 1998;Chen 2009;Chen and Huang 2014;Komori et al. 2014).In particular, Komori et al. (2014) performed audio-magnetotelluric (AMT) surveys under a thematic research project called "Understanding the Life of the Tatun Volcanic Group: Integrated Analysis of Geophysical, Geochemical and Geodetic Observations" (http://www.earth.sinica.edu.tw/res_project_e.php).The surveys showed the hydrother-mal activity in detail, with a combination of the resistivity structure and other geodetic and geochemical studies.While these studies targeted the hydrothermal system beneath the TVG in a "vertical" direction, the present study will first focus on the "lateral" variations in rock properties (i.e., weathering, alteration, and fracturing) and hydrothermal phenomena near the surface, by performing the dense VLF-MT survey described below.

VLF-MT suRVeys
VLF-MT is based on the principle of magnetotellurics (MT) (Cagniard 1953;Vozoff 1991;Simpson and Bahr 2005), which utilizes an induction phenomenon of VLF waves emitted from the center of a transmitter.The propagating wave principally induces a radial electric field and a tangential magnetic field with respect to the transmitter (Tabbagh et al. 1991), and these are measured by a receiver deployed in the far field region to calculate the apparent resistivity and phase.The VLF-MT is effective in estimating near surface formation resistivity at depths of tens to a hundred meters, making it possible to compare the calculated resistivity with surface geological features (Yamaguchi et al. 2001).Furthermore, VLF-MT survey can be conducted in a very short time (typically less than ten seconds per measurement, depending on data stacking), enabling spatiallydense surveys to be carried out compared to conventional MT observations.
In the present study, the VLF-MT surveys were performed using an electromagnetic wave with a singular frequency of 22.1 kHz, emitted by the Ebino VLF transmitter (Miyazaki Prefecture, Japan), approximately 1200 km away from the study area.The VL-101 equipment manufactured by Tierra Technica Ltd. was used as the receiver.An induction coil and two potential electrodes were deployed to measure the tangential magnetic field and the radial electric field, respectively.Notably, such a configuration is usually adopted in order to maintain a high signal to noise (S/N) ratio (Scott 1975).In addition, the phases were used to check the quality of the data, noting that most of their values (mean: 47°, standard deviation: 13°) were considered to be within the acceptable range.
Figure 1 shows the survey area and the survey points.A total of 324 measurements were carried out in 2009, 2012, and 2013 in the mountain area of the TVG and additionally in the Chinshan District.The measurements were basically done along mountain roads, while avoiding artificial structures and power lines.It is noteworthy that there was no significant change in the volcanic activity during the period (Murase et al. 2014).The calculated resistivity data were analyzed by a geostatistical approach described below, in order to provide a spatial distribution of the resistivity.

GeosTATIsTIcAL dATA ANALysIs
Geostatistical techniques are used to statistically predict unsampled values of variables from known datasets, by considering their spatial correlation (Deutsch and Journel 1998).The techniques have been widely applied to various field data to estimate the spatial distribution of variables such as lithologies (Koike and Matsuda 2005), ore grades (Journel 1974;Chatterjee et al. 2010), pollutants (Atteia and Dubois 1994;McGrath et al. 2004), structural features such as crack density and its orientation (Koike et al. 2015), hydrological parameters such as permeability (Cassiani et al. 1998;Fairley et al. 2003;Heffner and Fairley 2006;Anderson and Fairley 2008), and others.
The spatial correlation of data can be expressed using the "Variogram", and is defined as (Deutsch and Journel 1998): where: c is the semivariogram, Var{} is the variance, Z(u) and Z(u + h) are the random variables at the location u and u + h within A, respectively.Under an assumption of data stationarity (Deutsch and Journel 1998), the variogram c is a function of the separation distance h.In a case of small h, the variables are expected to have good correlation, exhibiting a low value of c, while a large |h| would give a high c.
Therefore, c allows us to evaluate the spatial continuity of the data values.The present study used the indicator semivariogram defined as (Goovaerts 1997;Deutsch and Journel 1998): where: N(h) is the number of data pairs, I(u i ) is the indicator transform, z(u i ) is the known value of the random variable Z at the location u i , and Cut k is the cutoff value.The indicator semivariogram is effective when extracting a specific characteristic, such as lithology, from numerical information (Shoji and Koike 2007;Koike and Shoji 2008).In the present study, nine cutoffs were defined in order to examine the relationship between data continuity and the spatial extent of surface geological features in the TVG.Figures 2a -i show the experimental indicator semivariograms for each cutoff value, noting that the cutoff values correspond to the cumulative frequency "p" of the resistivity shown in Fig. 2j.The variogram models are shown as the solid lines in the figures, while their parameters are given in Table 1, which were objectively obtained from the numerical code "AUTO-IK" (Goovaerts 2009).The accuracy and precision of the models are assured by the accuracy plot and the standardized PI-width plot, as shown in Figs.2k and l, respectively (see Goovaerts 2009).It was found that the last three cutoffs (320, 400, and 1000 Ωm) have ranges, a proxy for the continuity of data, that are about three times as long as the first three cutoffs (16, 32, and 40 Ωm).This indicates that the conductive area has a smaller spatial extent than the resistive one.In fact, low resistivities (< 40 Ωm) and high resistivities (> 300 Ωm) were observed mostly at the hydrothermal phenomena and the edifice of the E-W volcanic trend, respectively.According to the geological investigations (MRSO 1971), the former has a spatial extent of several hundred meters to one kilometer, while the latter has a much larger spatial extent.Therefore, the modeled ranges are quite consistent with the geological features.

sTochAsTIc sIMuLATIoN
The modeled variograms were used for an indicatorbased stochastic simulation processing, which is a tool to estimate the spatial distribution of variables by considering their uncertainties (Deutsch and Journel 1998).In the simulation, an indicator kriging estimates a conditional cumulative distribution function (CCDF) of an unsampled value, using known values of the adjacent points.The unsampled value is determined using the kriged CCDF and the cumulative probability selected by a Monte-Carlo method.The determined value is used together with the other known values for estimating the next unknown values.By iterating the above processes, the spatial distribution of the variables is estimated.The program sisim in GSLIB (Deutsch and Journel 1998) was used for the analyses in the present study.

ResuLTs
Figure 3a shows the estimated distribution of the resistivity, which was obtained from an average of 300 iterations of the simulation (called "E-type" estimation), while the topographical relief is overlain as contours.In the TVG, a high resistivity region (≥ 300 Ωm) is extending in E-W direction from Mt. Huangzuei, through Mt.Chishin, to Mt. Datun.On the other hand, the NE-SW volcanic trend has a moderate resistivity of ~100 Ωm.Several spots with low resistivity (≤ 30 Ωm) are distributed along valleys in the mountains, which correspond to hot spring and fumarolic areas, as denoted in Fig. 1.In the Chinshan District, low resistivities are aligned along the northern side of the Chinshan Fault, and are also present around the Chinshan hot spring area.In addition, the alluvial plain in the Chinshan District maintains resistivities in the range of 50 -100 Ωm.The variance map of the above simulations is shown in Fig. 3b, and could be useful for examining the relative reliability of the estimated resistivities.The variance tends to be small in areas where the surveys were densely performed.On the other hand, there are some blue-colored spots with relatively-large variances, corresponding to the presence of an abrupt change in the resistivity value.
Figure 4 shows a close-up map of the resistivity distribution focusing on an area of Mt.Chishin, which is overlain by lineaments and topographical features derived from a LiDAR digital elevation map (DEM) (Liu et al. 2007;Konstantinou et al. 2009;Belousov et al. 2010).Extremelyhigh resistivity areas (≥ 1000 Ωm) correspond to the lava of Mt.Chishin and Mt.Shamao, one of the youngest volcanic edifices in the TVG.Particularly in Mt.Chishin, the lava flowing in a S-E direction and the avalanche on the western flank are well traced by areas of high resistivity.In contrast, the lineament-bearing lava exhibit a resistivity of ~ a few hundred Ωm.In the area of Siao-you-keng, extremely-low resistivities (≤ 10 Ωm) are present in the western rim of the crater, and correspond to active fumarolic areas.These are also present in the Beitou, Da-you-keng, and Matsao areas.

dIscussIoN
It is widely accepted that resistivity is significantly decreased by clay minerals, saline fluids and enhanced connectivity (i.e., high effective porosity) (Revil et al. 1998(Revil et al. , 2002;;Komori et al. 2010Komori et al. , 2013a)).This section compares the resistivity map with previous geological and geophysical studies in order to examine the properties of the surface lithologies.

Lavas in the TVG
Lava of the E-W volcanic trend is characterized by resistivities higher than 1000 Ωm, and is generally interpreted as fresh, dense volcanic rocks, which have low-porosity and minor clay content due to minimal alteration (Komori et al. 2013b).The result is consistent with the fact that the volcanic trend has a young age (Belousov et al. 2010).On the other hand, the lineament-bearing lava has lower resistivity, because lineaments increase the water content of the lava as well as its connectivity (Shankland and Waff 1974).Besides, weathering promotes the progressively increasing interaction of water and rock alteration that can result in the moderate resistivity of the lava in the NE-SW volcanic trend.

hydrothermal phenomena in the TVG
The low and extremely low resistivities in the hot spring and fumarolic areas are considered to be due to high temperature saline fluids and a large amount of clay minerals.In fact, the lowest resistivity of 0.5 Ωm was measured on the hot spring water in the area of Da-you-keng in the present study.Such low fluid resistivities are supported by several geochemical studies that report that the fluid has a high saline content (~20000 ppm of TDS), high temperature (~80°C), and low pH (~1) (CGS 2009).Furthermore, MRSO (1973) described the surface of the hydrothermal area as being composed of altered clays, including smectite and kaolinite.Therefore, the geological and geochemical aspects of the area can adequately explain the resistivity characteristics of the area.
In the TVG, Ohsawa et al. (2013) and Komori et al. (2014) considered that hydrothermal fluids mix with meteoric water in the shallow groundwater system during their ascent from depths, and then flow along topographical features.Therefore, the low resistivities observed along the valleys in   the present study are considered to be due to the seepage of hot springs as a consequence of topography-driven flow.

chinshan district
Yang and Chen (1989) performed detailed transient electromagnetic (TEM) surveys in the area around the Chinshan Fault, in order to reveal the distribution of onedimensional resistivity structures.The investigators showed that a sediment layer with a resistivity less than a few Ωm occurred at shallow depths on the northern side of the Chinshan fault, due to the southward-dipping normal faulting.In the present study, therefore, the low resistivities found on the northern side of the fault could be interpreted as conductive sediments.

coNcLusIoNs
Lineaments and rock alteration accompanying weathering and hydrothermal activity, commonly destabilize strata, causing landslides (Reid et al. 2001).In the present study, such properties in the rocks were successfully represented by resistivities of below a few hundred Ωm.Therefore, near-surface resistivity should allow us to identify hazardous areas in the TVG.On the basis of the spatial distribution of resistivity rendered here, further geophysical investigations, including electromagnetic observations with multiple frequencies and seismic reflection surveys, should reveal the shallow structural characteristics of the above hazards (Bichler et al. 2004;Chiang et al. 2015).
Today, airborne surveys are becoming the mainstream method of conducting surface investigations, by virtue of their capability to rapidly scan extensive tracts of the Earth's surface with high sampling densities (Siemon et al. 2009;Komazawa et al. 2010;Okuma et al. 2014).However, ground-based surveys are still the more robust tool, because they are overwhelmingly cost-effective and safe, and are barely affected by weather.We expect that our methodology and results can be used effectively for the mitigation of volcanic hazards in the TVG.

Fig. 1 .
Fig. 1.Location map of the TVG, Northern Taiwan, showing representative geological features and mountain peaks.Andesitic lava of the TVG covers the northern tip of the island.VLF-MT measurement points are denoted by solid circles, while the mountain peaks are denoted by red triangles.

Fig. 2 .
Fig. 2. Results of geostatistical data processing using the AUTO-IK program.(a) -(i) Indicator semivarigrams for each cutoff value.Broken and solid lines represent the experimental variograms from the VLF-MT data and the modeled variograms from Table 1, respectively.(j) Frequency distribution of the resistivity plotted as a histogram (grey bars) and as cumulative frequency "p" (red line).(k) Result of the accuracy plot.(l) Result of the standardized PI-width plot.

Fig. 3 .
Fig. 3. Results of the geostatistical data analyses and indicator-based stochastic simulations.(a) Spatial distribution of the interpolated resistivity, with topographical relief shown by overlain contours.(b) Variance map of the simulations.

Fig. 4 .
Fig. 4. Detail map of the resistivity distribution, focused on Mt.Chishin.Topographical features and lineaments traced from the LiDAR DEM are overlain on the distribution.

Table 1 .
Variogram Models of Resistivity Variability at the TVG.