Three-dimensional resistivity modeling of GREATEM survey data from Ontake Volcano, northwest Japan

Ontake Volcano is located in central Japan, 200 km northwest of Tokyo and erupted on September 27, 2014. To study the structure of Ontake Volcano and discuss the process of its phreatic eruption, which can help in future eruptions mitigation, airborne electromagnetic (AEM) surveys using the grounded electrical-source airborne transient electromagnetic (GREATEM) system were conducted over Ontake Volcano. Field measurements and data analysis were done by OYO Company under the Sabo project managed by the Ministry of Land, Infrastructure, Transport and Tourism. Processed data and 1D resistivity models were provided by this project. We performed numerical forward modeling to generate a three-dimensional (3D) resistivity structure model that fits the GREATEM data where a composite of 1D resistivity models was used as the starting model. A 3D electromagnetic forward-modeling scheme based on a staggered-grid finite-difference method was modified and used to calculate the response of the 3D resistivity model along each survey line. We verified the model by examining the fit of magnetic-transient responses between the field data and 3D forward-model computed data. The preferred 3D resistivity models show that a moderately resistive structure (30–200 Ω m) is characteristic of most of the volcano, and were able to delineate a hydrothermal zone within the volcanic edifice. This hydrothermal zone may be caused by a previous large sector collapse.


Introduction
Electrical conductivity of the subsurface which is sensitive to porosity, water content, salinity, temperature and enhance by the presence of clay minerals together with considerable heterogeneity that exists both vertically and laterally are important for imaging and monitoring active volcanoes and hydrothermal systems (Revil et al. 2002). Volcanic edifices tend to have a permeable structure that is readily infiltrated by meteoric water. Hydrothermal systems develop where descending meteoric water encounters hot rock and magmatic fluids (Aizawa et al. 2009). Revil et al. (2010) constructed a 3D resistivity tomography of Vulcan's Forge, Vulcano Island, southern Italy, using the results of DC resistivity measurements along 9 profiles crossing the volcanic edifice. Their results showed a very conductive body (>10 Ω m) inside the Pietre Cotte crater with conductive volumes that are in agreement with the position of temperature and CO 2 anomalies at the ground surface. They interpreted this conductive body as the main hydrothermal system. Auken et al. (2009) used a helicopter-borne transient electromagnetic (TEM) method (SkyTEM) to survey the volcanic terrain of the Galápagos Islands, Ecuador. Their results revealed three resistivity structures: highly resistive basalt, a thin layer with intermediate resistivity on the leeward mountain side which they interpreted as a colluvial deposit and a low-resistivity layer found around the island and at intermediate depth on the windward mountain side. These low-resistivity structures are in the range of freshwater-saturated volcanic rocks.
A number of helicopter-borne time-domain EM systems have been developed in the past decades, such as AeroTEM (Boyko et al. 2000), ExplorHEM (Leggatt et al. 2000), VTEM (Witherly et al. 2004), NEWTEM (Eaton et al. 2004), HeliGEOTEM (Fountain et al. 2005), SKYTEM (Sørensen and Auken 2004) and ORAGS-TEM (Beard et al. 2004). These systems are generally more expensive to operate than fixed-wing EM systems, but offer better spatial resolution. Time-domain methods offer advantages over frequency-domain methods, such as greater depth of investigation and detail, as well as more accurate mapping of freshwater/saltwater boundaries (Steuer et al. 2009).
GREATEM survey techniques have been introduced in the fields of engineering and environmental science, particularly for studies involving active volcanoes (Mogi et al. 2009). To enhance the depth of exploration, we used a large-moment grounded electrical source as a transmitter and towed a magnetic sensor through the air as a receiver. The method increases the survey depth of penetration because it galvanically injects electric current into the ground using a large-moment source and a long transmitter-receiver distance (Abd Allah et al. 2013a).
The Ontake Volcano is more than 3000 m high at its summit, and a grounded electrical-source airborne transient electromagnetic (GREATEM) system is the only way to survey it because the transmitter and receiver are too heavy to fly at such a high altitude using a helicopter, as has been done in other airborne surveys.
Ontake Volcano is a stratovolcano located at the southern tip of the Norikura Volcano Range (Japan Meteorological Agency 1990). The activity of the Older Ontake Volcano began in the middle of the Pleistocene (0.78 Ma) with the eruption of mainly andesitic lava flows. This activity continued intermittently until 0.25 Ma (Kimura and Yoshida 1999). A new eruptive stage began at 0.12 Ma, and the style of the activity changed to eruption of rhyodacitic Plinian pumice, with formation of a dacite lava dome, caldera collapse and a growing stratovolcano.
The previous phreatic eruption occurred in 1979. Subsequently, on September 14, 1984, an M6.8 earthquake, the Western Nagano Earthquake, caused four major landslides and flank failures at and around Ontake Volcano. Despite phreatic eruptions, no activity had been considered to have produced fresh magma during the past 20,000 years (Kobayashi 1993;Kimura and Yoshida 1999), but more recent research shows that four magmatic eruptions have occurred during the past 10,000 years (Oikawa and Okuno 2009). Phreatic eruptions have occurred every several hundred years, and some of these have been large enough to leave deposits. No historical records of phreatic eruptions before 1979 have been discovered, but fumarolic activity has been ongoing for several hundred years at Jigokudani to the southwest of the summit (Oikawa 2008).
At 11:52 a.m. JST (UTC + 9) on September 27, 2014, Ontake Volcano produced another phreatic (steam-type) eruption that sprayed volcanic ash, gas and debris on surrounding areas. Based on the volume of ejecta and the maximum height of the ash column of ~7000 m, the volcanic explosivity index (VEI) of this eruption was 2 (Japan Meteorological Agency 2014). The 2014 phreatic eruption had features similar to those of the 1979 eruption, including the eruption style and the approximate positions of vents, and it is estimated that a pyroclastic flow cascaded down the mountain in a southwesterly direction for a distance of more than 3 km (Kato et al. 2015).
The current study examines the characteristics of Ontake Volcano and its geological setting. Field measurement and data analysis were done by OYO Company under the Sabo project managed by the Ministry of Land, Infrastructure, Transport and Tourism.

Description of survey area
Ontake Volcano is a composite volcano in central Japan (Fig. 1). The edifice contains two overlapping volcanoes with a volume of 80 km 3 , while the caldera is nearly buried (Matsumoto and Kobayashi 1995). The volcanic edifice was built during early and late activity separated by a period of erosion. During the beginning of the late stage, a caldera formed. More recent activity filled in the caldera and the radiating valley, forming the roughly conical shape that exists today (Kobayashi 1993). During the most recent stage of activity, several small andesitic stratovolcanoes, extending from north to south, formed at the summit. Several of their craters are now crater lakes (Mimura et al. 1988). After the caldera was almost completely filled in and the current summit group that stretches north-south formed, the Kiso River lahar deposits formed from a large debris avalanche and debris flow that occurred approximately 50,000 years ago. This flow extended approximately 150 km along the Kiso River. During the last 20,000 years, activity has consisted mainly of phreatic explosions. At the southeastern foot of the volcano, earthquakes have been frequent since 1978 (Kimura and Yoshida 1999). A steam plume has been observed since the eruption of 1979, and a small phreatic eruption occurred in March 2007 (Oikawa 2008).

Geological setting
The basement of the Ontake area is composed of a Jurassic to earliest Cretaceous subduction complex (sandstone, mudstone and chert of the Mino Belt) and Late Cretaceous to Paleogene igneous rocks (felsic welded tuff of the Nohi rhyolites). Contemporaneous with, but spatially distinct from, the Ueno basaltic rocks are andesitic to dacitic volcanics that formed mainly to the east of Ontake Volcano (Ujike and Stix 2000). Takeuchi (1991) concluded that the stress field of the study area was compressional between the late Miocene and the early Pliocene, non-compressional or possibly extensional between 2.7 and 1.4 Ma and moderately compressional after the middle Pleistocene. The Ontake volcanism, which can be divided into active periods (the Older and Younger Ontake Periods) and inactive periods, is summarized by Yamada and Kobayashi (1988) (Fig. 1).
The Older Ontake Period began in the middle Pleistocene. Lavas and volcanoclastic deposits built a large stratovolcano with a volume of erupted products of up to 80 km 3 (Kobayashi 1974). That active period was followed by a long period of inactivity before the onset of activity that formed the Younger Ontake Volcano. The central part of the Older Ontake Volcano was destroyed during subsequent caldera formation in the Younger Ontake Period. The Younger Ontake Period is subdivided into two stages by a sudden change in the nature of the products, without any evident cessation of activity. The upper parts of the cones are composed mainly of dense welded agglutinate with small amounts of lava. A late volcanic unit of the Marishiten Volcano Group (MrVG), Yonnoike Volcano, underlies a widespread fallout tephra and Aira-Tn (AT) ash, whose 14 C age was estimated to be 21-22 ka, including the most recent phreatic explosion of 1979 (Machida and Arai 1983).
Studying one of the earliest products of the younger Ontake Mamahahadake Volcano Group (MmVG), the Pm-I Pumice Fall, Machida and Suzuki (1971) repeatedly determined the fission-track ages of zircon phenocrysts and estimated the age to be ~70-90 ka. A prolonged cessation of volcanism between the two stages has been suggested, and a recent geochronological study using K-Ar dating by Kobayashi (1995, 1999) revealed the detailed history of the Older and the Younger Ontake. The Older Ontake Volcano formed between 750 and 420 ka.

Methods
The method used in this study is the same as that described by Mogi et al. (1998), Verma et al. (2010Verma et al. ( , 2013 Matsumoto and Kobayashi (1995) and Abd Allah el al. (2013a, b). Here, we provide a short description about of our method and its fundamentals. The GREATEM system (Fig. 2) uses a 2-to 5-km-long grounded electrical dipole as a transmitting source and a three-component induction coil magnetometer that measures the time derivative magnetic field in a towed bird as a detector. Mogi et al. (1998) illustrated theoretical transient responses of magnetic fields for horizontally layered structures and noted several features of the GREATEM response, such as the depth of investigation permitted by the large source-receiver distance. They highlighted the advantages of the system for investigating deeper structures and the possibility of identifying resistivity structure responses from sensor heights of 100 m or more, which is important for improving the safety of AEM surveys. With a grounded source, a large-moment source (3.8 × 10 4 a.m.) can be applied, and with a long transmitter-receiver distance, a depth of investigation as deep as 1 km can be achieved in a limited survey area. In addition, GREATEM greatly simplifies the helicopter logistics, as only a receiver bird needs to be flown beneath the helicopter. Other advantages include a nonsignificant amplitude reduction with flight altitude-because the pulse is directly injected into the earth by a ground source-that enables higher-altitude measurements. Thus, the system is safer than conventional heli-TEM systems such as SkyTEM (Sørensen and Auken 2004), VTEM (Witherly et al. 2004) and AeroTEM (Boyko et al. 2000) that carry large and heavy under-slung transmitter loops (sometimes exceeding 30 m in diameter and weighing more than 500 kg) with a terrain clearance of only about 30 m.

Data acquisition
This experiment is part of the larger Sabo project that is managed by the Ministry of Land, Infrastructure, Transport and Tourism to survey the potential risk of landslides on the volcano, and the field operation, data processing and one-dimensional (1D) modeling were conducted by the OYO Corporation under this survey project (OYO Co. 2014). In the current study, helicopter airborne surveys using two sources were conducted over the northern and eastern parts of the Ontake Volcano area with a flight speed of 50-100 km/h, as shown in Fig. 3. The source signals were a time-varying current of approximately 7.2 A transmitted to the subsurface by two grounded electrical dipole sources with a length of approximately 5.3 km. Steel electrodes implanted in the ground to a depth of 20 cm, with 300 electrodes on each side of the cable. The cables were snaked to decrease contact resistance and ensure good coupling to the ground. Three components of the secondary magnetic field related to the subsurface resistivity were recorded in the air by a sensor mounted in the bird being towed by the helicopter. These data were used in the moving and directional corrections, as described in the next section. After these corrections, the corrected vertical time derivative of the electromagnetic field component (Hz) used in both 1D and three-dimensional (3D) modeling was obtained. Magnetic field responses were recorded with the current both 'off ' and 'on. ' The waveforms were digitized through a 24-bit AD converter at a rate of 80 microseconds (µs), and 10,000 EM responses were recorded during one cycle of 0.8 s. The bird was flown at a height of about 100 m, and the measured data were subsequently processed to eliminate noise and obtain a clearer signal.
The survey used two grounded sources, one set to the north and used to survey the northern to central area of the volcano and another set to the east and used to the survey the eastern to central area of the volcano. The northern source had wider coverage than the eastern source in the summit area, and here we use only the data of the northern source. The flight line that received the magnetic field due to the northern source is shown by the blue line in Fig. 3a.

Data processing
The data reduction procedures used in this study are the same as those described in detail by Mogi et al. (2009), Okazaki et al. (2011 and Abd Allah el al. (2013a, b). Their results after data processing were found to be consistent with previous results from other ground geophysical methods in the same study areas, which support the accuracy of the GREATEM data correction and interpretation. Here, we summarize parts of the procedure, as shown in Fig. 4.

Movement correction
Movement correction was performed using a transfer function relating variations in the recorded magnetic field to magnetometer movement, as monitored by a fiber-optic gyro (Japan Aviation Electronics Industry, JCS7401). We used a highly accurate gyro to correct for the movement of the sensor, and this correction was checked in previous surveys (Mogi et al. 2009;Okazaki et al. 2011). The transfer function was calibrated at a time when no signal was applied. Corrections were made by subtracting predicted magnetic field variations using the transfer function based on the movement measured by the gyro from the observed magnetic field variations, yielding data free of motion noise. Figure 4a shows an example of a raw waveform, and Fig. 4b shows an example of magnetic field data before (green lines) and after (red lines) the movement correction. Small transient responses remain but are clearly distinguishable in the motion-corrected time series. The magnetic field data after movement correction at increased resolutions are shown in Fig. 4c.

Coordinate transformation
Transformation of the magnetic field components from bird-based coordinates to geographical coordinates was based on the directional sensor data. The coordinate transformation correction was made using a triaxis orthogonal coordinate transformation. The accuracy of this correction is order of 0.1°.

Removal of local noise
Magnetic field data obtained from the ground magnetometer were used to remove natural and artificial noise, as necessary. Multiple stacking provided the precise waveform induced by the transmitted signal at the ground-monitor site. When necessary, data stacking was used to minimize random noise. Because the helicopter was moving at 50 km/h, or 14 m/s, the sensor moved about 22 m during one 1.6-s cycle. A transient signal after current cutoff was present twice in each cycle, so one transient response could be inferred for each 11 m of survey distance. After stacking n times, the stacked data covered 11 n meters. Three to six transient signals were summed to reduce the noise waveform. Examples of a signal after the 24th stacking and after the 48th stacking are shown in Fig. 4d, e, respectively.

Inversion of assumed layer structure
After the above corrections were made, the noise-free transient response X(t), as shown in Fig. 4f, was inverted to determine a 1D resistivity structure. We used the part of the transient response with less dispersion due to noise and abandoned the later time of the transient curve contaminated with noise using a threshold of order of 10 −3 . The transient response curve of dhz/dt was used to model the waveform. To reduce the error on the curve, we used 4-12 times stacking. The number of data points depends on the noise conditions. Usually, 20-40 data points were collected and spaced at every 44-132 m. The inversion was made by comparing field data with theoretical responses (Ward and Hohmann 1988; Jomori 1991), as shown in Fig. 4g. Usually, 1D inversion is stable if we assume the number of layers and determine only resistivity of each layer. Based on the resolution of the current study eight horizontal layers were assumed and a nonlinear least-squares method (e.g., Strack 1992) was applied. Practically, the depth of investigation depends on how long the clear transient curve is observed. If the transient is buried in the noise, we determine the end of clear transient and the depth of investigation is defined as the depth corresponding to the transient response time at which response data become difficult to fit to a theoretical curve because of scatter caused by noise (Abd Allah et al. 2013a). Figure 5 shows a west-east cross section of the 1D resistivity structure (OYO 2014) along y = 5000 m (a), 6000 m (b) and 7000 m (c). These results reveal Fig. 4 Flowchart of the data processing: a is an example of a raw waveform, and b is three-component time-domain magnetic field data. Lines denote the data before (green) and after (red) motion correction. c The data after motion correction. This shows improved resolution of the magnetic fields shown in the red lines of b. d The stacking effect on the signal after the 24th stacking, and e stacking effect on the signal after the 48th stacking. f Is an example of the final noise-free signal, and g shows a one-dimensional (1D) inversion comparison of the transient field data (red circles) and the calculated, layered-model transient curve (blue line). After Abd Allah et al. (2013b) low-resistivity structures (<25 Ω m) in both the west area between x = ~1500-3000 m and the east area between x = ~5000-7500 m at elevations of ~2000-2500 m. A moderate-resistivity structure (30-200 Ω m) prevails over most of the model and is characteristic of the volcano. In addition, a high-resistivity structure (300-1200 Ω m) prevails over most of the west and central area at y = 5000 m at an elevation from 1000-1500 m, as shown in Fig. 5a.

One-dimensional inversion results
The 1D models based on horizontal layers are adequate in many exploration situations, but there are numerous cases, such as large resistivity contrast boundaries, for which 3D modeling is required (Hördt et al. 1992). In this study, the results of the 1D inversion assuming a horizontal layered structure were found to fit the field data; thus, the 1D results are correct under this assumption; however, subsurface structures are 3D by nature. Therefore, the motivation of this study was to overcome the problem of 1D inversion results by generating a 3D resistivity model that reflects the real dimensionality and tends to have a larger conductors (and conductivity contrast) to explain the observed values, as described in the following sections.

Three-dimensional modeling of GREATEM data
We endeavored to construct a 3D model including topography by forward modeling. A 3D inversion including topography remained problematic and difficult to apply for the present study because the numerical modeling approach was designed for flat surfaces originally, but in the general case, the surface is surrounded by an air layer overlying topography. To overcome this obstacle, we need to apply a special technique to stabilize the solution (Fomenko and Mogi 2002). Large resistivity contrasts between air and topography are a serious issue in surveys of high-relief areas, and horizontally layered resistivity structures might be distorted in this situation. To overcome this problem, we endeavored to construct a 3D resistivity model that considered large changes in lateral resistivity. The 3D EM forward algorithm using a staggered-grid finite-difference method developed by Fomenko and Mogi (2002) was modified for the case of GREATEM by adding a routine for computing electrical and magnetic field components resulting from a grounded finite-wire source ). To compute a finite-length electrical dipole source term, we used the formulae for a long grounded wire on a layered earth (shown as Eqs. 4.184-4.190) described in Ward and Hohmann (1988).
Frequency-domain computation was executed at the frequencies of five equal logarithm spacings in one decade in the frequency range of 10 5 -10 −2 Hz. A timedomain transient response in the time range of 0.1 ms to 1 s was obtained using Fourier transformation.
Using the staggered-grid cell method (Fig. 6a), which computes the tangential components of the electrical field at the edges of the grid, we obtained an accurate result for a high-resistivity contrast model, and topography could be handled as an anomaly in the air layer which is set as a pure insulator with a resistivity of 10 8 Ω m. The z-axis was oriented downward, the center of the ground finite source was located at the point (x, y, 0), and a cable of length L was oriented along the x-axis or y-axis. The modeling area is set sufficiently far from the source line that its effect is negligible. The 3D forward algorithm is accurate and stable through a wide range of frequencies (10 5 -10 −2 Hz) and resistivity contrasts. It is suitable for computer calculation of dhz/dt in frequency domain, which can then be converted to time domain using the inverse Fourier transformation (Newman et al. 1986).

Resistivity modeling of field data
We calculated the EM response numerically to fit GREATEM data from our survey area and used these responses to create 3D resistivity models. As shown in 3D modeling, we chose the area around the volcano summit, which covers an area of about 2.5 × 5 km, as shown in Fig. 3. The total model size was 26 × 30 × 18 km, and a 3D model domain of 14 × 16 × 2.5 km divided into 85 × 90 × 35 grids. In the vertical direction, the thickness of each grid cell varied from 5 m at the surface to 6400 m at the top and the bottom of the model. The horizontal size was approximately 100 m in the modeling area, and large-size (800-6400 m) grids were added at the circumference of the modeling area to avoid edge effects (Fig. 6b). The 3D resistivity model was based on an initial model of 1D resistivity structures inverted from the GREATEM field survey data (Fig. 6c). The 3D forward-modeling scheme developed by Fomenko and Mogi (2002) (Appendix) was used to calculate the response of the 3D resistivity model at each corresponding survey line. The field data were stacked within each cell, and the transient response curves of the field data were fitted to the transient response curve of the 3D-computed data obtained at the center of the cell. We used a trial-anderror method, in which the resistivity values of the initial model were changed repeatedly to obtain a good fit between the 3D synthetic response and the field data. In the field survey, 10,000 samples of data digitized at a rate of 80 μs were obtained during one cycle of 0.8 s, but in the numerical modeling, 40 transient responses computed at 40 time channels and 10 time-equal spaces of a log scale in one decade were obtained in the range of 100 μs to 1 s (four decades), so we used only the data from corresponding time channels (30 samples) and computed the root mean square error (misfit) between the normalized intensity by the maximum amplitude (dhz/dt) of both the field data and modeling responses for different blocks of the 3D model. The root mean square misfit error (RMSE) was defined as [|| ∑(d obs − d cal ) || 2 /N] 1/2 × 10 2 . Figure 7 shows the RMSE obtained after the 5th trial, 25th trial and 50th trial models in comparison with the initial model. As the figure shows, the RMSE misfit improved with repeated calculation of tens of trial models. The best RMSE misfit of 1.1 % was obtained after the 50th trial model. Figure 8 shows examples of fits to the field data for different blocks over the modeling area for RMSE values of 1.1 % (a), 2 % (b), 3.7 % (c), 5.3 % (d), 7 % (e) and 15 % (f ). However, for all of the 3D models, there were a few blocks with RMSE misfit of 15-20 % and others with RMSE misfit of 5-10 %; the best fit to the field data for most blocks was an RMSE misfit of ~1-5 %, as shown in Fig. 9.   Fig. 6 a A staggered-grid cell V (i, j, k). Tangential electric and normal magnetic field components are sampled at the center of edges (E) and at the center of interfaces (H), and all of them are continuous. Δ i , Δ j and Δ k denote the sizes of the cells (Fomenko and Mogi 2002). b Sketch of the grid cells used for the three-dimensional (3D) EM modeling in this study. c The initial model used for the 3D EM numerical modeling computation. The model consisted of a topography layer with different block conductivities and a high-resistivity air layer (10 6 Ω m) Figure 10 shows a bird's-eye view of the 3D resistivity structure model, indicating the locations of the resistivity plan views and cross sections that are shown in Figs. 11 and 12. Figure 11 shows west-east resistivity sections picked from the 3D resistivity model along three profiles. The results reveal that low-resistivity structures (<25 Ω m) prevail in the area x = 5000-6000 m at an elevation of 1800-2500 m and in the area between x = 3500-3800 m at an elevation of 1900-2300 (Fig. 11a). A moderateresistivity structure (40-300 Ω m) covers both flanks of the volcano. A high-resistivity structure (300-1300 Ω m) occurs at elevations of 800-1200 m between x = 4000-6000 m (Fig. 11a) and between x = 5000-6000 m (Fig. 11b, c). Figure 12 shows plan views of the 3D resistivity model along the (a) D-Dʹ, (b) E-Eʹ, (c) F-Fʹ and (d) G-Gʹ sections at different elevations (asl) shown in Fig. 10: z = 1500, 2000, 2300 and 2700 m. These results reveal that at an elevation of 1500 m, the west area is dominated by a moderate-resistivity structure (40-200 Ω m) whereas the east area is dominated by a by high-resistivity structures (>800 Ω m). A highly conductive region  Fig. 11 Northwest-southeast cross section of the 3D resistivity model (Fig. 10), a along the A-Aʹ profile, which corresponds to y = 5000 m, b the B-Bʹ profile, which corresponds to y = 6000 m, and c the C-Cʹ profile, which corresponds to y = 7000 m 1900-2300 m. Overall, these figures show similar resistivity distributions, but with different resistivity boundaries and values. The main difference between the 1D and 3D inversion results is that the 1D computation assumed a layer structure and did not consider lateral resistivity change. Because the computation processes of the 1D inversion assumed horizontal layers that extend to infinity, the conductor thickness and its contrast are too small to recover the observed point-by-point data. The 3D model reflects actual sizes, tends to have larger conductor size and has resistivity contrasts that better reflect the observed values (Abd Allah et al. 2013b). In the 1D model, layers in the volcanic edifice tend to be parallel to the topography. This result was caused by point-by-point inversion from the surface. Layers in the 3D model tend to be horizontal, depending on the grid configuration. Kasaya et al. (2002) conducted magnetotelluric measurements in the frequency range of 384-0.0005 Hz at 25 sites around the central region of Ontake Volcano using the Phoenix V5. Their resistivity model indicated two highresistivity blocks with resistivities higher than 700 Ω m. These blocks were separated by a relatively low-resistivity region in the central part of the model. The shallower regions above the high-resistivity blocks showed lower resistivity of 300-500 Ω m. At depths greater than 8 km, the resistivity gradually decreased. Two conductors were detected in shallower (1-3 km) and deeper portions (>8 km) in the northwestern half of the model. Nurhasan et al. (2006) showed that a good conductor (0.5-20 Ω m) in a volcanic areas can be related to either the presence of hydrothermal water or altered clay minerals. Altered clay minerals can act as a relatively low permeability zones, whereas actively circulating fluids probably exist in fractures. Aizawa et al. (2009) conducted electric selfpotential and magnetotelluric surveys in five large inland stratovolcanoes in Japan, including Ontake Volcano. Their result shows surface geothermal activity localized above the dominant zone of high conductivity, which is inconsistent with the hypothesis that the highly conductive zones consist only of altered clay, without hot saline fluid. Therefore, they interpreted the conductive zones (0.5-20 Ω m) as corresponding to zone of active hydrothermal circulation in which hot saline fluids and altered rocks coexist.

Discussion
Comparing the results of the present study with reasonable interpretations of the results of the previous studies, the following geological features can be clarified: 1. The highly conductive zone (<25 Ω m) can be interpreted as being due to hydrothermal zone that may consist of hot saline water and altered rock within the volcanic edifice. 2. A moderately resistive zone (50-200 Ω m), which can be interpreted as volcanic sediments deposited from erosion of the older Ontake, covers both flanks of the volcano and is characteristic of most of the volcano. 3. The high-resistivity structure (300-1300 Ω m) that covers the model area at elevations of 800-1200 m between x = 3500-6000 (Fig. 11a) and x = 5000-6000 (Fig. 11b, c) can be interpreted as being due to cool, dense and dry unaltered rhyolite.
Low-frequency conductivity depends on pore water conductivity, the outer component of the electrical double layer coating the surface of the minerals, the electrical conductivity of clayey materials, cation exchange capacity and saturation (Revil 2013). To better understand the hydrothermal processes in Ontake Volcano and to what extent the resistivity structure is related to hot saline water or/and altered clay layer, further studies may still be needed.

Conclusions
In the current study, we obtained a reasonable 3D resistivity structure around Ontake Volcano. Our results were able to delineate hydrothermal circulation within the edifice. Taking into account the locations of surface geothermal activity, we interpret the good conductors to roughly correspond to a hydrothermal zone that consists of hot saline water and altered rock. Volcanic sediments with embedded high-resistivity granite and rhyolite cover both flanks of the volcano area and are characteristic of most of the volcano area. The 3D resistivity model can provide more accurate and reliable data, especially in areas of great topographic relief. The GREATEM system can be a useful tool for studying the characteristics of volcanic areas with complex topography to depths of several hundred meters. The information obtained can be used to aid mitigation of natural disasters such as volcanic eruptions and earthquakes and for development and planning in volcanic areas.
In future studies, we will consider about addressing the non-uniqueness of modeling to be able to image the geometry of anomalies both laterally and vertically. This will make our preferred model more robust and well supported.
Authors' contributions SA carried out the 3D modeling and computation and also contributed to preparing the results and the manuscript drafting. TM contributed to the 3D modeling and computation and helped draft and revise the manuscript. Both authors read and approved the final manuscript.
in which µ is the magnetic permeability and σ is the conductivity of a normal grid, ds correspond to the normal grid cells, and ds ′ correspond to the central grid cells. After SFD approximation of Eqs. (1-4), we obtain six scalar equations shown as (5-10) in which The H field is calculated using the curl operator in Eqs. (8-10). Tangential electric and normal magnetic field components are continuous at the staggered grid due to the selection of grid nodes. It helps to stably calculate the magnetic Hz component in the center of cell interfaces X i+1/2 , Y j+1/2 at each level Z = Z k by using Eq. (10).
The original system describing Eqs. (5-10) is a pair of coupled first-order equations for E and H. In order to decrease the number of equations, the H field was eliminated from Eqs. (5-7) using the right-hand side of Eqs. (8-10), and the second-order SFD equation was introduced as follows: (σ ) z i−1,j−1,k = σ i,j,k � i � j + σ i−1,j,k � i−1 � j + σ i,j−1,k � i � j−1 The matrix Ã is neither Hermitian nor symmetric, but it can be transformed to the complex symmetric matrix Â by multiplying the obtained equations by a factor ∆ i for E x members in the right side, ∆ j for E y and ∆ k for E z .