Three-dimensional inversion of audio-magnetotelluric data acquired from the crater area of Mt. Tokachidake, Japan

Subvolcanic hydrothermal systems can lead to hydrothermal eruptions as well as unrest phenomena without an eruptive event. Historical eruptions and recent unrest events, including ground inflation, demagnetization, and a gradual decrease in the plume height, at Mt. Tokachidake, central Hokkaido, Japan, are related to such a subvolcanic hydrothermal system. This study investigates the three-dimensional (3-D) resistivity structure of Mt. Tokachidake to image its subvolcanic hydrothermal system. A 3-D inversion of the magnetotelluric data, acquired at 22 sites around the crater area, was performed while accounting for the topography. Our resistivity model was characterized by a high-resistivity layer at a shallow depth (50–100 m) and two conductors near the active crater and dormant crater. The high-resistivity layer was interpreted to be composed of dense lava, which acts as a caprock surrounding the conductor. The high conductivity beneath the active crater can be explained by the presence of hydrothermal fluid in fractured or leached zones within the low-permeability lava layer, as the sources of ground inflation and demagnetization were identified within the conductive zone immediately beneath the resistive layer. The resistivity structure was used to estimate the volume of hydrothermal fluid within the pore space. The minimum volume of hydrothermal fluid beneath the active crater that can explain the resistivity structure was estimated to be 3 × 106 m3. This estimate is comparable to the water volume that was associated with the long runout and highly fluidized lahar in 1926. The resistivity structure and volume of hydrothermal fluid presented in this study can be used as a reference for further numerical simulations, which aim to reveal the mechanisms of recent unrest events and assess the risk of hazards, such as lahar.


Introduction
Mt. Tokachidake belongs to the Tokachidake volcanic group, which comprises several volcanic edifices located in central Hokkaido, Japan (Fig. 1). Based on petrological studies, the activity of the Tokachidake volcanic group can be divided into older (1.0-0.5 Ma), middle (300-70 ka), and younger stages (since 60-50 ka) (Ishizuka et al. 2010). Massive dense lava flows, collectively known as the Tairagatake lava (Ta lava), sourced from Mt. Tairagatake in the middle stage, are widely distributed around Mt. Tokachidake. Volcanic activities during the younger stage in the Holocene resulted in the Ground, Suribachi, Kitamuki, Central, and 62 other craters on the northwestern side of the summit of Mt. Tokachidake.
Major historical eruptions of Mt. Tokachidake occurred in 1926, 1962, and 1988-1989(Katsui et al. 1990). These events have all been linked to the hydrothermal system beneath Mt. Tokachidake (Katsui et al. 1963(Katsui et al. , 1990Uesawa 2014). The 1926 eruption was initiated Open Access *Correspondence: r_tanaka@sci.hokudai.ac.jp 1 Institute of Seismology and Volcanology, Faculty of Science, Hokkaido University, N 10 W8, Sapporo 060-0810, Japan Full list of author information is available at the end of the article a c b Fig. 1 Topographic map of the study area and computational mesh of the central part used in the three-dimensional (3-D) resistivity inversion. a Maps showing the location of Mt. Tokachidake (inset) and local topography (red relief image map). b Topography of the area around the core region of the computational domain and spatial distribution of the magnetotelluric measurement sites used in this study. The contour interval is 50 m. The black circles, red stars, green diamonds, blue inverted triangles, and black crosses indicate the 2009, 2014, 2015, and 2016 audio-frequency magnetotelluric (AMT) observations, as well as the 2016 WMT observations, respectively. c Surface mesh around the core region from the Central crater cone, and culminated in the sector collapse of the cone (now known as the Taisho crater) and the subsequent Taisho lahar that killed 144 people at the foot of the volcano (Tada and Tsuya 1927; Committee for Prevention of the Natural Disaster of Hokkaido 1971). Uesawa (2014) suggested that hydrothermal fluids played a significant role in the long runout and highly fluidized nature of the lahar. The subvolcanic hydrothermal system also plays a key role in the most recent unrest event at Mt. Tokachidake. Fumarolic activity at the 62-2 and Taisho craters, along with elevated seismicity in the shallow edifice of Mt. Tokachidake, has continued since the 1988-1989 eruption. Subsequently, since 2000, there has been a gradual decrease in the plume height and temperature of the 62-2 crater. Intermittent localized ground deformation indicates inflation at a rate of 10 4 m 3 /year in a shallow area (several hundred meters deep) beneath the 62-2 crater (Regional Volcanic Observation and Warning Center, Sapporo Regional Headquarters, Japan Meteorological Agency (RVOWC, SRH, JMA) 2018). Annual magnetic total field measurement campaigns, which began in 2008, have revealed continuous demagnetization, indicative of heating at a shallow depth of approximately 150 m beneath the crater . Based on numerical simulations, Tanaka et al. (2017) showed that conduit obstruction, which reduces the permeability, can explain such an unrest episode accompanying demagnetization and inflation, as well as the gradual cooling of fumaroles and declining plume emissions. Moreover, they found that a shallow caprock may be a key structure that controls the time and spatial scales of the ground deformation and magnetic field changes. Takahashi and Yahata (2018), from a geological perspective, pointed out that the dense and less porous Ta lava with a thickness of 500 m, behaves as a caprock.
Based on these findings, characterizing the structure of the hydrothermal system is therefore important, i.e., locating the spatial extent of fluid-rich zones and the overlying caprock, to better understand the processes responsible for the most recent unrest phenomena and evaluate the risks of hydrothermal eruptions and lahars.
As the electrical resistivity is sensitive to temperature, as well as the presence of saline water and certain altered minerals, such as smectite, resistivity is commonly used to investigate the structure of a volcanic hydrothermal system. Several studies on other volcanoes have recently interpreted volcanic activities based on monitoring data in relation to subsurface resistivity structures (Seki et al. 2016;Tsukamoto et al. 2018;Yoshimura et al. 2018). Yamaya et al. (2010) and Takahashi et al. (2017) have also described the resistivity structure of the shallow hydrothermal system at Mt. Tokachidake. Yamaya et al. (2010) imaged the shallow resistivity structure beneath the active crater area through a trial-and-error 3-D forward model incorporating the topography based on audio-frequency magnetotelluric (AMT) data obtained at five sites. They reported a highly conductive zone below the 62-2 crater. However, the horizontal boundary, as well as the lower limit of the conductor, remains poorly constrained owing to the sparsely distributed and small number of measurement points. Takahashi et al. (2017) later improved the resistivity structure via a 3-D inversion with 17 additional sites. However, their model did not account for the topography. The resistivity of rock is many orders of magnitude lower than that of the air. Therefore, topographic highs act as relative conducting bodies compared with the surrounding air space, which can result in a distortion of the induction vectors (Parkinson 1962) surrounding topographic undulations. Such a topographic effect should be incorporated into magnetotelluric (MT) modeling for volcanoes with steep topography (Nam et al. 2007).
In this study, we delineate the shallow hydrothermal system of Mt. Tokachidake based on 3-D resistivity inversion modeling incorporating the topography. The resistivity structure is used to interpret geophysical monitoring data, as well as to estimate the volume of hydrothermal fluid beneath the craters. This information can contribute to a better understanding of the subsurface processes' mechanisms, forecasts of the volcanic activity sequences, and assessments of possible hazards, such as syn-eruptive spouted type lahar.

Magnetotelluric data acquisition
AMT surveys were conducted at 22 sites around the active craters of Mt. Tokachidake in 2009, 2014, 2015, and 2016. Yamaya et al. (2010 used the data acquired in 2009 for their forward modeling. Furthermore, Takahashi et al. (2017) used all of the data for inversion in their modeling. In this study, we selected the AMT datasets obtained from 19 sites (Fig. 1b, c) after excluding bad data points and subsets overlapping with the results of the wideband MT (WMT) survey described below. Two components of the horizontal electric field and three components of the magnetic field were measured using the MTU-5A system (Phoenix Geophysics Ltd.). Remote reference processing (Gamble et al. 1979) was applied to reduce local noise using the reference data obtained approximately 3 km northwest of the survey area. Impedance responses in the frequency range of 10 400-0.35 Hz were calculated from the time series data through frequency analysis using the SSMT2000 software (Phoenix Geophysics Ltd.). The WMT data points were also included for three sites (Fig. 1b, c), acquired in 2016, to constrain the deeper part of the hydrothermal system. In the WMT surveys, ADU-07e instruments (Metronix Geophysics Ltd.) were utilized to record the same five components in the time series as the AMT surveys. Similar to the Phoenix system, the remote reference processing was applied to the data sampled at 32 and 1024 Hz using the reference data set acquired at a continuous recording station (located approximately 600 km from the target area) provided by GERD Co. Ltd. The data set sampled at 3.2 kHz was applied to the cross-reference processing of the data series recorded on the same day at two sites. Impedance responses in the frequency range of 12 288-0.000488 Hz were calculated from the time series data through frequency analysis using the BIRRP code (Chave and Thomson 2004).
As data acquired in four different years were used in this study, verifying possible changes in the subsurface structure over these time periods was necessary. For this reason, the surveys were performed repeatedly in different years at two sites, i.e., at E03 (AMT surveys operated in 2009, 2014, 2015, and 2016) and C07 (AMT and WMT surveys operated in 2009 and 2016). No significant differences were observed among the data series acquired during each survey at E03 and C07. Therefore, all data were inverted together, assuming an absence of changes in the resistivity structure that could be detected by the MT method. In the following inversion process, the bestquality data acquired in 2014 and the WMT data were selected as the datasets for E03 and C07, respectively.

3-D inversion
3-D inversion was performed using the ModEM software (Egbert and Kelbert 2012) to construct a suitable resistivity model considering the topography. The computational domain (20 × 20 × 35 km: NS, EW, UD) was gridded in rectangular coordinates; the smallest horizontal mesh size was 40 m for the central target area while the smallest vertical mesh size for the upper part was 10 m. The grid increased in size toward the model borders, resulting in a mesh of 44 × 46 × 90 blocks. The topography was replicated using the 10-m mesh digital elevation model developed by the Geospatial Information Authority of Japan. For blocks corresponding to air, a fixed resistivity of 10 6 Ωm was utilized. Figure 1b, c shows the plan view and 3-D model of the grid surrounding the central part of the computational area, respectively. The observation sites were relocated by up to 40 m to the center of the block, closest to their original positions, to avoid any instability caused by the boundary between the ground and air blocks. Both the real and imaginary parts of the impedance tensors at 19 AMT sites (13 frequencies between 900 and 0.7 Hz) and three WMT sites (14 frequencies between 256 and 0.03 Hz) were inverted. Error floors of 5 and 10% were used for the off-diagonal and diagonal components, respectively, to avoid overweighting during the inversion process.
The inversion performed in this study consisted of two-steps. The first inversion was conducted using the initial model with a uniform resistivity value for land. In the second inversion, the initial model was derived from the results of the first inversion. We reset the damping parameter (Eq. 1 in Egbert and Kelbert 2012) at the start of the second inversion to avoid falling into local minima. The inversions for the resistivity values of the first initial model were conducted and the best initial value (100 Ωm) was determined from the minimum normalized root mean square (nRMS, Egbert and Kelbert 2012) of the final model. As a result, the nRMS value decreased from 37.4 to 2.07 after 129 iterations of the first and second inversions. (2) An underlying conductive zone (C1, 10-20 Ωm) exists from the ground surface to a depth of 700 m (1000 m a.s.l). The center of C1 appears immediately beneath the 62-2 crater, but the conductive zone extends northward up to Taisho crater. The ceil of C1 partially touches the ground surface at these craters (panels (b), (f ), and (g)-(i)). (3) Another conductive center (C2, < 20 Ωm) is located at a distance of ~ 300 m northeast of the 62-2 crater. The C2 conductive zone lies beneath the Ground crater and has approximately the same dimension as C1. C2, unlike C1, is accompanied by a relatively thick (100-150 m) cover of a surface resistive layer, especially at the Ground crater (panels (c)-(e), (g), and (h)). Comparing the resistivity structure obtained in this study with that in Takahashi   (2017), there are the following differences: (1) in the resistivity structure obtained in this study, C1 extends to a shallower area and is imaged with lower resistivity value; (2) C2, which was only vaguely imaged in Takahashi et al. (2017), is imaged as the other low resistivity structure. Figure 3 compares the apparent resistivities and phases at the selected sites calculated from the field data and inverted final model (Additional file 1: Figure S1 presents comparisons of the responses at other sites). The field data and inverted model exhibit good agreements in terms of these responses.

Results
Additional file 1: Figures S2 and S3 show the results of the sensitivity checks in the vertical direction for C1 and C2. Additional file 1: Tables S1 and S2 list the results of the sensitivity checks in the horizontal direction for C1 and C2. Although it is possible to measure the sensitivity by a sensitivity matrix (Schwalenberg et al. 2002;Usui 2015), we determined the sensitivity by observing the changes in the nRMS while varying the resistivity in and around the conductors. The sensitivity in the vertical direction was measured by changing the resistivity in each conductor dividing it into three parts. We divided the regions in the following way. (1) As shown in Fig. 2, C1 and C2 have their resistivity peaks on each. Therefore, we divided the two in the middle. (2) The external boundary of the two conductors was assumed as 30 Ωm.
(3) Layers were divided so that the thickness of regions 1 and 2 was to be approximately 200 m. Additional file 1: Figures S2 and S3 show the region where the resistivity is changed and the results of the nRMS calculations for the sensitivity checks of the conductors C1 and C2, respectively. Note that each part was changed including deeper parts, that is, when we changed the resistivity of Region 1, we also changed that of Regions 2 and 3. While the nRMS values are sensitive to the resistivity changes in Regions 1 and 2, the sensitivity for the region below 1,450 m a.s.l. (above sea level) (region 3 in Additional file 1: Figure S1) is considered low. Therefore, the existence of a part of the C1 conductor above 1,450 m a.s.l. (regions 1 and 2 in Additional file 1: Figure S2) can be established based on high sensitivity. Similarly, the sensitivity check suggested the existence of a part of the C2 conductor above 1,200 m a.s.l. (regions 1 and 2 in Additional file 1: Figure S3), which can be established based on the high sensitivity.
The corresponding resistivity change was assumed to be negligible if the change in nRMS was less than 0.1; hence, the resistivity range can be characterized by ± 20 and ± 50% as imaged by the final model in Regions 1 and 2 of the C1 conductor and − 20 to + 40% and from − 40 to + 100% in Regions 1 and 2 for the C2 conductor, Fig. 3 Apparent resistivity and impedance phase determined via observations (symbols) and calculations (solid lines). Apparent resistivities (upper panel) and impedance phases (lower panel) obtained for the five sites in the study area are shown. The off-diagonal components, Z XY and Z YX , are shown in red and blue, respectively. The diagonal components, Z XX and Z YY , are shown in green, and pink, respectively respectively. As the typical resistivity of the two conductors in the final model was 10-20 Ωm, the resistivity range of C1, determined by considering the sensitivity (± 50%), was 5-30 Ωm. In a similar manner, the resistivity range of C2, determined by considering the sensitivity (− 40 to + 100%), was 5-40 Ωm. The sensitivity in the horizontal direction for the conductors was measured by changing the horizontal boundary of the shallow part (regions 1 and 2 in Additional file 1: Figures S2 and S3) in C1 and C2. Specifically, the change in the nRMS was observed by shifting the boundary of the conductors through several grids to the east, west, south, and north. The boundary was determined based on 30 Ωm from the result of the sensitivity check described above. The boundary shift was performed by replacing the resistivity of the adjacent cell to the boundary with the resistivity of the cell considered as the boundary. The sensitivity check in the horizontal direction of conductor C1 (Additional file 1: Table S1) showed that the shift in the boundary to the north, east, and west direction caused a larger change than the threshold (0.1) in the nRMS, whereas the expansion of one grid to the south caused a small change. Therefore, conductor C1 was well constrained in the north, east, and west, and had an uncertainty of one cell (40 m) to the south. The sensitivity check in the horizontal direction of conductor C2 (Additional file 1: Table S2) showed small changes in the nRMS for all directions. Therefore, conductor C2 had an uncertainty of 3-5 cells (150-250 m) in all directions.

Relationship between resistivity model and volcanic activity
Here, we discuss a possible physical entity for conductor C1 beneath the 62-2 crater imaged by 3-D inversion. Figure 4 shows the locations of the inflation and demagnetization sources, as well as the micro-earthquake hypocenters on the cross-sections of our resistivity model. The inflation sources, based on GNSS campaigns, are determined by the dashed ellipses in Fig. 4 (RVOWC, SRH, JMA 2018). Although their locations slightly vary with time, they are located within the conductive zone. The demagnetization source, identified exclusively based on changes observed between 2008 and 2009, is located in the shallow part of C1 at an altitude of approximately 1600 m. Seismicity is low inside this conductive zone and is concentrated in the deep area. These spatial correspondence between the deformation/demagnetization sources and the conductor imply that the low resistivity of C1 (5-30 Ωm) represents pores or fracture filled with thermal water. We consider that the voids are partly interconnected, otherwise the bulk resistivity must be much higher because isolated conductive fluids in a resistive matrix cannot significantly contribute to enhancing the electrical connectivity. Another possible candidate that can explain the low resistivity is surface conduction due to clay minerals, such as smectite (Takakura 2009). However, according to microscopic observations and altered mineral assemblages in rock samples reported in Takahashi and Yahata (2018), the hydrothermal environment beneath the active craters of Mt. Tokachidake can be characterized as a high-temperature magmatic-hydrothermal type accompanied by the intense leaching of minerals other than silica due to the percolation of sulfuric acid thermal water. As the altered products are mainly composed of sulfates (Takahashi and Yahata 2018), clay minerals, such as smectite, are unlikely to occur in conductor C1. Therefore, the low resistivity of C1 is most a b Fig. 4 Vertical cross-section of the three-dimensional (3-D) resistivity structure with the corresponding sources of geophysical phenomena. a Cross-section of the study area in the E-W direction at NS = − 250 m. b Cross-section of the study area in the N-S direction at EW = 0 m. The hypocenters located within ± 250 m of each profile are denoted by black dots. The cross symbol indicates the demagnetization source, as identified by Hashimoto et al. (2010). The region enclosed by the black dashed line denotes the area including the ground inflation source located by the Regional Volcanic Observation and Warning Center, Sapporo Regional Headquarters, JMA (2018) likely caused by the presence of interconnected thermal water rather than clay minerals. Although it is difficult to interpret C2 due to the lack of information from other observations, as the resistivity of C2 is almost identical to that of C1, C2 may also contain hydrothermal fluid.
We then interpreted the relationship between the resistivity structure beneath the active craters and recent volcanic activity. A shallow caprock is a key structure for controlling the spatial and time scales of ground deformation and thermal demagnetization as shown by the numerical simulations in Tanaka et al. (2017). The Ta lava geologically corresponds to such a caprock layer (Takahashi and Yahata 2018). The resistivity model presented in this study supports this interpretation. In our resistivity model, the C1 conductive zone, suggestive of interconnected pores filled with high-salinity thermal water, is covered with the superficial high-resistivity layer. We speculate that this layer corresponds to the Ta lava because it is characterized by a high density and low porosity (Takahashi and Yahata 2018). Conductor C1 can then be interpreted as a hydrothermal leaching zone in the middle part of the Ta lava. Leaching of minerals may increase the pores filled with the hydrothermal water and result in the reduction of the bulk resistivity. Such a physically/geologically interpreted resistivity model can also be useful as a reference and/or constraint for the electrical and hydrological properties of the modeled space in future numerical hydrothermal flow simulation studies. Such a simulation study can more precisely investigate the behavior of the hydrothermal fluid in leached zones surrounded by a low-permeability host rock and the resulting ground deformation and demagnetization.

Assessment of hydrothermal fluid volume
The porosity of the leached zones provides useful information to assess a fluid volume in a subvolcanic hydrothermal system. Our resistivity model can be used to estimate the porosity of these zones. Porosity is related to the electrical conductivity based on an empirical formula reported in Archie (1942). Hashin and Shtrikman (1962) gives a theoretical basis for determining the upper and lower bounds of macroscopic parameters for mixtures of different materials using variational principles. Waff (1974) presents explicit representations for the electrical conductivity (i.e., the reciprocal of the resistivity). In this study, we estimated the porosity distribution based on the HSc model (the upper bound of the conductivity proposed by Hashin and Shtrikman 1962), where the pore-fluid is assumed to be perfectly interconnected. Assuming that the pore is filled with only a liquid phase, the HSc model provides the minimum pore volume to explain the bulk resistivity. The volume fraction of the liquid (i.e., porosity) in the HSc model can be expressed as follows: where σ represents the conductivity. The subscripts w, r, and eff indicate the liquid, rock (soil), and effective (i.e., bulk) conductivities, respectively. Here, σ eff refers to our resistivity model. Although we do not have direct measurements of σ r , we assumed a sufficiently small value of 10 -6 S/m for solid rock because this value does not significantly affect the porosity calculation. The resistivity of the liquid was assumed to be 1.2 Ω m based on the resistivity of thermal water sampled at the crater area of Mt. Tokachidake (Takahashi et al. 2017). Figure 5a shows the porosity distribution converted from our resistivity model through this calculation. The region of relatively high porosity (mostly 10-20%) was found beneath the 62-2 crater and Ground crater, corresponding to the high conductivities in C1 and C2. This result appears to be smaller than the laboratory measurements by Takahashi and Yahata (2018), in which the effective porosity of the hydrothermally altered rocks sampled at Mt. Tokachidake range from 10 to 70%, predominantly from 30 to 40%. We speculate that the difference between our estimations and the laboratory measurements results from our assumptions that the intergranular liquid has a uniform resistivity and the pore space is entirely filled with a liquid phase. Figure 5b shows the porosity distribution when the resistivity of the liquid was assumed to be 2.0 Ω m. The estimated porosity of the region beneath the 62-2 crater and Ground crater become larger and close to the measurement of the rock samples. In contrast, although the porosity must be, in principle, less than 100%, the porosity beneath Taisho crater locally exceeds the limit. This is a direct result from the very low bulk resistivity of this particular component in our resistivity model. The empirical equation (Sen and Goode 1992) shows that a 1.2 Ω m value for saline water at 100 ºC becomes 0.7 Ω m at 200 ºC. Assuming that the resistivity of the liquid phase is 0.7 Ω m (Fig. 5c), the porosity beneath Taisho crater does not exceed the limit (i.e., 1.0). However, the porosity beneath the 62-2 crater and Ground crater becomes too low and contradicts the sample measurements. Therefore, a more appropriate analysis requires the spatial distribution of the gas fraction and the liquid resistivity due to the temperature and salinity of the liquid.
The calculated porosity distribution allows the estimation of the minimum volume of hydrothermal fluid that can explain the estimated resistivity structure. Assuming that hydrothermal fluids (1.2 Ω m) exist within the conductor C1, the volume of hydrothermal fluid is estimated to be 5 × 10 6 m 3 . Considering the uncertainty (± 50%) of (1) χ = 3σ w (σ eff − σ r ) (2σ w + σ eff )(σ w − σ r ) , the resistivity value based on the sensitivity check, the minimum volume of hydrothermal fluid is 3 × 10 6 m 3 . As the bulk volume of the region where the porosity is larger than 0.7 is less than 5 × 10 4 m 3 , the porosity of the region was replaced with 0.7 based on the observations in Takahashi and Yahata (2018). In a similar manner, the minimum volume of hydrothermal fluid that can explain the estimated resistivity of conductor C2 was estimated as 8 × 10 5 m 3 . Uesawa (2014) suggests the involvement of hydrothermal fluid (> 1 × 10 6 m 3 ) as a possible catalyst for the long runout and highly fluidized characteristics of the 1926 Taisho lahar. Although what percentage of the hydrothermal fluid beneath the craters was involved with the lahar remains a point of uncertainty, our estimation is comparable to the water volume proposed by Uesawa (2014). Therefore, we speculate that a similar sector collapse in the future, if it occurs, could yield the expulsion of a substantial amount of thermal water from the subvolcanic system, similar to the 1926 Taisho lahar.

Conclusions
3-D resistivity inversion, incorporating topography, was performed for the region beneath the active craters at Mt. Tokachidake. Our resistivity model was characterized by a high-resistivity layer covering a shallow depth and two distinct conductive zones located beneath the