論文 / 著書情報
Article / Book Information
Title

Anatomy of active volcanic edifice at the Kusatsu‒Shirane volcano,Japan, by magnetotellurics: hydrothermal implications for volcanic
unrests

Authors

Kuo Hsuan Tseng, Yasuo Ogawa, Nurhasan, Sabri Bülent Tank, Naoto
Ujihara, Yoshimori Honkura, Akihiko Terada, Yoshiya Usui, Wataru
Kanda

Citation

Earth, Planets and Space, Volume 72, ,Pub. date

2020, 10

DOI

https://doi.org/10.1186/s40623-020-01283-2

Creative Commons

The information is in the article.

Powered by T2R2 (Science Tokyo Research Repository)

(2020) 72:161
Tseng et al. Earth, Planets and Space
https://doi.org/10.1186/s40623-020-01283-2

Open Access

EXPRESS LETTER

Anatomy of active volcanic edifice
at the Kusatsu–Shirane volcano, Japan,
by magnetotellurics: hydrothermal implications
for volcanic unrests
Kuo Hsuan Tseng1, Yasuo Ogawa2* , Nurhasan1,3, Sabri Bülent Tank1,4, Naoto Ujihara1,5, Yoshimori Honkura2,
Akihiko Terada2, Yoshiya Usui1,6 and Wataru Kanda2

Abstract
We aimed to perform three-dimensional imaging of the underlying geothermal system to a depth of 2 km using
magnetotellurics (MT) at around the Yugama crater, the Kusatsu–Shirane Volcano, Japan, which is known to havefrequent phreatic eruptions. We deployed 91 MT sites focusing around the peak area of 2 km × 2 km with typicalspacings of 200 m. The full tensor impedances and the magnetic transfer functions were inverted, using an unstructured tetrahedral finite element code to include the topographic effect. The final model showed (1) low-permeability
bell-shaped clay cap (C1) as the near-surface conductor, (2) brine reservoir as a deep conductor (C3) at a depth of
1.5 km from the surface, and (3) a vertical conductor (C2) connecting the deep conductor to the clay cap which
implies an established fluid path. The columnar high-seismicity distribution to the east of the C2 conductor implies
that the flushed vapor and magmatic gas was released from the brine reservoir by breaking the silica cap at the brittle–ductile transition. The past magnetization/demagnetization sources and the inflation source of the 2014 unrest
are located just below the clay cap, consistent with the clay capped geothermal model underlain by brine reservoir.
The resistivity model showed the architecture of the magmatic–hydrothermal system, which can explain the episodic
volcanic unrest.
Keywords: Phreatic eruption, Magnetotelluric method, Clay cap, Brine, Unrest
Introduction
The Kusatsu–Shirane volcano is an active andesitic–
dacitic quaternary volcano in Central Japan (Fig. 1a, b).
The activity of the volcano started at 0.5–0.6 Ma withmagmatic eruptions and pyroclastic flows. After a longdormancy of 0.2–0.3 Ma, magmatic activity restarted
from 16 Ka and ended at 1.5 Ka (Terada 2018). The subsequent historically known activities of the Kusatsu–Shirane volcano since 1805 were phreatic eruptions around*Correspondence: oga@ksvo.titech.ac.jp
2
Volcanic Fluid Research Center, Tokyo Institute of Technology, Tokyo,
Japan
Full list of author information is available at the end of the article

the peak of the volcano, especially near the Yugama crater (Y in Fig. 1c).
Phreatic eruptions are local-scale volcanic phenomenathat include explosions of confined volumes of gas andsteam, but exclude those directly involved with magma(Stix and de Moor 2018). The cause of phreatic eruptions is related to the changes in transient pressure in thesubsurface hydrothermal systems, which is mainly due
to the intense temperature increase of shallow aquifers(Kobayashi et al. 2018). A phreatic eruption has a comparatively small impact radius; however, the concentrated
gases can be emitted, which spread over a wide area
(e.g., Tazieff 1989), and the mechanical energy of eruptions can also provoke fatal ejecta near the crater (e.g.,© The Author(s) 2020. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) andthe source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material
in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material
is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds thepermitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat​iveco​
mmons​.org/licen​ses/by/4.0/.

Tseng et al. Earth, Planets and Space(2020) 72:161

Page 2 of 11

Fig. 1 a Location of Kusatsu–Shirane Volcano, Japan. b Topography map of Kusatsu–Shirane Volcano with MT/AMT stations. The squares andtriangles denote the MT and AMT stations, respectively. The red circles and open curves show the crater rims of past eruptions. The purple closed
curves show the fumarole zones. The blue curves show the national road 292 as a location reference. SG, MS, and MZ denote the locations ofSesshougawara fumarolic area, Mt. Motoshirane, and Manza hot spring field, respectively. The long orange line denotes a section in Additional file 2:
Figure S23. c A peak area [the red rectangle in (b)] of Kusatsu–Shirane Volcano with MT/AMT stations. The symbols are the same as in (b). Moreover,the black squares with red and with white edges denote the MT station with and without magnetic field data, respectively. KS denotes Mt. Shirane,
and the red triangle denotes its peak. Y denotes the Yugama crater lake

Kaneko et al. 2016) or large lahars (e.g., Herd 1986). It isfundamentally important to know the geothermal system
underneath the volcano edifice to understand the phreatic eruption mechanism and to mitigate volcanic risks.
In this study, we used the magnetotelluric method (MT),an electromagnetic field method, to image the subsurfacehydrothermal system in terms of the electrical resistivity,
which is sensitive to the distribution of fluids and lowpermeability clay cap (e.g., Aizawa et al. 2009; Heise et al.
2008).
The hydrothermal system of the Kusatsu–Shirane volcano has been intensively studied by performing geochemical analysis at the Yugama crater lake and thesurrounding fumaroles and hot springs (Sano et al. 1994;
Ohba et al. 2019). Ohba et al. (2019) showed a model of a
two-phase gas–liquid reservoir underlain by a magmatic

fluid reservoir based on the isotope ration of hydrogenand oxygen in relation to recent volcanic unrest.
During the 1990s, which was an active decade in termsof seismicity, peculiar seismic tremor events were monitored, and underlying horizontal crack models with a
two-phase reservoir were proposed for the system underneath the Yugama crater (Nakano et al. 2003). However, the overall seismic image of the geothermal system
underneath the volcano remains unclear.
Temporal changes due to high-temperature fluid
migration were also detected by total magnetic intensity
monitoring around the Yugama crater during both active
and inactive periods from 1982 to 2012 (Takahashi andFujii 2014).
For imaging the overall geothermal system, Nurhasan
et al. (2006) carried out a 15 km long profiling of

Tseng et al. Earth, Planets and Space(2020) 72:161

audio-magnetotellurics measurements and showed a
two-dimensional resistivity section across the volcano
down to a depth of 2 km along the orange line in Fig. 1b,c. Recently, Matsunaga et al. (2020) carried out magnetotelluric measurements 2 km to the south of the Kusatsu–
Shirane volcano and showed a three-dimensional
structure to a depth of 10 km with 1–2 km horizontal
resolutions. Those studies; however, lack detailed sitespacings and images of the geothermal system aroundthe active peak area around the Yugama crater. In thisstudy, we focus on the peak area of the Kusatsu–Shirane
volcano and will present a three-dimensional resistivity distribution from a densely sampled magnetotelluric
data.

Magnetotelluric data
We compiled a dataset consisting of magnetotelluric(MT) and audio-magnetotelluric (AMT) data obtained
during 2001–2005. We have 91 sites in total (Fig. 1b, c).
See Additional file 1: Table S1 for details. Out of these,
41 sites used an MT system covering a frequency rangeof 320 Hz–0.5 mHz. Most of them measured only telluric
fields onsite, simultaneously with another site with magnetic field measurements (Additional file 1: Table S1).
The remaining 50 sites used an AMT system covering afrequency range of 10 kHz–0.1 Hz, which measured two
horizontal components of telluric fields and three components of magnetic fields at the same locations. For the
MT data which had no magnetic fields (Additional file 1:
Table S1), the impedances were calculated using the telluric fields at the site and magnetic fields at another site,
assuming that the magnetic fields have a smooth spatial
variation. For those impedances using the interstation
data (e.g. Martí et al. 2020), the different locations of
the telluric, and magnetic fields were considered in the
inverse modeling.
In this study, we focus on a 2 km × 2 km area (Fig. 1c)and to a depth of few kilometers in the model interpretations; we also include the NW–SE profile AMT data by
Nurhasan et al. (2006) for modeling. We chose 12 frequencies in a range from 320 Hz to 0.1 Hz, which equally
spaced in logarithm for inversion. This frequency range
could cover the available data in both the MT and AMT
measurements.
The observed data are summarized in Fig. 2a as phasetensor ellipses (Caldwell et al. 2004) for 5 representative
frequencies. The distribution of ellipses directly shows
the directionality and dimensionality of the structure
without distortions due to the near-surface anomalies.
The ellipses were filled with tan−1 (2 ), where 2 is thegeometric mean of the diagonal components of the phasetensor components (max and min ). The 2 phase is arotational invariant and indicates the changes in electric

Page 3 of 11

conductivity with depth. The phases of higher and lower
than 45° respectively indicate the existence of more conductive and more resistive bodies with increasing depth.
The real and imaginary parts of the induction vectors arealso plotted in Fig. 2a for the AMT sites, with Parkinson’s
convention by which the real vectors point to conductors.
We see more circular-shaped ellipses at higher frequencies, which mean almost one-dimensional resistivity distribution at the near surface. At the frequencies
from 70 to 4.4 Hz, the high phase (~ 75°) spreads aroundthe Yugama crater lake (Y in Fig. 2a), and the phases
are higher to the east of the Yugama crater. These indicate widespread shallow conductors. The induction vectors point towards the central peak around Yugama lake,
which also includes topographic effects. At a frequencyof 1.1 Hz, the relatively high-phase area is shifted to thenorth. The shapes of the ellipses are more circular near
the Yugama crater and become more E–W elongated atsites away from the Yugama crater. This feature is an indication of the existence of another deep conductor underneath the Yugama. At a frequency of 0.3 Hz, all the sites
show consistent E-W elongated ellipses, which indicatethat the peak area has a deep conductor in the N–S strikeat depth.

Inverse modeling
We carried out 3-D inversion using impedance data andmagnetic transfer function data from 280 to 0.1 Hz. The
sites with AMT data utilized impedance data and magnetic transfer function data; Conversely, the sites withMT data utilized only the impedance data (Additional
file 1: Table S1). The error floor of impedances was set√
as 5% of 2 |Zxy∗Zyx| (Usui 2015); The error floor of the
tipper was set to 0.05 (see Additional file 2: Text S1 fordetails). As the undulation of the topography has a comparable scale to our target, we used a 3-D inversion codedeveloped by Usui (2015) and Usui et al. (2017), which
can incorporate the sharp topography using unstructured
tetrahedral meshes. It also allows the inclusion of local
galvanic distortions as parameters.
The initial model for the inversion was uniform earth
having 100 Ωm resistivity with topography and ­109 Ωmair. The Yugama crater lake water was modeled with
0.2 Ωm resistivity after the water sample measurement(Terada et al. 2018). The resistivity of the lake water and
air was fixed during the iterative calculation. The modeled horizontal electric fields were calculated in such a
way that the electric fields were tangential to the groundsurface.
The inversion program FEMTIC minimizes the data
misfit along with the two constraints, namely the modelroughness and the norm of galvanic distortions. Thus,
there are two tradeoff hyper-parameters: α 2 for model

Tseng et al. Earth, Planets and Space(2020) 72:161

Page 4 of 11

Fig. 2 The phase tensor ellipses map with a the observed and b the calculated data from the final model. The phase tensor ellipses are filled by­tan−1 2, where 2 is the geometric mean of maximum and minimum of phase tensor. The blue and green arrows show the real and imaginaryparts of the induction vectors with Parkinson convention. The black curve, red triangle, and blue circle with Y denote the national road 292, the peakof Mt. Kusatsu–Shirane and the Yugama crater lake, respectively

roughness and β 2 for the norm of galvanic distortions.
We chose a value of β 2=1 to adequately control the galvanic distortion. For determining the best model, we
tested seven different values for α 2 (100, 31.6, 10, 3.16,
1, 0.316, and 0.1) and found the best model with α 2=1

using the L-curve method (Additional file 2: Figure S1).
With the preferred tradeoff parameters, the model converged at the 20th iteration with the root-mean-square(RMS) changes of less than 1% from that of the previous iteration. The RMS decreased from 32.5 to 1.24 from

Tseng et al. Earth, Planets and Space(2020) 72:161

Page 5 of 11

the initial model to the final model. Comparisons of theobserved and calculated phase tensors are presented inFig. 2a, b for five representative frequencies. Both thephase tensor ellipses and its tan−1 (2 ) phase value show
good agreement. The induction vectors also had a reasonable fit. Additional file 1: Table S1 shows the locationsof the MT and AMT sites, and the Additional file 2: Figures S2–16 show the comparison of observed and calculated responses for apparent resistivity, impedance phase,and tipper for each site.

Fig. 3 Depth slices at a 1800 m ASL, b 1400 m ASL, c 1000 m ASL, d 500 m
ASL and e 0 m ASL from the 3-D resistivity inversion at the peak area ofKusatsu–Shirane volcano. Seismic hypocenters are mapped as dots on each
map within 200 m tolerance in depth. The white, black, and gray dots denote
hypocenters for three periods before, during and after the 2014 unrest, namelyfrom August 2013 to February 2014, from March 2014 to May 2014, and from
June 2014 to February 2019. The white and black starts denote horizontal
locations of magnetizations and demagnetizations corresponding to volcanic
events from 1978 to 2012 after Takahashi & Fujii (2014). The red circles and open
curves show the crater rims of past eruptions. The rectangle in (b) denotes thefitting between the resistor R1 and the locations of historical eruptions. Major
fumaroles (F1 and F2) are shown by white closed curves

Results
The best-fit modeling results are shown as a set of depth
slices in Fig. 3. Here we describe the major features of themodel.
Figure 3a shows the resistivity distribution at 1800 m
above sea level (ASL), which is approximately 400 m
below the peak. The area to the east of the Yugama crater is covered with a widespread conductor (C1), whichis consistent with the high phase distribution at 70 to
17.5 Hz (Fig. 2). At 1400 m ASL (Fig. 3b), the conductor (C1) still exists but, the resistive anomalies R1 and R2
appear to the southeast and northwest of the area. The
C1 shows the semi-circular distribution, and inside thecircle, the resistor R1 appears in NE-SW directions andbreaking the circular distribution of C1 conductor in thesouthwestern directions. Interestingly, R1 follows the distributions of the historical phreatic eruption zones from1902 to 1942 (Terada 2018). At further depth, the circular distribution of C1 is disappearing towards 1000 m
ASL (Fig. 3c). The hypocenter of recent seismic activities (see Additional file 2: Text S2 for more details) andmagnetization/demagnetization sources (Takahashi andFujii 2014) are located at the north dge of the R1 resistor. From 1400 to 1000 m ASL, a columnar conductor(C2) is located at the western edge of the Yugama crater.
This conductive structure is a newly discovered feature,different from the previous research (e.g., Nurhasan et al.
2006) that we named it C2. There is a seismicity swarm
associated with the C2 to the east. At 500 m ASL andbelow (Figs. 3d, e), a deep conductor (C3) with a resistivity less than 1 Ωm appears in the northwestern part of
the area. The Yugama crater is located at the edge of theC3 conductor. The existence of the C3 conductor corresponds to the high phase distribution to the north of the
Yugama crater at 1.1 Hz (Fig. 2a).
We can also summarize the characteristics of theresistivity structure by sections (Fig. 4). Figure 4ashows the depth slice at 1400 m ASL, which is the sameas Fig. 3b with locations for the EW (Fig. 4b) and NS(Fig. 4c) sections. From Fig. 4b and c, the existence ofthe two conductors is evident. One is the shallow conductor (1–10 Ωm) covering the eastern half of the peakTseng et al. Earth, Planets and Space(2020) 72:161

Page 6 of 11

Fig. 4 a Depth slice at 1400 m ASL. The locations for W–E and S–N profiles for sections are shown by broken lines. b W–E and c S–N sections. The
seismic hypocenters are mapped as dots on each map within 200 m tolerance normal to the sections. The white, black, and gray dots denote
hypocenters before, during and after the 2014 unrest, respectively. The shallow seismicity (black dots) going through the clay cap (C1) is evident
during 2014 unrest. The white and black stars denote the same as in Fig. 3. The blue X denote the point pressure d The 3-D view of the three major
conductors. The blue surface shows the topography. The dots denote seismicity as in (a–c). source locations for explaining the radial inflation of theedifice during the 2014 unrest (Terada et al., 2015). Two white masks with irregular shape in a denote the fumarole zone at the north of Yugama
crater lake. The white line in c denotes the projection of fumarole zone F1 at 1400 m ASL

area (C1). It has a bell shape whose center is shifted
400 m to the east from the center of the Yugama crater (Y), as seen in the EW and NS resistivity sections.
The most active fumaroles in the northern edifice (F1and F2) are located at the inner rim of the C1 conductor; however, these fumaroles do not show any otherdeeper features. The other conductor (C3) is a deep one
located below 400 m ASL, which is beneath the Yugama
crater and extends to the northwest (see also Fig. 3e).
The Yugama crater is located at the edge of the C3. We
can also see the subvertical conductor C2 in Fig. 4b at
the western margin of the Yugama crater, which connects the deep (C3) and the shallow bell-shape (C1)
conductors.

The resistive block (R1) starts to appear clearly at
1400 m (Fig. 4a). The shape of the R1 is elongated in the
NE–SW direction to the east of the Yugama crater (Y),and the C1 conductor does not fully show an axially symmetrical structure. This feature means that the clay cap isthin along with the R1 resistor.
Figure 4b section shows that the R1 resistor is covered by the clay cap (C1). We also see the resistive block(R3) to the east at depth, which needs to exist in order
to explain the regional elongated phase tensor ellipse at
low frequencies around 0.3 Hz (Fig. 2). In the NS section(Fig. 4c), we can also see the shallow bell-shaped conductor (C1) above 1400 m ASL, and the deep conductor

Tseng et al. Earth, Planets and Space(2020) 72:161(C3) below 400 m ASL. The seismicity is located between
these deep and shallow conductors.
We tested the existence of C2 conductor, which connects C1 and C3 to the west of the Yugama crater(Additional file 2: Figure S17). The dotted block volume(ModC2 in Additional file 2: Figure S17) was replaced by
100 Ωm and resistivity contrast between C2, R1 and R2
are lost. The test model (ModC2) disconnects the deep
conductor and the shallow bell-shaped conductor. The
final RMS increased from 1.20 to 1.38. The phase tensors
at 4.4 Hz clearly show a significantly worse fit at the sites
at the southwest of the Yugama crater (Y) (Additional
file 2: Figure S18c). This change in the phase tensor plot ismore evident than the sounding curves (Additional file 2:
Figures S19 and S20), but the deviation of apparent resistivity is still distinguishable. Therefore, the existence of
C2 between the bell-shaped conductor (C1) and the deep
conductor (C3) was confirmed.
Next, we tested the existence of the C3 conductor
by replacing the resistivity in the dotted block volume(ModC3 in Additional file 2: Figure S17) with 50 Ωm.
The RMS increased from 1.20 to 1.34, and a significantly
worse fit is also seen in Additional file 2: Figure S18d at
1.1 Hz, where the calculated phase tensors from the test
model response show lower tan−1 (2 ) phases at thenorthwest of Yugama crater than those observed (Additional file 2: Figure S18a). The response sounding curves
represent a significant divergence between the best final
model and test model (Additional file 2: Figures S21,
S22). Thus the existence of C3 is confirmed.

Discussion
Clay–cap conductor (C1) and inner resistor (R1)

The distribution of the shallow conductor (C1) was also
imaged previously by Nurhasan et al. (2006) and Matsunaga et al. (2020) on the eastern slope of the volcano.
By reference to the drill hole data, they interpreted the
conductor as a layer containing clay (smectite), which
has low permeability as well as high electrical conductivity. In this paper, the distribution of the clay cap aroundthe Yugama crater was imaged in detail, and it was foundto be nearly bell shaped. These capping structures, thehydrothermal alteration products, are commonly foundin geothermal areas (e.g., Heise et al. 2008; Yamaya et al.
2013; Seki et al. 2015; Piña-Varas et al. 2018; Yoshimura
et al. 2018; Kanda et al. 2019). The bottom of the C1
smectite clay conductor is interpreted as the isotherm of
200 °C. The R1 resistor below C1 implies the smectite–
illite transformation above 200 °C (Ussher et al. 2000).
We can compare the distribution of C1 with seismicity. The hypocenters were relocated using the best fitting
seismic velocity model with the seismic monitoring data,
including three borehole seismograms surrounding thePage 7 of 11

Yugama crater within 500 m distances (See Additional
file 2: Text S2–S4for details). As seen in the resistivity
sections (Fig. 4b, c), most of the hypocenters are locatedbelow C1 conductor. This feature implies that the rising
fluids increase the pore pressures and induce the microseismicity. However, the flow is inhibited by the low-permeability clay cap below the Yugama crater (Fig. 4b, c).
Few earthquakes occurred inside the clay cap (C1) during
the 2014 unrest. This may imply that fluids invaded into
the clay cap during the 2014 unrest and supplied volcanic
fluids into the crater lake Yugama.
Takahashi and Fujii (2014) analyzed the temporalvariations of total magnetic intensity data from 1982 to2012 and inverted the corresponding magnetization anddemagnetization dipole locations to explain the temporal variation of temperature changes underneath the volcano. The magnetization and demagnetization source
locations are also plotted in Fig. 4a, b by black and whitestars, respectively. In the sections, the black and whitestars are located near the R1 and C1 boundaries. This feature indicates that in the active period of the volcano, thehigh-temperature fluid reaches the bottom of the lowpermeability clay cap (C1) and may partly intrude into
the clay cap during the past unrests.
Nakano et al. (2003) modeled the long-period volcanic
tremors during volcanic unrest (during the 1990s) at theKusatsu–Shirane volcano and inferred a horizontal crack
resonator which is located 300 m below Yugama crater.
This inferred resonator is also located consistently withthe demagnetization/magnetization sources, just below
the clay cap at the top of R1 resistor.
During the 2014 unrest, we observed the radial inflation of the edifice by using three borehole tilt-meterssurrounding the Yugama crater, simultaneously with theonset of increased seismicity (Terada et al. 2015). The
location of the causative point source of inflation is indicated by the green “X” in Fig. 4a–c (Terada et al. 2015).
The inflation source is located just below the clay cap at
the northeastern edge of the Yugama crater. The inflation
source will be interpreted by a supply of high-pressure
volcanic hydrothermal fluid. Tsukamoto et al. (2018) also
show a similar high-pressure source at Kirishima volcano, which lead the precursor of a phreatic eruption.
As seen in Fig. 3b, the resistive area of R1 starts to
spread horizontally inside the C1 conductor as theradius of C1 increases with depth until the C1 becomesless significant at 1000 m ASL. The resistor implies ahigher temperature zone where the smectite will become
unstable and be converted to illite, which is electricallymore resistive (Eberl 1978; Inoue et al. 1992). R1 is not
circularly distributed symmetrically with Yugama but
is elongated in the NE to SW direction. R1 follows the
distribution of the historical phreatic eruption craters

Tseng et al. Earth, Planets and Space(2020) 72:161

of 1902, 1932, and 1942, whose shapes are sketched in
Figs. 3b and 4a (Terada, 2018). R1 breaks the circularsymmetry of C1 conductor. The past phreatic eruptions
may have destroyed the symmetric clay cap structure.
Brine conductor (C3)

A significant conductor (C3) with a resistivity less than
1 Ωm appears at 400 m ASL below the Yugama craterand extends to the northwest (Fig. 3d, e and Fig. 4b, c).
The cutoff depth of microseismicity coincides with thebottom of C2, and the top of C3 conductor corresponds
to brittle–ductile transition at around 400 °C. This conductor will be interpreted as brine above 400 °C, which
will be separated from the magmatic fluids (Fournier
1999; Afanasyev et al. 2018). The top of the C3 fluids will
be self-sealed by silica (Saishu et al. 2014) or by halite(Afanasyev et al. 2018), which can last for several hundred thousand years. The seal on top of C3 will separate
the hydrostatic domain and underlying the lithostatic
domain. The episodic supply of magmatic fluids or gas
to C3 from below will increase the internal pressure of
C3. The part of the seal will break and released fluids will
increase the pore pressure in the overlying hydrostatic
domain and will induce microseismicity.
Such conductors below the seismogenic zones arecommonly found in volcanic regions and also under theseismically active regions (e.g., Ogawa et al. 2001, 2014;
Wannamaker et al. 2009; Ichihara et al. 2014, 2016).
Vertical conductor (C2) connecting the brine (C3)to the surface (C1)

We also found a vertical conductor (C2) that connects
the deep brine conductor (C3) to the clay cap rock conductor (C1). It can be regarded as a part of the deep
extension of C1 conductor, but at 1000 m ASL, C1 conductor almost disappears, but C2 exists. It is also connected to the deeper conductor C3. Thus it is difficult tointerpret C2 as the extension of smectite clay cap (C1),which must be lsss than 200 °C. This vertical conductor
appears at around 1000–1400 m ASL (Fig. 3b, c), which
is significantly deeper than C1 clay cap, and is associatedwith the microseismic hypocenters.
The microseismic activities can be interpreted as thefluid‐induced swarm earthquakes when the overpressured fluid was invading into the fracture in the structureand diffusing, which is a typical feature in the volcanichydrothermal system (Yukutake et al. 2011). The fracturein this structure can provide the fluid pathway and can
reduce the bulk resistivity, although the rock porosity isdecreased (e.g., Vidal et al. 2018). The vertical conductor C2 is interpreted as a conduit pathway of fluids fromthe deep brine reservoir to the surface. Similar verticalconductors connecting the deep reservoir to the surface

Page 8 of 11

were also suggested for other volcanoes (Seki et al. 2015;
Tsukamoto et al. 2018). The invaded and diffused fluid
made the structure seismogenic. When the volume of
upward-invaded fluid was significant, or the pressure
was relatively high, it induced the detectable microseismic activity (in this study, the threshold magnitude was− 3.0).
Conversely, when the fluid diffused in the structurewith relatively lower volume and pressure, the microseismic activity could be too small to be detected. Thisdifference could explain why the conductor C2 is not
fully overlapping with the swarm plot of microseismic
hypocenters in Fig. 4b. It also explains the occurrence ofshallow seismic activity during 2014 unrest, because thevolume of upward-released fluid increased then the fluid‐
induced microseismic activity became detectable.
Figure 4d represents the 3-D cartoon of the hydrothermal system of Kusatsu–Shirane volcano. The clay cap
C1 located around the Yugama crater lake appeared in abell shape and overlying the microseismic hypocenters.
The microseismic hypocenters that indicated the major
fluid pathway invaded upward. The fluid-filled fractured
structure C2 is underlain by the clay cap C1. The brinereservoir C3 provides the filling hydrothermal fluid.
The fluid released from this deep reservoir flow upward
through the pathway conduit C2 to the surface supported
the whole hydrothermal system of the Kusatsu–Shirane
volcano.
Figure 5 demonstrates the phase transition of thereleased fluid using the diagram of Hutnak et al. (2009).
We utilized the depth of microseismic swarm to speculate the bottom of the fluid pathway in the brittle–ductile transition zone (1400 m depth, 400 °C). We alsoutilized the bottom of C1 clay cap (400 m depth, 200 °C)(Fig. 5a). The phase transition between these two points
is represented in the H­ 2O–CO2 coexisting phase diagrams (Fig. 5b). We think that C2 is an established fluid
pathway from the brine. In contrast, the seismicity above
the bine may be induced by breaking the silica cap dueto the over-pressure caused by an episodic large flux ofmagmatic gas, such as ­CO2. This interpretation follows
the porphyry model of Blundy et al. (2015). When thesilica cap is broken, the pressure will be decreased fromlithostatic to hydrostatic and brine will be first vaporized(Weatherley and Henley, 2013). In a hydrostatic domain,
the phase of fluid will transit from pure vapor to liquid–gas-mixed fluid during upward migration. The sudden
increase in liquid–gas-mixed fluid underlain by the clay
cap C1 can lead the phreatic eruption. Therefore, theresistivity of the fluid pathway conduit C2 and its surrounding space should be monitored for predicting thephreatic eruption in the future.

Tseng et al. Earth, Planets and Space(2020) 72:161

Page 9 of 11

demagnetization sources of the past unrests and theinflation point source of the 2014 unrest are locatedjust below the clay cap (C1). This evidence supportsthat C1 works as the cap structure of the hydrothermal
system.
The resistor (R1) underneath the clay cap (C1) runs in
NE–SW directions and destroys the radial symmetry of
the clay cap, which suggests that the past phreatic eruptions may have broken part of the clay cap. The fluid pathfrom the C3 reservoir to the surface was imaged by a vertical conductor (C2). In contrast, the episodic magmaticgas flux may break the silica cap over the brine reservoir
and vapor and the magmatic gas will be flushed into thehydrostatic domain and will induce seismicity during
migration to the bottom of the clay cap.
The resistivity model provides the comprehension of
the dynamics of the magmatic hydrothermal systembeneath Kusatsu–Shirane volcano.

Supplementary information
Supplementary information accompanies this paper at https​://doi.
org/10.1186/s4062​3-020-01283​-2.
Additional file 1: Table S1. The site locations, altitude and magnetic fieldsite which was used for the calculations of the impedances.
Additional file 2: Text S1. Definition of the error floor. Text S2. Information of microseismicity data. Text S3. Observed and calculated data with
its location on the map (Figure S2–S16). Text S4. Sensitivity tests of the
conductors C2 and C3 (Figure S17–S22). Figure S23. Deeper conductive
feature beneath Kusatsu–Shirane volcano. Figure S24–S26. The seismic
hypocenters before 2014 unrest mapping with a Depth slice at 1400m
ASL b W–E and c S–N sections.

Fig. 5 a Simplified W–E resistivity section in Fig. 4b. The rectangles
on the terrain denote the observation sites. The white and blackstars and the green “X” symbols are the same as in Fig. 4. The swarm
dots with bigger and with smaller circles denote the hypocenters of
microseismic events during 2014 unrest, and out of this active period,
respectively. b Pressure–temperature projection of hydrothermal fluid­ O2 (after Hutnak et al. 2009). Two curves
phase diagram for ­H2O and C
which end with solid circles denote the liquid–vapor coexistence
curves of pure ­CO2 and pure H­ 2O with their critical points (CP). C2
may be an established fluid path from the supercritical fluid. In
contrast, if episodic large gas flux breaks the silica cap, breakage ofthe silica cap may flush the fluid upward in the hydrostatic domain

Conclusions
The magmatic hydrothermal system of Kusatsu–Shirane volcano was imaged by the three-dimensional
magnetotelluric modeling. The final model showed the
low-permeability clay cap as near-surface conductor(C1) and the brine reservoir as a deep conductor (C3) ata depth of 1.4 km from the surface. The magnetization/

Abbreviations
MT: Magnetotelluric; AMT: Audio-magnetotelluric; RMS: Root-mean-square;
ASL: Above sea level.
Acknowledgements
We thank Grant Caldwell for fruitful discussions for modeling. The constructive comments from Koki Aizawa and an anonymous reviewer improved themanuscript.
Authors’ contributions
YO designed the magnetotelluric field survey and supervised the research.
KHT performed the 3-D inversion. KHT and YO made interpretations. N and
SBT made preliminary magnetotelluric modelings. YO, N, SBT, NU, YH carriedout magnetotelluric data acquisitions. AT provided seismic data, YU and WK
contributed in the 3-D inversion. All authors read and approved the final
manuscript.
Funding
This study was supported by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan, under its Earthquake and Volcano
Hazards Observation and Research Program. Iw was also supported by Grantin-Aid for Scientific Research KAKENHI 20H01992.
Availability of data and materials
The data is available upon request to the corresponding author.
Competing interests
The authors declare that they have no competing interests.

Tseng et al. Earth, Planets and Space(2020) 72:161

Author details
1
Department of Earth and Planetary Sciences, Tokyo Institute of Technology,
Tokyo, Japan. 2 Volcanic Fluid Research Center, Tokyo Institute of Technology, Tokyo, Japan. 3 Present Address: Physics Department, Bandung Instituteof Technology, Bandung, Indonesia. 4 Present Address: Boğaziçi University,
Kandilli Obs. & E.R.I., Çengelköy, İstanbul, Turkey. 5 Present Address: Hydrographic and Oceanographic Department, Japan Coast Guard, Tokyo, Japan.
6
Present Address: Earthquake Research Institute, The University of Tokyo,
Tokyo, Japan.
Received: 15 February 2020 Accepted: 28 September 2020

References
Aizawa K, Ogawa Y, Ishido T (2009) Groundwater flow and hydrothermal
systems within volcanic edifices: Delineation by electric self-potential and magnetotellurics. J Geophys Res Solid Earth. https​://doi.
org/10.1029/2008J​B0059​10
Afanasyev A, Blundy J, Melnik O, Sparks S (2018) Formation of magmatic brinelenses via focussed fluid-flow beneath volcanoes. Earth Planet Sci Lett
486:119–128. https​://doi.org/10.1016/j.epsl.2018.01.013
Blundy J, Mavrogenes J, Tattitch B, Sparks S, Gilmer A (2015) Generation ofporphyry copper deposits by gas–brine reaction in volcanic arcs. Nature
Geosci 8:235–240. https​://doi.org/10.1038/ngeo2​351
Caldwell TG, Bibby HM, Brown C (2004) The magnetotelluric phase tensor. Geophys J Int 158:457–469. https​://doi.org/10.1111/j.1365246X.2004.02281​.x
Eberl D (1978) Reaction series for dioctahedral smectites. Clays Clay Miner
26:327–340. https​://doi.org/10.1346/CCMN.1978.02605​03
Fournier RO (1999) Hydrothermal processes related to movement of fluid fromplastic into brittle rock in the magmatic–epithermal environment. Econ
Geol 94:1193–1211. https​://doi.org/10.2113/gseco​ngeo.94.8.1193
Heise W, Caldwell TG, Bibby HM, Bannister SC (2008) Three-dimensionalmodelling of magnetotelluric data from the Rotokawa geothermal field,
Taupo Volcanic Zone, New Zealand. Geophys J Int 173:740–750. https​://doi.org/10.1111/j.1365-246X.2008.03737​.x
Herd DG (1986) The 1985 Ruiz Volcano disaster. Eos Trans Am Geophys Union
67:457. https​://doi.org/10.1029/EO067​i019p​00457​-03
Hutnak M, Hurwitz S, Ingebritsen SE, Hsieh PA (2009) Numerical models of caldera deformation: effects of multiphase and multicomponent hydrothermal fluid flow. J Geophys Res 114:B04411. https​://doi.org/10.1029/2008J​
B0061​51
Ichihara H, Mogi T, Tanimoto K, Yamaya Y, Hashimoto T, Uyeshima M, Ogawa
Y (2016) Crustal structure and fluid distribution beneath the southernpart of the Hidaka collision zone revealed by 3-D electrical resistivity modeling. Geochem Geophys Geosyst 17:1480–1491. https​://doi.
org/10.1002/2015G​C0062​22
Ichihara H, Sakanaka S, Mishina M, Uyeshima M, Nishitani T, Ogawa Y, Yamaya Y,
Mogi T, Amita K, Miura T (2014) A 3-D electrical resistivity model beneaththe focal zone of the 2008 Iwate-Miyagi Nairiku earthquake (M 7.2). Earth
Planets Space 66:50. https​://doi.org/10.1186/1880-5981-66-50
Inoue A, Utada M, Wakita K (1992) Smectite-to-illite conversion innatural hydrothermal systems. Appl Clay Sci 7:131–145. https​://doi.
org/10.1016/0169-1317(92)90035​-L
Kanda W, Utsugi M, Takakura S, Inoue H (2019) Hydrothermal system of theactive crater of Aso volcano (Japan) inferred from a three-dimensional
resistivity structure model. Earth Planets Space 71:37. https​://doi.
org/10.1186/s4062​3-019-1017-7
Kaneko T, Maeno F, Nakada S (2016) 2014 Mount Ontake eruption: characteristics of the phreatic eruption as inferred from aerial observations. Earth
Planets Space 68:72. https​://doi.org/10.1186/s4062​3-016-0452-y
Kobayashi T, Morishita Y, Munekane H (2018) First detection of precursory
ground inflation of a small phreatic eruption by InSAR. Earth Planet Sci
Lett 491:244–254. https​://doi.org/10.1016/j.epsl.2018.03.041
Martí A, Queralt P, Marcuello A, Ledo J, Rodríguez-Escudero E, Martínez-Díaz
JJ, Campanyà J, Meqbel N (2020) Magnetotelluric characterizationof the Alhama de Murcia Fault (Eastern Betics, Spain) and study ofmagnetotelluric interstation impedance inversion. Earth Planets Space
72:16. https​://doi.org/10.1186/s4062​3-020-1143-2

Page 10 of 11

Matsunaga Y, Kanda W, Takakura S, Koyama T, Saito Z, Seki K, Suzuki A,
Kishita T, Kinoshita Y, Ogawa Y (2020) Magmatic hydrothermal system
inferred from the resistivity structure of Kusatsu–Shirane Volcano. J
Volcanol Geotherm Res 390:106742. https​://doi.org/10.1016/j.jvolg​eores​.2019.10674​2
Nakano M, Kumagai H, Chouet BA (2003) Source mechanism of long-period
events at Kusatsu–Shirane Volcano, Japan, inferred from waveform
inversion of the effective excitation functions. J Volcanol Geotherm Res
122:149–164. https​://doi.org/10.1016/S0377​-0273(02)00499​-7
Nurhasan, Ogawa Y, Ujihara N, Tank SB, Honkura Y, Onizawa S, Mori T, Makino
M (2006) Two electrical conductors beneath Kusatsu–Shirane volcano,Japan, imaged by audiomagnetotellurics, and their implications for thehydrothermal system. Earth Planets Space 58:1053–1059. https​://doi.
org/10.1186/BF033​52610​
Ogawa Y, Mishina M, Goto T, Satoh H, Oshiman N, Kasaya T, Takahashi
Y, Nishitani T, Sakanaka S, Uyeshima M, Takahashi Y, Honkura Y,
Matsushima M (2001) Magnetotelluric imaging of fluids in intraplate
earthquake zones, NE Japan Back Arc. Geophys Res Lett 28:3741–3744.
https​://doi.org/10.1029/2001G​L0132​69
Ogawa Y, Ichiki M, Kanda W, Mishina M, Asamori K (2014) Three-dimensional
magnetotelluric imaging of crustal fluids and seismicity around Naruko
volcano NE Japan. Earth Planets Space 66:158. https​://doi.org/10.1186/
s4062​3-014-0158-y
Ohba T, Yaguchi M, Nishino K, Numanami N, Tsunogai U, Ito M, Shingubara
R (2019) Time variation in the chemical and isotopic composition of
fumarolic gasses at Kusatsu–Shirane Volcano, Japan. Front Earth Sci.
https​://doi.org/10.3389/feart​.2019.00249​
Piña-Varas P, Ledo J, Queralt P, Marcuello A, Perez N (2018) On the detectability of Teide volcano magma chambers (Tenerife, Canary Islands)
with magnetotelluric data. Earth Planets Space 70:14. https​://doi.
org/10.1186/s4062​3-018-0783-y
Saishu H, Okamoto A, Tsuchiya N (2014) The significance of silica precipitation on the formation of the permeable–impermeable boundary
within Earth’s crust. Terra Nova 26:253–259. https​://doi.org/10.1111/
ter.12093​
Sano Y, Hirabayashi J-I, Oba T, Gamo T (1994) Carbon and helium isotopic ratios
at Kusatsu–Shirane Volcano, Japan. Appl Geochem 9:371–377. https​://doi.
org/10.1016/0883-2927(94)90059​-0
Seki K, Kanda W, Ogawa Y, Tanbo T, Kobayashi T, Hino Y, Hase H (2015) Imaging the hydrothermal system beneath the Jigokudani valley, Tateyama
volcano, Japan: implications for structures controlling repeated phreatic
eruptions from an audio-frequency magnetotelluric survey. Earth Planets
Space 67:6. https​://doi.org/10.1186/s4062​3-014-0169-8
Stix J, de Moor JM (2018) Understanding and forecasting phreatic eruptions
driven by magmatic degassing. Earth Planets Space 70:83. https​://doi.
org/10.1186/s4062​3-018-0855-z
Takahashi K, Fujii I (2014) Long-term thermal activity revealed by magnetic
measurements at Kusatsu–Shirane volcano, Japan. J Volcanol Geotherm
Res 285:180–194. https​://doi.org/10.1016/j.jvolg​eores​.2014.08.014
Tazieff H (1989) Mechanisms of the Nyos carbon dioxide disaster and of socalled phreatic steam eruptions. J Volcanol Geotherm Res 39:109–116.
https​://doi.org/10.1016/0377-0273(89)90051​-6
Terada A (2018) Kusatsu–Shirane volcano as a site of phreatic eruptions. J Geol
Soc Japan 124:251–270. https​://doi.org/10.5575/geoso​c.2017.0060
Terada A, Morita Y, Hashimoto T, Mori T, Ohba T, Yaguchi M, Kanda W (2018)
Water sampling using a drone at Yugama crater lake, Kusatsu–Shirane
volcano, Japan. Earth Planets Space 70:64. https​://doi.org/10.1186/s4062​
3-018-0835-3
Terada A, Ohkura T, Kanda W, Ogawa Y (2015) Ground deformation causedby an accumulation of hydrothermal water beneath hot crater lake at
Kusatsu– Shirane volcano. In: Abstracts of Japan geoscience union meeting 2015, Makuhari, Chiba, Japan, 24-28 May 2015.
Tsukamoto K, Aizawa K, Chiba K, Kanda W, Uyeshima M, Koyama T, Utsugi M,
Seki K, Kishita T (2018) Three-dimensional resistivity structure of Iwo-Yama
Volcano, Kirishima Volcanic Complex, Japan: relationship to shallow
seismicity, surface uplift, and a small phreatic eruption. Geophys Res Lett.
https​://doi.org/10.1029/2018G​L0802​02
Ussher G, Harvey C, Johnstone R, Anderson E (2000) Understanding theresistivities observed in geothermal systems. In: Proceedings worldgeothermal congress, Kyushu-Tohoku, Japan, 28 May–10 June 2000, pp
1915–1920

Tseng et al. Earth, Planets and Space(2020) 72:161

Usui Y (2015) 3-D inversion of magnetotelluric data using unstructured tetrahedral elements: applicability to data affected by topography. Geophys J
Int 202:828–849. https​://doi.org/10.1093/gji/ggv18​6
Usui Y, Ogawa Y, Aizawa K, Kanda W, Hashimoto T, Koyama T, Yamaya Y,
Kagiyama T (2017) Three-dimensional resistivity structure of Asama
Volcano revealed by data-space magnetotelluric inversion using unstructured tetrahedral elements. Geophys J Int 208:1359–1372. https​://doi.
org/10.1093/gji/ggw45​9Vidal J, Patrier P, Genter A, Beaufort D, Dezayes C, Glaas C, Lerouge C, Sanjuan
B (2018) Clay minerals related to the circulation of geothermal fluids
in boreholes at Rittershoffen (Alsace, France). J Volcanol Geotherm Res
349:192–204. https​://doi.org/10.1016/j.jvolg​eores​.2017.10.019
Wannamaker PE, Caldwell TG, Jiracek GR, Maris V, Hill GJ, Ogawa Y, Bibby HM,
Bennie SL, Heise W (2009) Fluid and deformation regime of an advancing
subduction system at Marlborough, New Zealand. Nature 460:733–736.
https​://doi.org/10.1038/natur​e0820​4
Weatherley D, Henley R (2013) Flash vaporization during earthquakes
evidenced by gold deposits. Nature Geosci 6:294–298. https​://doi.
org/10.1038/ngeo1​759

Page 11 of 11

Yamaya Y, Alanis PKB, Takeuchi A, Cordon JM, Mogi T, Hashimoto T, Sasai Y,
Nagao T (2013) A large hydrothermal reservoir beneath Taal Volcano (Philippines) revealed by magnetotelluric resistivity survey: 2D resistivity modeling. Bull Volcanol 75:729. https​://doi.org/10.1007/s0044​5-013-0729-y
Yoshimura R, Ogawa Y, Yukutake Y, Kanda W, Komori S, Hase H, Goto T, Honda
R, Harada M, Yamazaki T, Kamo M, Kawasaki S, Higa T, Suzuki T, Yasuda Y,
Tani M, Usui Y (2018) Resistivity characterization of Hakone volcano, Central Japan, by three-dimensional magnetotelluric inversion. Earth Planets
Space 70:66. https​://doi.org/10.1186/s4062​3-018-0848-y
Yukutake Y, Ito H, Honda R, Harada M, Tanada T, Yoshida A (2011) Fluid-induced
swarm earthquake sequence revealed by precisely determined hypocenters and focal mechanisms in the 2009 activity at Hakone volcano,Japan. J Geophys Res 116:B04308. https​://doi.org/10.1029/2010J​B0080​3

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.