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

The 2018 phreatic eruption at Mt. Motoshirane of Kusatsu‒Shirane
volcano, Japan: Eruption and intrusion of hydrothermal fluid observed
by a borehole tiltmeter network

Authors

Akihiko Terada, Wataru Kanda, Yasuo Ogawa, Taishi Yamada, Mare
Yamamoto, Takahiro Ohkura, Hiroshi Aoyama, Tomoki Tsutsui,
Shin ya Onizawa

Citation

Earth, Planets and Space, Vol. 73, ,Pub. date

2021, 7

DOI

https://doi.org/10.1186/s40623-021-01475-4

Creative Commons

The information is in the article.

Powered by T2R2 (Science Tokyo Research Repository)

(2021) 73:157
Terada et al. Earth, Planets and Space
https://doi.org/10.1186/s40623-021-01475-4

Open Access

FULL PAPER

The 2018 phreatic eruption at Mt.
Motoshirane of Kusatsu–Shirane volcano, Japan:
eruption and intrusion of hydrothermal fluid
observed by a borehole tiltmeter network
Akihiko Terada1* , Wataru Kanda1 , Yasuo Ogawa1 , Taishi Yamada2 , Mare Yamamoto3, Takahiro Ohkura4,
Hiroshi Aoyama5 , Tomoki Tsutsui2 and Shin’ya Onizawa6

Abstract
We estimate the mass and energy budgets for the 2018 phreatic eruption of Mt. Motoshirane on Kusatsu–Shirane
volcano, Japan, based on data obtained from a network of eight tiltmeters and weather radar echoes. The tilt records
can be explained by a subvertical crack model. Small craters that were formed by previous eruptions are aligned
WNW–ESE, which is consistent with the strike of the crack modeled in this study. The direction of maximum compressive stress in this region is horizontal and oriented WNW–ESE, allowing fluid to intrude from depth through acrack with this orientation. Based on the crack model, hypocenter distribution, and MT resistivity structure, we inferthat fluid from a hydrothermal reservoir at a depth of 2 km below Kusatsu–Shirane volcano has repeatedly ascended
through a pre-existing subvertical crack. The inflation and deflation volumes during the 2018 eruption are estimatedto have been 5.1 × ­105 and 3.6 × ­105 ­m3, respectively, meaning that 1.5 × ­105 ­m3 of expanded volume formed
underground. The total heat associated with the expanded volume is estimated to have been ≥ ­1014 J, similar to or
exceeding the annual heat released from Yugama Crater Lake of Mt. Shirane and that from the largest eruption during
the past 130 year. Although the ejecta mass of the 2018 phreatic eruption was small, the eruption at Mt. Motoshirane
was not negligible in terms of the energy budget of Kusatsu–Shirane volcano. A water mass of 0.1–2.0 × ­107 kg wasdischarged as a volcanic cloud, based on weather radar echoes, which is smaller than the mass associated with the
deflation. We suggest that underground water acted as a buffer against the sudden intrusion of hydrothermal fluids,
absorbing some of the fluid that ascended through the crack.
Keywords: Phreatic eruption, Tiltmeter, Kusatsu–Shirane volcano, Hydrothermal system, Tensile crack, Brittle–ductile
transition zone
Introduction
In 2018, a phreatic eruption occurred on the Kagamiikekita and Kagamiike pyroclastic cone (PC) of the Mt.
Motoshirane Pyroclastic Cones Group (MPCG) on thesouthern part of Kusatsu–Shirane volcano, Japan (Ogawa*Correspondence: terada@ksvo.titech.ac.jp
1
Volcanic Fluid Research Center, School of Science, Tokyo Instituteof Technology, 2‑12‑1 Ookayama, Meguro, Tokyo 152‑8551, Japan
Full list of author information is available at the end of the article

et al. 2018; Yamada et al. 2021). The total mass of the fall
deposit is estimated to have been ~ 3.0 × ­104 tonne (Kametani et al. 2021). The eruption products were mainly clay
minerals derived from a high-temperature acid alteration
zone beneath the volcano; no juvenile material has beenfound (Yaguchi et al. 2019).
Until several minutes before the occurrence of theeruption, no significant precursors (e.g., increased seismicity, seismic velocity change, unusual ground deformation and unusual thermal activity) (Barberi et al.

© The Author(s) 2021. Open Access 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 theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother 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 the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/.

Terada et al. Earth, Planets and Space(2021) 73:157

1992) were detected around the MPCG (Yamamoto
et al. 2018). However, the Shirane Pyroclastic Cone
Group (SPCG), located 1.5 km north of the eruption
site, exhibited an increase in seismicity, accompanied
by a shallow inflation and changes in water chemistryof Yugama Crater Lake (YCL) and nearby fumarolicgas from 2014 to 2016 (Ohba et al. 2019a). Moreover,
3 months after the 2018 MPCG eruption, a shallow
inflation and a marked increase in seismicity occurredat the SPCG. Therefore, Kusatsu–Shirane volcano
can be categorized as a “persistently restless volcano”(Roman et al. 2019) or a “volcano with slow unrest”(Stix 2018). Small craters formed by past phreatic eruptions occur throughout the northern part of the MPCG,
such as the Kagamiike-kita and Kagamiike PC (Ishizaki et al. 2020). Thus, understanding the relationshipbetween the hydrothermal systems beneath the MPCG

Page 2 of 17

and SPCG can provide insights into the risk and processes associated with a small phreatic eruption.
At the onset of the 2018 MPCG eruption, eight tiltmeters (including six borehole-type; Fig. 1), including KSS
located 1.2 km from the main crater (MC; Figs. 1 and2), were in operation at Kusatsu–Shirane volcano. Sucha tiltmeter network with high temporal resolution anddense spatial distribution provides details of progressive ground deformation during an eruption (Ueda et al.
2013; Aloisi et al. 2020; Zobin et al. 2020). Moreover,the volcanic ash cloud was monitored by weather radar(Meteorological Research Institute 2018). The volcanic
plume track allows the heat discharge rate of an ash
cloud to be estimated, which can then be converted to amass flux of water (Briggs 1969; Kagiyama 1981; Terada
and Sudo 2012; Narita et al. 2019), although surveillancecameras and photographs taken from the ground did not

Fig. 1 Location and topography of Kusatsu–Shirane volcano. a Location map of Kusatsu–Shirane volcano in Japan (star). Blue and red circles
represent the locations of the Nagano weather radar site and an aerological observatories using a balloon and radiosonde, respectively (bothoperated by the Japan Meteorological Agency). b Topographic map of Kusatsu–Shirane volcano. Solid lines are topographic contours at intervalsof 30 m. Red squares and circles indicate the locations of tiltmeters installed at the bottom of boreholes and at a depth of 1–2 m from theground surface, respectively (Table 1). Yellow squares (MC, WC, and SC) indicate the locations of craters that formed during the 2018 eruption.
PC = pyroclastic cone; MPCG = Motoshirane Pyroclastic Cone Group; SPCG = Shirane Pyroclastic Cone Group; YCL = Yugama Crater Lake; MC = main
crater; WC = west crater; SC = south crater

Terada et al. Earth, Planets and Space(2021) 73:157

Page 3 of 17

Fig. 2 a Aerial photograph of Kusatsu–Shirane volcano (taken on 5 November 2015). PC = pyroclastic cone; MPCG = Motoshirane Pyroclastic
Cone Group; SPCG = Shirane Pyroclastic Cone Group; YCL = Yugama Crater Lake. b Eruption cloud and ballistic block falls hitting the ski slope,(photograph taken by the Kusatsu Tourism Cooperation at 10:02 JST on 23 January 2018). c Weak volcanic cloud rising from the Kagamiike-kita
PC at about 10:20 JST on 23 January 2018. d Aerial photograph of the new main crater (MC) at the Kagamiike-kita PC and the south crater (SC) at
the Kagamiike PC formed during the 2018 eruption (photograph taken on 25 June 2018). The west crater (WC) is outside the photographconstrain the height of the ash cloud at the climax of theeruption.
In this paper, we present a simple model of inflation–
deflation during the 2018 MPCG eruption based on thetilt changes recorded at eight tiltmeter stations. In addition, we estimate the mass of water discharged as an
ash cloud. On the basis of these estimates, we discuss
the mass and energy budgets during the eruption. At
Kusatsu–Shirane volcano, the resistivity structure has
been investigated by Magnetotelluric (MT) sounding
methods (Nurhasan et al. 2006; Matsunaga et al. 2020;
Tseng et al. 2020). Such a model of the undergroundstructure provides a constraint of a fluid pathway beneath
Kusatsu–Shirane volcano. On the basis of the underground structure, the observations during the eruption in2018 at the MPCG and the progression of the unrest of
the SPCG since 2014, we propose a schematic model ofthe 2018 MPCG eruption.

Kusatsu–Shirane volcano and the 2018 eruption
Kusatsu–Shirane volcano consists of three pyroclastic
cones (Figs. 1 and 2a), which are the SPCG, Ainomine
Pyroclastic Cone (APC; undated but probably Holocene in age), and MPCG that are aligned N–S over a
distance of 4 km (Ishizaki et al. 2020). Kusatsu–Shirane
volcano exhibits a persistent heat release of 110 MW(Terada 2018), which is transferred by a hyper-acidic,
hot crater lake (locally referred to as YCL (Yugama
Crater Lake); Figs. 1 and 2a; Ohba et al. 1994; Teradaand Hashimoto 2017), nearby fumaroles around the
SPCG (Ohba et al. 2019a), and hot-spring water emissions on the flank of the volcano. Phreatic eruptions
have repeatedly occurred around the SPCG and, asa result, the ground surface of the SPCG is covered
with whitish tephra consisting of clay and is unvegetated. The eruption series in 1937–1939 was the largest eruption in the last 130 yr, and produced 5 × ­109 kg

Terada et al. Earth, Planets and Space(2021) 73:157

Page 4 of 17

of tephra (Minakami 1939). The latest eruption of the
SPCG occurred at YCL in 1982–1983. Recently, microearthquake swarms accompanied by ground deformation and changes in the chemistry of YCL and nearby
fumaroles were detected in 1991–1992 and 2014–2016(Ohba et al. 2008, 2019b).
The MPCG was mostly covered with vegetation prior
to the 2018 eruption (Figs. 1 and 2a). No eruption has
been recorded in historic times, whereas continuousseismic activity has been detected around the APC
and MPCG (Mori et al. 2006). On 23 January 2018, aphreatic eruption occurred at the Kagamiike-kita PC(slightly younger than 1500 cal. year BP) and the Kagamiike PC (ca. 4800 cal. yr BP), which are in thenorthern part of the MPCG. The onset of eruption
was preceded by ~ 2 min of volcanic tremor located
1 km north of the MC (Fig. 2) at a depth of 0.5–1.0 km
from the surface (Yamada et al. 2021). The tremor was
accompanied by rapid changes in tilt.
During the 2018 event, ballistic blocks fell up to
0.5 km from the MC and resulted in one death. Volcanic ash released during the eruption was found up
to 25 km ENE of the volcano. The maximum heightof the volcanic cloud was estimated to be 5500 m asl
by weather radar (Meteorological Research Institute2018), corresponding to 3400 m in height relative to thealtitude of the MC (Fig. 2). Eyewitness reports and seismic records suggest that most of the ash emission had
finished within ~ 10 min. A photograph taken at around
10:20–10:30 (JST) from Kusatsu town (Fig. 1) showeda whitish plume and black tephra deposit aroundthe Kagamiike-kita PC, but no ash cloud (Fig. 2c). During the eruption, the MC, the west crater (WC), andother smaller craters were formed on the Kagamiikekita PC (Fig. 2d), which are aligned on an azimuth of
112°. In addition, the south crater (SC) was formed on
the Kagamiike-kita PC (Fig. 2c) (Kametani et al. 2021).

Tilt records during the 2018 eruption
The Kusatsu–Shirane Volcano Observatory, Tokyo Institute of Technology (KSVO) and the National ResearchInstitute for Earth Science and Disaster Resilience(NIED) operate force-balanced-type tiltmeters (Satoet al. 1980) at the present study site, installed at the bottom of boreholes at depths of 50–200 m (Table 1). The
output sensitivity is 20 mV/µrad at KSE and KSW, and5 mV/µrad at KSS, N.KSHV, N.KSYV, and N.KSNV (Tanada et al. 2017), which produces a high tilt resolution of­10–8–10–9 rad. Applied Geomechanics 701-2A dual-axis
tiltmeters with an output sensitivity of 1 mV/µrad were
deployed at depths of 2 and 1 m at JIE and YNW, respectively. The 701-2A contains electrolytic level sensors that
produce changes in resistance in response to a rotationof the sensor. Tilt records are sampled at 20 Hz (NIED)or 200 Hz (KSVO). In this study, acausal low-pass filterswith cut-off frequencies of 0.034 Hz and 0.021 Hz wereapplied to the tilt records of JTS-33 and 701-2A (Table 1),
respectively. No data were obtained for the N–S component at JIE owing to a mechanical problem. A tiltmeter at
V.KSAO, 1.6 km northeast of MC, operated by the Japan
Meteorological Agency was not operational when the2018 eruption occurred but did record the 2011 event.
Ground deformation detected by the tiltmeter network

We detected a simple and significant variation in thetilt records from 10:00 (all times are in Japan Standard Time, JST) on 23 January 2018 (Fig. 3a). The largest change in tilt of 16.3 µrad was observed at KSS
located 1.2 km north of the MC, followed by a change
at KSE of 12.7 µrad. Tilt changes at the foot of the
volcano (N.KSYV and N.KSNV) of 0.1–0.9 µrad were
somewhat complex, probably due to perturbation by
long-period events (Yamada et al. 2021). Vectors oftilt changes during the time period between 10:00 and
10:02 revealed a simple radial pattern from several
hundred meters north of the MC (Fig. 4a). After the

Table 1 Details of the instruments in the tiltmeter network
Station

Altitude
(m)

Depth
(m)

Sensor typeOperation

Installation(year)

KSS

1906

110

Akashi, JTS-33; Mitsutoyo, JTS-33

KSVO

1991; 2016

KSE

1893

50

Akashi, JTS-33

KSVO

2001

KSW

1883

220

Akashi, JTS-33

KSVO

2001

YNW

1985

1

Applied Geomechanics, 701-2A

KSVO

2016

JIE

1990

2

Applied Geomechanics, 701-2A

KSVO

2016

N.KSHV

1520

198

Mitsutoyo, JTS-33

NIED

2012

N.KSYV

1057

199

Mitsutoyo, JTS-33

NIED

2014

N.KSNV

845

200

Mitsutoyo, JTS-33

NIED

2014

V.KSAO

1776

95

Mitsutoyo, JTS-33

JMA

2011

Terada et al. Earth, Planets and Space(2021) 73:157(a)Page 5 of 17(b)Fig. 3 Tilt components recorded at 09:50–10:10 (JST) on 23 January 2018. a North–south and b east west components. Acausal low-pass filterswith cut-off frequencies of 0.034 and 0.021 Hz were applied to the tilt records of JTS-33 (KSE, KSS, KSW, N.KSHV, N.KSYV, and N.KSNV) and 701-2A(YNW and JIE), respectively. c Raw vertical ground velocity waveform at KSS. See Fig. 1 and Table 1 for station names and locations. We term theinflation and deflation phases as phase 1 and phase 2, respectively. The time of eruption onset was estimated by Yamada et al. (2021)

onset of the eruption at 10:02:09 (Yamada et al. 2021),
the vectors of the tilt change reversed. The speed of
the tilt changes gradually decreased, and thus the timing of cessation of the tilt changes is uncertain. We
assumed that tilt changes after 10:20 were negligible. In this study, we term the inflation and deflation
phases as phase 1 (10:00–10:02) and 2 (10:02–10:20),
respectively (Fig. 3).
The N–S component at KSS showed no significant
change prior to 10:01 (Fig. 4b), suggesting a somewhat
complex process during the initial stage of the inflation. However, details of this tilt change were obscuredby the following large tilt changes. We therefore analyzed the tilt vectors of phases 1 and 2 (Fig. 4a).

Modeling

To determine the volume changes during the inflation (phase 1) and deflation (phase 2), we modeled thedeformation source causing the tilt changes. Although
the tilt vectors exhibit a simple radial pattern, the tiltchange at N.KSHV was twice as large as that recorded
at N.KSYV, even though the two sites are located at asimilar distance from the center of the radial pattern(Fig. 4). The pattern of tilt change indicates the existence of a subvertical crack rather than a simple point
source. Furthermore, the new vent chain, including
MC, WC, and another smaller crater, formed over a
distance of 450 m and aligned along an azimuth of 112°
(Figs. 2d and 4). On the basis of these observations, weTerada et al. Earth, Planets and Space(2021) 73:157

Page 6 of 17

Fig. 4 Tilt vectors in the study area. Yellow rectangles show the localities of vents formed during the 2018 eruption. White, red, and gray circles
represent hypocenters for the periods from 1 January 2016 to 22 January 2018 (before the eruption), from 23 January to 22 April 2018 (after the
eruption), and from 22 to 23 April 2018 (representative of SPCG unrest in 2018), respectively, as determined by the Kusatsu–Shirane Volcano
Observatory, Volcano Fluid Research Center, Tokyo Institute of Technology. a Phase 1 (black arrows, Fig. 3) and phase 2 (white arrows, Fig. 3). b First
half of phase 1 (black arrows; 10:00–10:01 JST). Note that the vector scale is different from that in Fig. 4a

Terada et al. Earth, Planets and Space(2021) 73:157

Page 7 of 17

The strike of the crack (A) and dip angle (D) weresearched from 70° to 160° and 15° to 165°, respectively, atintervals of 1°.

propose that a subvertical tensile crack exists beneath
the Kagamiike-kita PC.
We modeled a finite rectangular tensile crack in anelastic and homogeneous half-space medium as proposedby Okada (1992). This model comprises the following
eight parameters: spatial location of the crack center (X,Y, and Z); horizontal length (L) and vertical width (W)of the crack; strike of the crack (A) relative to north; dipangle (D) relative to the ground surface; and dislocation(U) of the crack. Our model assumes Lamé’s constant  =µ and L/W = 2. We searched for the best-fit combinationof seven parameters (X, Y, Z, L, A, D, and U) that produced the smallest misfit (root mean square sum of thedifference between the observed and calculated values ofeach tilt vector) by the grid-search method.
We searched for the crack center (X, Y, and Z) within
500 m horizontally of the MC and at 500–2100 m in elevation (0–1600 m asl) at intervals of 50 m. The horizontal length of the crack (L) and dislocation were searchedfrom 200 to 2400 m and 0.01 to 2.60 m, respectively. The
search interval for L was 100 m, and for U was 0.01 m.

Results

Figure 5 compares the observed and modeled tilt vectors in phase 1. The N–S component at JIE was excluded
because of the absence of data at this site. The observed
tilt records are well fitted by a crack model with optimal
parameters. The optimal position of the center of thecrack is located beneath the center of the MC at 1000 m
depth (1100 m asl) with a strike (A) of 107°, which is consistent with the alignment of the MC (112°; Fig. 4). An
optimal dip angle (D) of 88°N suggests the crack dips
steeply to the north. Such a subvertical crack can cause
uplift on both its northern and southern sides. Upliftof the northern side is consistent with the tilt changes
recorded around YCL, while uplift on the southern sidecan explain the large tilt change at N.KSHV relative to
N.KSYV. The optimal horizontal length (L), vertical width(W), and dislocation (U) are 1700 m, 850 m (assuming

(b)

(a)

10

20

10
10

-10

0

10
20
20

0

2 .0
1.5

10

0

1

0.5

0

20

EW [km]EW [km]Fig. 5 Comparison of observed (red arrows) and modeled (blue arrows) tilt vectors of phase 1 (Fig. 3). The location of the modeled crack isrepresented by the yellow rectangle. Modeled vertical displacement is projected onto a virtual plane at 2100 m asl (corresponding to an altitudeof MC) at intervals of 5 mm. Solid and dashed contour lines correspond to vertical uplift and subsidence, respectively. White, red, and gray circles
represent hypocenters determined by the KSVO for the periods from 1 January 2016 to 22 January 2018 (before the eruption), from 23 January to
22 April 2018 (after the eruption), and from 22 to 23 April 2018 (representative of SPCG unrest in 2018), respectively. Optimal crack parameters are
summarized in the lower right of the figure. The square with the dashed line represents the area shown in Fig. 6

Terada et al. Earth, Planets and Space(2021) 73:157

this is half of the horizontal length), and 0.35 m, respectively, which resulted in a volumetric expansion of ~ 5.1
× ­105 ­m3.
To assess the uncertainty of our estimation of optimal
crack location, we computed the spatial distribution of
misfit values in horizontal and vertical cross-sectionsat intervals of 100 m (Fig. 6). For each panel, we fixed Z(horizontal cross-section), X (N–S cross-section), and Y(E–W cross-section) to optimum values of 1100 m asl,
36.629°(latitude), and 135.540°(longitude), respectively.
Using the same method as in Section “Modeling”, wesearched for the best-fit parameters L, A, D, U, and location (X and Y for horizontal cross-section, Y and Z for
N–S cross-section, and X and Z for (E–W cross-section)

Page 8 of 17

within the range described in Section “Modeling”. If we
accept a misfit of < 1 µrad, the central value of the uncertainty of the crack location is in the range of 600 m in the
N–S direction, 300 m in the E–W direction in the horizontal cross-section, and 500 m in the vertical cross-sections (Fig. 6). The corresponding uncertainty in volume
expansion is in the range of 2.8 × ­105 to 6.1 × ­105 ­m3.
To assess the degree of uncertainty arising from theassumption in L/W = 2, we searched for the best-fit
combination of five parameters (Z, L, A, D, and U) using
the same method as above. The horizontal location of
the crack center (X and Y) is fixed at the MC, whichcorresponds to the optimal value obtained assuming L/W = 2. Additional file 1: Figure S1 compares theFig. 6 Spatial distribution of misfit values in horizontal and vertical cross-sections at intervals of 100 m. For each panel, we fixed Z (horizontal
cross-section), X (N–S cross-section), and Y (E–W cross-section) to optimum values of 1100 m asl, 36.629° (latitude), and 135.540° (longitude),respectively. Using the same method as in Section “Modeling”, we searched for the best-fit parameters L, A, D, U, and location (X and Y for horizontal
cross-section, Y and Z for N–S cross-section, and X and Z for (E–W cross-section) within the range described in Section “Modeling”. Yellow lines are
topographic contours at intervals of 20 m. Yellow squares indicate the locations of craters formed during the 2018 eruption. White and red circles
represent hypocenters determined by the KSVO for the periods from 1 January 2016 to 22 January 2018 (before the eruption) and from 23 January
to 22 April 2018 (after the eruption), respectively

Terada et al. Earth, Planets and Space(2021) 73:157

Page 9 of 17

-0
.

-10

-5

-5

5

0

-10
-5

-10

-15

-1
5
-5

-10
-1

0

-0.5
-5
-5

0
-1

EW [km]0-1

EW [km]Fig. 7 Comparison of observed (red arrows) and modeled (blue arrows) tilt vectors of phase 2 (Fig. 3). The location of the modeled crack isrepresented by the yellow rectangle. Modeled vertical displacement is projected onto a virtual plane at 2100 m asl (corresponding to an altitudeof MC) at intervals of 1 mm. Solid and broken contour lines correspond to vertical uplift and subsidence, respectively. White, red, and gray circles
represent hypocenters determined by the KSVO for the period from 1 January 2016 to 22 January 2018 (before the eruption), from 23 January
to 22 April 2018 (after the eruption), and from 22 to 23 April 2018 (representative of SPCG unrest in 2018), respectively. Vertical N–S and E–Wcross-sections of hypocenters are also shown. The optimal dislocation of the crack is shown at lower right in the figure; other crack parameters are
as in Fig. 5

observed and modeled tilt vectors in phase 1, assuming
L/W = 4, 1, and 0.5. In the cases given here, the centerof the crack ranged between 400 and 1100 m asl. The
volumetric expansion is estimated to be 2.8 × ­105 to 6.4
× ­105 ­m3. Misfits were minimized with the assumptionof L/W = 2.
Tilt vectors of phase 2 are in the opposite direction to
those of phase 1, implying contraction of the crack that
expanded during phase 1. For simplicity, we assumed
the crack location, size, dip, and angle are fixed, andsearched for the optimal dislocation in phase 2. Figure 7 compares the observed and modeled tilt vectors
in phase 2. The tilt record at JIE was not used because
there are no data after 10:02 due to an electrical problem. Similar to the results for phase 1, the modeled tilt
vectors are consistent with the observations. The dislocation in phase 2 is estimated to be ~ – 0.25 m, implying a deflation volume ­VDef of ~ 3.6 × ­105 ­m3. Therefore,
the net volume inflation V­ Net during the 2018 eruption
was ~ 1.5 × ­105 ­m3.

Vapor release from the main crater
Volcanic ash cloud during the 2018 eruption

Although vigorous tephra and gas emission occurredat the beginning of the eruption (Fig. 2b), photographs(Fig. 2c), volcanic tremors, infrasound (Yamada et al.
2021), and tilt changes suggest that the eruption was
brief (i.e., ≤ 10–15 min). Surveillance video and photographs taken from the ground mostly focused on theeruption site, and did not record the entire volcanic
plume.
In such a situation, radar observations (Syarifuddin
et al. 2019; Marzano et al. 2020) can track an ash cloud.
Sato et al. (2018) presented a time-series of the ash cloudecho of the phreatic eruption at Aso volcano on 8 October 2016. The echo obtained from weather radar with awavelength of 5.6 cm (C-Band), operated by the Japan
Meteorological Agency (JMA), is consistent with the
distribution of the tephra deposit. The total mass of ash(3.2–7.5 × ­108 kg), as estimated from the cloud heightby the echoes (12,000 m) and duration of ash emission

Terada et al. Earth, Planets and Space(2021) 73:157

Page 10 of 17(160–220 s), is consistent with field observations (6.0–6.5
× ­108 kg; Sato et al. 2018).
During the 2018 MPCG eruption, C-Band weather
radar stations of the JMA clearly detected an echo of anash cloud at a range of 0–9 km from the MC at 10:05–
10:10 (Meteorological Research Institute 2018). The trackof the ash cloud echo coincides with the distribution of
the tephra deposit obtained by field geological surveys(Kametani et al. 2021). The most distal position of the ash
cloud was 9 km from MC at 10:05–10:10, which implies
a horizontal wind speed of ~ 20 m/s based on the onset
time of the eruption (10:02 JST). This estimate is consistent with speeds of 20–30 m/s recorded at heights of
3000–4000 m asl at two sites (Fig. 1a) by radiosondes at
09:00. Therefore, we estimated the change in height ofthe volcanic plume with distance from the MC (Fig. 7),
based on reflectivity images from the Nagano radar (bluecircle in Fig. 1a), located 66 km from Mt. Motoshirane.
The ash cloud was detected at a range of 2100–3000 m
asl at 2 km from the MC. The ash cloud ascended with
distance from the MC, reaching a height of 3100–5400 m
asl at 7 km from the MC.

Fig. 8 Tracks of the ash cloud. Red bars represent ash heightdetected by the Nagano weather radar (blue circle in Fig. 1a) at 2, 5,
and 7 km from the MC at 10:05–10:10 (JST) (Meteorological ResearchInstitute 2018). Geometries of a buoyant plume were calculated with
Eq. 1 and 2, assuming a wind speed of 20 m/s (solid lines) and 30 m/s(dashed lines). Values in the figure indicate the heat discharge rate q̇
in GW (­ 109 W), as defined in Eq. (2). The source of the plume is fixed
to the location of the MC

Qp = q̇t,where Qp is derived from the heat of vapor Qv and volcanic ash Qt:

Qp = Qv + Qt .

Simple plume model

To evaluate the mass of water vapor emitted from thevents during the 2018 eruption Mv (kg), we used a simple
model. Assuming a steady-state buoyant plume under aconstant horizontal wind, the height of a given positionin the plume increases with distance from the source tothe power of 2/3 (Briggs 1969; Kagiyama 1981):
2
3

(1)

h = Cx ,where h is the center of the height of the plume (m) and x
is the distance from the source (m). C can be estimated byleast-squares fitting of the observed the center of the ashheight against the distance from the vent. Alternatively, C
can be represented by an empirical equation using horizontal wind velocity and heat discharge rate, as follows:1

C = 3.8 × 10−5 u−1 q̇ 3 ,(2)

where u is horizontal wind speed (m/s) and q̇ is heat discharge rate (W).
We calculate plume height (h) as a function of distancefrom the MC (x) assuming a heat discharge rate (q̇ ) in therange of 1–500 GW under wind velocities (u) of 20 and
30 km/s, respectively (Fig. 8). Equations (1) and (2) areconsistent with the empirical observation that a plumewith strong buoyancy ascends rapidly, even in the case ofstrong horizontal wind (Briggs 1969).
If a steady-state plume is assumed for a duration of t ,the total heat energy emission Qp (J) can be calculated by

(3)

(4)

Assuming that most of the vapor condenses near the ventand releases latent heat to the plume, Qv is calculated as
the product of the total vapor mass Mv (in kg) and thespecific enthalpy ­Hv (kJ/kg) (Kagiyama et al. 1981; Terada
and Sudo 2012; Narita et al. 2019):

Qv = Hv Mv .

(5)

We used an ­Hv value of 2,675 J/kg for vapor at 100 °C
with a pressure of 1 × ­105 Pa (Wagner and Pruss 2002).
The temperature of the volcanic plume is unknown, but
a high temperature (e.g., more than several hundred
degrees Celsius) is unlikely, given the main constituent
minerals of volcanic ash emitted during the Motoshirane
eruption (Yaguchi et al., 2019). If the vapor temperature
is 300 °C, an ­Hv value of 3074 kJ/kg (Wagner and Pruss
2002) is used. We assume that the temperature dependence of ­Hv is sufficiently small to ensure the validity of the
discussion of water balance in Section “Discussion”.
Assuming that most of the heat of the tephra is converted to buoyancy near the vent, Qt can be calculated as

Qt = Ct TMt ,(6)

where ­Ct (J/kg/K) is the heat capacity averaged among
the tephras, T (K) is the temperature of the plume, andMt (kg) is the total weight of ejecta. If ­Hv, ­Ct, and T are
known or given, then q̇ can be estimated from Eq. (2),Terada et al. Earth, Planets and Space(2021) 73:157

and if t is assumed, then Mv can be calculated using theabove equations.
Evaluation of water vapor emission from the vents

Equations (1) and (2) were used to obtain plume shapes
for heat discharge rates. C can be estimated by leastsquares fitting of the plume track obtained by weather
radar (Section “Volcanic ash cloud during the 2018 eruption” ; red bars in Fig. 8). Assuming wind speed u of 20
and 30 m/s on the basis of measurements by radiosondes(red circles in Fig. 1a) at 09:00, the optimal heat discharge
rate is estimated to be 34 GW (u = 20 m/s) or 110 GW(u = 30 m/s). Considering the uncertainty in wind speed,
we evaluate the heat discharge rate q̇ to be 10–100 GW.
This value is > 1000 times higher than the representative
fumarolic heat discharge during non-eruptive periods at
active volcanoes in Japan (Kagiyama et al. 1981). Such alarge heat discharge of was estimated at Ontake volcano
24 h after its 2014 eruption (Yamaoka et al. 2016).
Assuming a duration for the emissions of 600 s, thetotal heat emission from the plume Qp is estimated to be
0.6–6.0 × ­1013 J (from Eq. 3). The total heat of the volcanic ash Qt is estimated to be 3 × ­1012 J, which was calculated by multiplying Mt = 3 × ­107 kg (Kametani et al.
2021) and a heat capacity of 1000 J/kg/K (i.e., representative of clay or volcanic rocks), with a temperature difference of 100 °C relative to ambient air. This leads to Qt = 3
× ­1012 J, which corresponds to 50% or 5% of Qp. Thus, themass of vapor emission is estimated to be 0.1–2 × ­107 kg,
based on Eqs. (3)–(6).

Discussion
During the 2018 MPCG eruption, we detected a progressive ground deformation by a borehole tiltmeter network.
Although satellite synthetic aperture radar (SAR) is useful for observing spatial surface deformation associatedwith a phreatic eruption (e.g., Kobayashi et al. 2018; Narita et al. 2020), a tiltmeter network with high temporal
resolution and dense spatial distribution provides detailsof progressive ground deformation during an eruption (e.g., Ueda et al. 2013; Aloisi et al. 2020; Zobin et al.
2020). Our crack model implies that inflation of 5.1 × ­105
­m3 beneath the main crater (MC) at the Kagamiike-kita
PC for 2 min was followed by deflation of 3.6 × ­105 ­m3.
We did not observe intense fumarolic activity or hightemperature volcanic gas emissions around the MC after
the eruption. The ejecta did not contain juvenile material(Yaguchi et al. 2019). Therefore, the rapid inflation/deflation was not caused by magma intrusion into the crack,
but was instead related to hydrothermal fluid. In this section, we discuss how the eruption occurred, based on aconceptual model of the hydrothermal system beneath
Kusatsu–Shirane volcano.

Page 11 of 17

Conceptual model of the hydrothermal system at Kusatsu–
Shirane volcano

A self-sealed zone beneath volcanoes, formed by theprecipitation of silica from aqueous fluids (Saishu et al.
2014), can play a key role in constraining the pore pressure of hydrothermal fluids (Fournier 1999). A self-sealed
zone can be located at the brittle–ductile (BD) transition
zone, because relatively low-permeability silica precipitates mainly at a temperature of 370–400 °C, as inferredfrom data from deep wells, quartz solubility calculations,and hydrothermal laboratory experiments (Fournier
1991, 1999; Saishu et al. 2014). This narrow self-sealed
zone separates the lithostatic pressure region from thearea where meteoric-derived hydrothermal fluids circulate through brittle crust at hydrostatic pressures (Fig. 9).
Such hydrothermal systems have been recognized basedon data including core and water samples from deep
wells in geothermal fields, such as at Kakkonda, Japan(Doi et al. 1998).
The bottom of a hypocenter distribution can delineate the BD transition zone, because the hypocenter distribution of volcanic–tectonic earthquakes can markthe boundary of quasistatic behavior (Fournier 1999;
Ingebritsen and Manning 2010; Castaldo et al. 2019). At
Kusatsu–Shirane volcano, most hypocenters are locatedshallower than 200 m asl (Figs. 5–7), and events below
sea level are extremely rare (Mori et al. 2006; Tseng et al.
2020; Yamada et al. 2021). Therefore, we propose that a
BD transition zone exists at 200 m asl beneath Kusatsu–
Shirane volcano (Fig. 9). Recent MT surveys havedetected a large and strong conductor between 500 m
asl and 1000 m below sea level (C2 of Matsunaga et al.
2020), which likely corresponds to a hydrothermal fluid
reservoir.
If a self-sealed BD transition zone episodicallybreaches, then sudden injection of supercritical hydrothermal fluid at a lithostatic pressure to the region ofhydrostatic pressure may cause failure of the overlying
brittle rock (Fournier 1999). The record of tilt changes(Figs. 5 and 7) suggests inflation of a subvertical crack
extending from above the BD transition zone to the surface, which provided a fluid pathway from depth.
The depth of a hydrothermal fluid reservoir associatedwith a phreatic eruption can be a key control on the magnitude of a phreatic eruption (Stix and de Moor 2018).
On the basis of the crack model and MT resistivity structure, we consider the fluid reservoir is located at 2 km
depth beneath the MPCG (Figs. 5 and 7). Compared withthe 2018 MPCG eruption, the 2014 Ontake eruption produced an order of magnitude more tephra (Maeno et al.
2016). Prior to the 2014 Ontake eruption, a tiltmeter andbroadband seismometer recorded rapid inflation of 1
× ­106 ­m3 at a depth of 1000 m, which was likely caused

Terada et al. Earth, Planets and Space(2021) 73:157

Page 12 of 17

Fig. 9 Schematic diagram of the hydrothermal system associated with the 2018 MPCG eruption. HF clusters represent hypocenter distributions ofhigh frequency volcanic–tectonic earthquakes. The depth of the self-sealed brittle–ductile transition zone is based on the hypocenter distribution(Figs. 5–7). 1 Matsunaga et al. (2020); 2 Tseng et al. (2020); 3 Yamada et al. (2021); 4 Ohba et al. (2019a, b); 5 Himematsu et al. (2020)

by water boiling (Maeda et al. 2017). The hydrothermal fluid associated with the 2014 Ontake eruption wasderived from 3–6 km beneath the ground surface (Katoet al. 2015; Narita et al. 2019), which is 2–4 km deeper
than for the 2018 MPCG eruption. Hydrothermal fluid
from a deeper reservoir is likely to have a higher specificenthalpy due to its closer proximity to magma, leading
more explosive eruption.
Deflation during phase 2

Tilt changes recorded during phase 2 indicate deflationof VDef = 3.6 × ­105 ­m3 over a period of 10–20 min (Fig. 7).
Assuming a density of 55.5 kg/m3, which is representativeof the vapor phase at the center of the crack at a hydrostatic pressure of 10 MPa (Wagner and Pruss 2002), thecorresponding mass is estimated to be 2 × ­107 kg, whichis comparable to the maximum mass estimate based on
the plume track (0.1–2 × ­107 kg). However, the fluid volume in the crack can be larger than VDef because of thecompressibility difference between hydrothermal fluid
and surrounding rock (e.g., Rivalta and Segall 2008).
Poroelasticity of the surrounding rock also affect the estimation of the actual fluid volume (e.g., Juncu et al. 2019).
Moreover, the fluid can be much denser than assumed
above, if the fluid is supercritical or a liquid phase. Given
these uncertainties, it is difficult to constrain the possible
mass based on the deflation volume. Therefore, the massof 2 × ­107 kg calculated above is a minimum estimate.
By comparing the mass Mv emitted as a plume andthe minimum deflation mass estimated from the deflation volume VDef, we infer some of mass associated withthe deflation was dispersed underground (Fig. 9). We

Terada et al. Earth, Planets and Space(2021) 73:157

Page 13 of 17

suggest that underground water acted as a buffer against
the impact of the sudden intrusion of hydrothermal
fluid, which absorbed the ascending hydrothermal fluidin the pre-existing crack. Gradual subsidence aroundthe Kagamiike PC after the 2018 MPCG eruption (Himenatsu et al. 2020) was likely caused by a relaxation ofgroundwater.
On 27 May 2011 (i.e., 2 months after the 2011 Mw 9.0
Tohoku Earthquake located 300–600 km from Kusatsu–
Shirane volcano; Ozawa et al. 2011), volcanic tremor(Yamada et al. 2021) and rapid inflation followed by
deflation (Fig. 10) were observed at the MPCG, similar to the ground deformation associated with the 2018
Motoshirane eruption. However, a phreatic eruption andunusual thermal activity did not occur in 2011. The tilt
vectors of the 2011 event (Fig. 10) are remarkably similar
to those of the 2018 MPCG eruption (Fig. 4). Although
only four tiltmeters were operating, we estimated thevolume changes during the inflation–deflation event in2011 using the same crack as for the 2018 MPCG eruption. As a result, we obtained an inflation of 5.8 × ­104 ­m3followed by deflation of 1.4 × ­104 ­m3. In this case we consider that the subvertical crack opened, but an eruption
failed to eventuate. Similar rapid inflation followed by
deflation without eruption was reported at Ontake volcano (Nakamichi et al. 2009) and Tokachidake volcano(Aoyama et al. 2020). On the basis of the analysis of a

VLP event in 2007, Nakamichi et al. (2009) proposed thatan enhanced flux of hot gases released from the magma
heated a hydrothermal system at 2000 m depth, whichcaused the vaporization of 8.8 × ­104 ­m3 (1.75 × ­106 kg) of
water, followed by discharge of vapor from the crack.
Heat associated with the 2018 eruption

Kusatsu–Shirane volcano is one of the most active volcanoes in Japan in terms of its persistent release of heat.
To assess the impact of the 2018 eruption on the energy
budget of Kusatsu–Shirane volcano, we estimated theheat released by the 2018 eruption and other thermal
activity (Fig. 9). The surface heat release from YCL is estimated to be 20 MW (Ohba et al. 1994), corresponding to
6 × ­1014 J/year. The largest phreatic eruption in the last
130 year, which occurred during 1937–1939, produced
heat of 5 × ­1014 J, as calculated from a tephra deposit of
5 × ­109 kg (Minakami 1939) based on an assumption ofa temperature difference of 100 °C and a heat capacity of1000 J/kg/K.
During the 2018 eruption, the heat discharge associated with the plume is estimated to have been 0.6–6.0
× ­1013 J (Section “Vapor release from the main crater”).
In addition, the tiltmeters suggest a net volume increase
VNet after the eruption of 1.5 × ­105 ­m3, which is large
relative to the heat release from the plume. Although thecompressibility of the thermal fluid filling the crack and

(a)

(b)

1

1

-0.35

2

-0
.

70

0

0
1

0

1

EW [km]

0

-0.35

0.35

-0.7

0

-0.3

5

EW [km]Fig. 10 Comparison of observed (red arrows) and modeled (blue arrows) tilt vectors on 27 May 2011. a Inflation phase and b deflation phase. The
location of the modeled crack is represented by the yellow rectangle. Modeled vertical displacement is projected onto a virtual plane at 2100 m asl(corresponding to an altitude of MC) at intervals of 0.2 mm and 0.07 mm, respectively. Solid and dashed contour lines correspond to vertical upliftand subsidence, respectively. Crack parameters such as X, Y, Z, L, A, and D are the same as for Fig. 5

Terada et al. Earth, Planets and Space(2021) 73:157

surrounding rocks is unknown (Rivalta and Segall 2008),
the deflation volume V­ Def provides a minimum estimate.
Assuming the thermal fluid is at boiling temperature at apressure of 10 MPa, which corresponds to the center ofthe crack at 1 km depth under hydrostatic pressure, thedensity and specific enthalpy of liquid water are 688 kg/
m3 and 1,420 kJ/kg, respectively (Wagner and Pruss
2002). As such, VNet is converted to a minimum mass of
1.0 × ­108 kg and minimum energy of 1.5 × ­1014 J (Fig. 9).
Although the atmospheric mass emission is small relative to past phreatic eruptions at Kusatsu–Shirane volcano, the heat energy of the 2018 MPCG eruption iscomparable to or larger than the annual heat output of
YCL or the largest eruption in the past 130 year. Therefore, the 2018 MPCG eruption is not negligible in termsof the energy budget of Kusatsu–Shirane volcano.
Implications for Kusatsu–Shirane volcano hazards and risks

Aligned small craters formed by past eruptions that are
10–30 m in horizontal scale are found throughout the
MPCG, suggesting that small phreatic eruptions haveoccurred repeatedly (Ishizaki et al. 2020). These small
craters are aligned NW–SE or WNW–ESE, which is consistent with the strike of the crack modeled in this study.
The direction of maximum compressive stress in thisregion is horizontal and oriented NW–SE or WNW–
ESE (Yoshida et al. 2012), allowing fluid to intrude from
depth through a crack with this orientation. We proposethat the pre-existing crack has repeatedly acted as a fluid
pathway from the self-sealed BD transition zone to thesurface for many years. Indeed, similar rapid inflation (5.8× ­104 ­m3) followed by deflation (1.4 × ­104 ­m3) without aphreatic eruption was observed on 27 May 2011 (Fig. 10).
The uncompensated deflation volume can be dispersedor absorbed by surrounding low-enthalpy undergroundwater, which acts as a buffer and suppresses explosions.
We infer that the inflation of hydrothermal fluid in thecrack in 2011 was too small to breach the surface rocks
above the crack.
Yamada et al. (2021) reported that the source locationof the volcanic tremor that occurred prior to the eruption
was the eastern side of APC (Fig. 1), 1 km north of the
MC on the Kagamiike-kita PC. The volcanic tremor was
likely caused by small shear fractures induced by a sudden intrusion of fluid (Yamada et al. 2021). The tilt vectorof KSS before 10:01 (Fig. 4b) suggests a pressure source
close to the tremor source near APC. We hypothesizethat two cracks, beneath the eastern side of APC andbeneath the Kagamiike-kita PC, were opened at the sametime, but an eruption of APC failed to eventuate.
Aligned small craters associated with past small phreatic eruptions are common at Kusatsu–Shirane volcano

Page 14 of 17(Ishizaki et al. 2020), but have been rarely documented
and monitored, particularly at the MPCG. Therefore,
these events pose a future risk to nearby tourists visiting
the area.
Relationship between the MPCG and SPCG

From March 2014 to 2016, the SPCG experienced amicroearthquake swarm, ground deformation indicating
an increase in volume of 1.2 × ­105 ­m3 over 20 months(Additional file 2: Fig. S2), an increase in the Cl and S­ O4
concentration of YCL, and an increase in the ­CO2/H2S
ratio of nearby fumaroles (SPCG unrest in 2014) (Ohbaet al. 2019b). Moreover, 3 months after the 2018 MPCG
eruption, a shallow inflation and a marked increase inseismicity occurred at SPCG from April 2018 (MPCGSPCG unrest in 2018). There is no clear evidence of dike
intrusion at shallow depths, but GNSS measurements
have shown that a sill-like source located a few kilometers
northwest of the SPCG at a depth of 5 km was inflated
during 2014–2016 and 2018–2020 (Munekane 2021). The
periods of inflations coincide with the period of SPCG
unrest in 2014 and MPCG–SPCG unrest in 2018. Thesefindings suggest that the supply of magmatic fluid fromthe sill-like source to shallow depths enhanced the volcanic activity of Kusatsu–Shirane volcano (Munekane2021).
A large fluid reservoir (C2 conductive zone) underlies
the region beneath the MPCG and SPCG (Matsunaga
et al. 2020), suggesting that magmatic fluid supply to theC2 conductive zone can affect the activity of the MPCG
as well as the SPCG. Tilt records have revealed a subvertical crack extending from the C2 conductive zone to thesurface (Fig. 5); therefore, we propose that the hydrothermal fluid that was supplied from the sill-like source
triggered the 2018 MPCG eruption, which highlights theclassification of Kusatsu–Shirane volcano as a “persistently restless volcano” (Roman et al. 2019) or “volcano
with slow unrest” (Stix 2018).
One possible mechanism for the triggering of the 2018
MPCG eruption is self-sealing of the fluid pathway to
the SPCG during the unrest. Temporary sealing prior toa phreatic eruption has been reported for many volcanoes (Christenson et al. 2010; Geirsson et al. 2014; Ohba
et al. 2019a; de Moor et al. 2019). It is plausible that thefluid pathway to the SPCG was sealed during the period
of quiescence after SPCG unrest in 2014, leading to theeruption at the MPCG in 2018 in response to the increasein magmatic fluid supplied from depth.
Three months after the 2018 MPCG eruption, the
SPCG underwent shallow inflation and a marked
increase in seismicity. Fumarolic activity and seismicity
in the MPCG showed a decline after the eruption, while

Terada et al. Earth, Planets and Space(2021) 73:157

Page 15 of 17

the sill-like source inflation continued, indicating that thefluid path to the MPCG had been tightly sealed after the
eruption. However, there are no observations to suggest
the mechanism or physicochemical process of this sealing. Sealing is one possible scenario, but further observations, such as defining the subsurface structure andmonitoring the C2 conductor, will be necessary to more
fully understand the conditions that determine whether
an eruption occurs, as well as the possible location of aphreatic eruption.

subvertical crack between the reservoir and groundsurface, causing the 2018 MPCG eruption.

Conclusions
A tiltmeter network recorded rapid inflation followed
by deflation during a phreatic eruption in 2018 at the
MPCG. The tilt records can be explained by a subvertical crack model centered at 1100 m asl beneath the
new crater. The strike of the crack is consistent withthe direction of chains of small craters formed by past
eruptions and the orientation of maximum compressive
stress in the region, allowing fluid to intrude from depth
through the crack. The inflation and deflation volumes
associated with the 2018 MPCG eruption are estimatedto be 5.1 × ­105 and 3.6 × ­105 ­m3, respectively, meaning that 1.5 × ­105 ­m3 of expanded volume remained
underground. We estimated the discharged water massin the volcanic plume to be 0.1–2.0 × ­107 kg, which issmaller than the mass associated with the deflation.
We propose that underground water acted as a buffer
against the sudden injection of hydrothermal fluid,which absorbed some of the hydrothermal fluid ascending through the crack. Gradual subsidence around theKagamiike PC after the 2018 MPCG eruption inferredfrom SAR was likely caused by a relaxation of intruded
groundwater at shallow depths. This is the second timethat rapid inflation followed by deflation has occurredat Kusatsu–Shirane volcano since the first tiltmeter was
established there in 1991. We suggest that similar small
phreatic or failed eruptions affected by groundwater
have repeatedly occurred in the MPCG. Although theejecta mass of the 2018 phreatic eruption is small, thetotal heat associated with the ground deformation isestimated to be ≥ ­1014 J, which is comparable with orexceeds the estimated annual heat release from YCL
and that from the largest eruption during the past
130 yr. Therefore, the 2018 MPCG eruption was not
negligible in terms of the energy budget of Kusatsu–
Shirane volcano. Prior to the 2018 eruption, precursorsof the eruption were not observed at MPCG. However,the SPCG had exhibited volcanic unrest since 2014. We
propose that enhanced amounts of magmatic fluid supplied from depth into the hydrothermal reservoir 2 km
beneath the MPCG and SPCG opened a pre-existing

Supplementary Information

Abbreviations
APC: Ainomine Pyroclastic Cone; Asl: Above sea level; BD: Brittle–ductile; HF:
High frequency; JMA: Japan Meteorological Agency; JST: Japan Standard Time;
KSVO: Kusatsu–Shirane Volcano Observatory, Tokyo Institute of Technology;
MC, WC, SC: Main, west, and south Crater of the 2018 eruption; MPCG: Motoshirane Pyroclastic Cone Group; MT: Magnetotelluric; NIED: National ResearchInstitute for Earth Science and Disaster Resilience; PC: Pyroclastic cone; SAR:
Synthetic aperture radar; SPCG: Shirane Pyroclastic Cone Group; YCL: Yugama
Crater Lake.

The online version contains supplementary material available at https://​doi.​
org/​10.​1186/​s40623-​021-​01475-4.
Additional file 1: Figure S1.Comparison of observed (red arrows) andmodeled (blue arrows) tilt vectors of phase 1 (Fig. 3), assuming L/W = a
0.5, b 1, and c 4. The location of the modeled crack in each case is represented by yellow rectangles. Modeled vertical displacement is projected
onto a virtual plane at 2100 m asl (corresponding to the altitude of theMC). Solid and dashed contours correspond to vertical uplift and subsidence, respectively. Optimum crack parameters are summarized at lower
right in each diagram. The optimal length (L), width (W), and dislocation(U) resulted in calculated volumetric expansions of ~ 2.8 × ­105 ­m3 (L/W =0.5), ~ 6.4 × ­105 ­m3 (L/W = 1), and ~ 4.2 × ­105 ­m3 (L/W = 4), respectively.
Additional file 2: Figure S2.Inflation source during the 2014 unrest of
the SPCG from March 2014 to October 2015. Optimal crack parameters
are summarized at lower right in the figure. We searched for the best-fit
combination of four parameters (Z, L, A, and U), assuming a sill-like source(D = 0°) and L/W = 2. The latitude and longitude of the center of the sill
were fixed at 36.645° and 138.538°, respectively, representing the centerof Mizugama crater. In this analysis, the smallest misfit (root mean squaresum of the difference between the observed and calculated values ofeach tilt vector, and each uplift at KSYG (44 mm) and KSE (25 mm) measured by GNSS) was calculated by the grid-search method as described in
Section “Modeling”. The optimal length (L), width (W), and dislocation (U)were 480 m, 240 m (assumed as half of the horizontal length), and 1.00 m,
respectively, which resulted in a calculated volumetric expansion of ~1.2
× ­105 ­m3.

Acknowledgements
We thank the National Research Institute for Earth Science and Disaster Resilience (NIED) for sharing the tilt records of N.KSHV, N.KSYV, and N.KSNV. The
Japan Meteorological Agency provided the tilt records from V.KSAO for the
event in 2011. The town of Kusatsu, Agatsuma District Forest Office, Ministry ofAgriculture, Forestry and Fisheries, and Ministry of Environment provided theapproval for our observations. We are grateful to the Kusatsu Tourism Cooperation for allowing us to use photographs. An early version of the tiltmeter
network at Kusatsu–Shirane volcano was established by Jun’ichi Hirabayashi inthe 1990s. We thank two anonymous reviewers for their valuable comments.
Some figures were drawn with Generic Mapping Tools (Wessel and Smith
1998). We used a 10-m mesh digital elevation model from the Geospatial
Information Authority of Japan (GSI) to draw the topographic maps.
Authors’ contributions
AT constructed the crack model, estimated water discharge, and drafted themanuscript. WK, YO, TY, MY, and SO helped to interpret the results of analyses
through discussions. TO, HA, and TT assisted the observation and helped to
improve the manuscript. 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, the Integrated Program for Next

Terada et al. Earth, Planets and Space(2021) 73:157

Generation Volcano Research and Human Resource Development and Grantin-Aid for scientific research (No. 17K20141).
Availability of data and materials
Tilt data from KSVO are available from the corresponding author on reasonable request. Tilt data from NIED and JMA can be downloaded from https://​www.​vnet.​bosai.​go.​jp.

Declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Author details
1
Volcanic Fluid Research Center, School of Science, Tokyo Institute of Technology, 2‑12‑1 Ookayama, Meguro, Tokyo 152‑8551, Japan. 2 Sakurajima Volcano
Research Center, Disaster Prevention Research Institute, Kyoto University,
1722‑19 Sakurajima‑Yokoyama, Kagoshima 891‑1419, Japan. 3 Research Center
for Prediction of Earthquakes and Volcanic Eruptions, Graduate School of Science, Tohoku University, 6‑6 Aramaki‑aza‑Aoba, Aoba‑ku, Sendai 980‑8578,
Japan. 4 Aso Volcanological Laboratory, Kyoto University, 5280 Kawayo,
Minami‑Aso, Kumamoto 869‑1404, Japan. 5 Institute of Seismology and Volcanology, Faculty of Science, Hokkaido University, N10W8 Kita‑ku, Sapporo,
Hokkaido 060‑0810, Japan. 6 Volcanology Research Department, Meteorological Research Institute, Japan Meteorological Agency, 1‑1 Nagamine, Tsukuba,
Ibaraki 305‑0052, Japan.
Received: 28 February 2021 Accepted: 1 July 2021

References
Aloisi M, Bonaccorso A, Cannavò F, Currenti G, Gambino S (2020) The 24
December 2018 eruptive intrusion at Etna volcano as revealed bymultidisciplinary continuous deformation networks (CGPS, boreholestrainmeters and tiltmeters). J Geophys Res Solid Earth. https://​doi.​org/​
10.​1029/​2019J​B0191​17
Aoyama H, Tanaka R, Hashimoto T, Murakami M, Narita S (2020) Seismic andgeodetic monitoring near the active crater at Tokachidake volcano(preliminary report). Geophys Bull Hokkaido Univ 83:25–48. https://​doi.​
org/​10.​14943/​gbhu.​83.​25
Barberi F, Bertagnini A, Landi P, Principe C (1992) A review on phreatic eruptions and their precursors. J Volcanol Geotherm Res 52:231–246. https://​doi.​org/​10.​1016/​0377-​0273(92)​90046-G
Briggs GA (1969) Plume Rise. Critical Review Series, Rep. TID-25075, Atomic
Energy Commission, Washington, DC, USA
Castaldo R, D’Auria L, Pepe S, Solaro G, Novellis VD, Tizzani P (2019) The impactof crustal rheology on natural seismicity: Campi Flegrei caldera casestudy. Geosci Front 10:453–466. https://​doi.​org/​10.​1016/j.​gsf.​2018.​02.​003
Christenson BW, Reyes AG, Young R, Moebis A, Sherburn S, Cole-Baker J, Britten K (2010) Cyclic processes and factors leading to phreatic eruption
events: Insights from the 25 September 2007 eruption through Ruapehu
Crater Lake, New Zealand. J Volcanol Geotherm Res 191:15–32. https://​doi.​org/​10.​1016/j.​jvolg​eores.​2010.​01.​008
de Moor JM, Stix J, Avard G, Muller C, Corrales E, Diaz JA, Alan A, Brenes J,
Pacheco J, Aiuppa A, Fischer TP (2019) Insights on hydrothermal-magmatic interactions and eruptive processes at Poás volcano (Costa Rica)from high-frequency gas monitoring and drone measurements. Geophys
Res Lett 46:1293–1302. https://​doi.​org/​10.​1029/​2018G​L0803​01
Doi N, Kato O, Ikeuchi K, Komatsu R, Miyazaki S, Akaku K, Uchida T (1998)
Genesis of the plutonic-hydrothermal system around quaternary granite
in the Kakkonda geothermal system, Japan. Geothermics 27:663–690.
https://​doi.​org/​10.​1016/​S0375-​6505(98)​00039-X

Page 16 of 17

Fournier RO (1991) The transition from hydrostatic to greater than hydrostatic
fluid pressure in presently active continental hydrothermal systems in
crystalline rock. Geophys Res Lett 18:955–958. https://​doi.​org/​10.​1029/​
91GL0​0966
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
Geirsson H, Rodgers M, LaFemina P, Witter M, Roman D, Muñoz A, Tenorio V,
Alvarez J, Jacobo VC, Nilsson D, Galle B, Feineman MD, Furman T, Morales
A (2014) Multidisciplinary observations of the 2011 explosive eruption of
Telica volcano, Nicaragua: Implications for the dynamics of low-explosivity ash eruptions. J Volcanol Geotherm Res 271:55–69. https://​doi.​org/​10.​1016/j.​jvolg​eores.​2013.​11.​009
Himematsu Y, Ozawa T, Aoki Y (2020) Coeruptive crustal deformation associated with the 2018 Kusatsu-Shirane phreatic eruption based on PALSAR-2
time-series analysis. Earth Planets Space 72:116. https://​doi.​org/​10.​1186/​
s40623-​020-​01247-6
Ingebritsen SE, Manning CE (2010) Permeability of the continental crust:
dynamic variations inferred from seismicity and metamorphism. Geofluids 10:193–205. https://​doi.​org/​10.​1111/j.​1468-​8123.​2010.​00278.x
Ishizaki Y, Nigorikawa A, Kametani N, Yoshimoto M, Terada A (2020) Geology
and eruption history of the Motoshirane Pyroclastic Cone Group, KusatsuShirane Volcano, central Japan. J Geol Soc Japan 126:473–491. https://​doi.​
org/​10.​5575/​geosoc.​2020.​0022
Juncu D, Árnadóttir Th, Geirsson H, Gunnarsson G (2019) The effect of fluid
compressibility and elastic rock properties on deformation of geothermal
reservoirs. Geophys Int J 217:122–134. https://​doi.​org/​10.​1093/​gji/​ggz011
Kagiyama T (1981) Evaluation methods of heat discharge and their applications to the major active volcanoes in Japan. J Volcanol Geotherm Res
9:87–97. https://​doi.​org/​10.​1016/​0377-​0273(81)​90016-0
Kametani N, Ishizaki Y, Yoshimoto M, Maeno F, Terada A, Furukawa R, Honda
R, Ishizuka Y, Komori J, Nagai M, Takarada S (2021) Total mass estimate of
the January 23, 2018 phreatic eruption of the Kusatsu–Shirane Volcano,central Japan. Earth Planets Space 73:141. https://​doi.​org/​10.​1186/​
s40623-​021-​01468-3
Kato A, Terakawa T, Yamanaka Y, Maeda Y, Horikawa S, Matsuhiro K, Okuda
T (2015) Preparatory and precursory processes leading up to the 2014
phreatic eruption of Mount Ontake Japan. Earth Planets Space 67:111.
https://​doi.​org/​10.​1186/​s40623-​015-​0288-x
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
Maeda Y, Kato A, Yamanaka Y (2017) Modeling the dynamics of a phreatic
eruption based on a tilt observation: Barrier breakage leading to the 2014
eruption of Mount Ontake, Japan. J Geophys Res Solid Earth 122:1007–
1024. https://​doi.​org/​10.​1002/​2016J​B0137​39
Maeno F, Nakada S, Oikawa T, Yoshimoto M, Komori J, Ishizuka Y, Takeshita Y,
Shimano T, Kaneko T, Nagai M (2016) Reconstruction of a phreatic eruption on 27 September 2014 at Ontake Central Japan, based on proximal
pyroclastic density current and fallout deposits. Earth Planets Space
68:82. https://​doi.​org/​10.​1186/​s40623-​016-​0449-6
Marzano FS, Mereu L, Scollo S, Donnadieu F, Bonadonna C (2020) Tephra mass
eruption rate from ground-based X-Band and L-Band microwave radars
during the November 23, 2013, Etna paroxysm. IEEE Trans Geosci Remote
Sens 58:3314–3327. https://​doi.​org/​10.​1109/​TGRS.​2019.​29531​67
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 fromthe resistivity structure of Kusatsu-Shirane Volcano. J Volcanol Geotherm
Res 390:106742. https://​doi.​org/​10.​1016/j.​jvolg​eores.​2019.​106742
Meteorological Research Institute (2018) Eruption cloud echoes from Mt.
Kusatsu-Shirane on January 23rd, 2018, observed by JMA weather radars.
Reports of Coordinating Committee for Prediction of Volcanic Eruption
129: 78–82.
Minakami T (1939) Explosive activities of volcano Kusatsu–Sirane during 1937
and 1938. Bull Earthq Res Inst 17:590–623
Mori T, Hirabayashi J, Nogami K, Onizawa S (2006) A new seismic observation
system at the Kusatsu-Shirane volcano. Bull Volcanol Soc Japan 53:41–47.
https://​doi.​org/​10.​18940/​kazan.​51.1_​41
Munekane H (2021) Modeling long-term volcanic deformation at the Kusatsu–
Shirane and Asama volcanoes, Japan using the GNSS coordinate time

Terada et al. Earth, Planets and Space(2021) 73:157

series. 21 April 2021, PREPRINT (Version 1) available at Research Square.
https://​doi.​org/​10.​21203/​rs.3.​rs-​433272/​v1
Nakamichi H, Kumagai H, Nakano M, Okubo M, Kimata F, Ito Y, Obara K (2009)
Source mechanism of a very-long-period event at Mt Ontake, central
Japan: response of a hydrothermal system to magma intrusion beneaththe summit. J Volcanol Geotherm Res 187:167–177. https://​doi.​org/​10.​1016/j.​jvolg​eores.​2009.​09.​006
Narita S, Murakami M, Tanaka R (2019) Quantitative relationship between
plume emission and multiple deflations after the 2014 phreatic eruption
at Ontake volcano Japan. Earth Planets Space 71:145. https://​doi.​org/​10.​
1186/​s40623-​019-​1124-5
Narita S, Ozawa T, Aoki Y, Shimada M, Furuya M, Takada Y, Murakami M (2020)
Precursory ground deformation of the 2018 phreatic eruption on IwoYama volcano, revealed by four-dimensional joint analysis of airborneand spaceborne InSAR. Earth Planets Space 72:145. https://​doi.​org/​10.​
1186/​s40623-​020-​01280-5
Nurhasan OY, Ujihara N, Tank SB, Honkura Y, Onizawa S, Mori T, Makino M(2006) Two electrical conductors beneath Kusatsu-Shirane volcano,Japan, imaging by audiomagnetotellurics and their implications forhydrothermal system. Earth Planet Space 58:1053–1059. https://​doi.​org/​
10.​1186/​BF033​52610
Ogawa Y, Aoyama H, Yamamoto M, Tsutsui T, Terada A, Ohkura T, Kanda W, Koyama T, Kaneko T, Ominato T, Ishizaki Y, Yoshimoto M, Ishimine Y, Nogami
K, Mori T, Kikawada Y, Kataoka K, Matsumoto T, Kamiisi I, Yamaguchi S, Ito
Y, Tsunematsu K (2018) Comprehensive survey of 2018 Kusatsu-Shirane
Eruption. Proc Symp Nat Disaster Sci 55:25–30
Ohba T, Hirabayashi J, Nogami K (1994) Water, heat and chloride budgets ofthe crater lake Yugama at Kusatsu-Shirane volcano, Japan. Geochem J
28:217–231. https://​doi.​org/​10.​2343/​geoch​emj.​28.​217
Ohba T, Hirabayashi J, Nogami K (2008) Temporal changes in the chemistryof lake water within Yugama Crater, Kusatsu-Shirane volcano, Japan:
implications for the evolution of the magmatic hydrothermal system.
J Volcanol Geotherm Res 178:131–144. https://​doi.​org/​10.​1016/j.​jvolg​eores.​2008.​06.​015
Ohba T, Yaguchi M, Nishino K, Numanami N, Daita Y, Sukigara C, Ito M, Tsunogai U (2019a) Time variations in the chemical and isotopic composition of fumarolic gases at Hakone volcano, Honshu Island, Japan, overthe earthquake swarm and eruption in 2015, interpreted by magma
sealing model. Earth Planets Space 71:48. https://​doi.​org/​10.​1186/​
s40623-​019-​1027-5
Ohba T, Yaguchi M, Nishino K, Numanami N, Tsunogai U, Ito M, Shingubara
R (2019b) Time variation in the chemical and isotopic composition of
fumarolic gasses at Kusatsu-Shirane Volcano Japan. Front Earth Sci 7:249.
https://​doi.​org/​10.​3389/​feart.​2019.​00249
Okada Y (1992) Internal deformation due to shear and tensile faults in a halfspace. Bull Seism Soc Am 82:1018–1040
Ozawa S, Nishimura T, Suito H, Kobayashi T, Tobita M, Imakiire T (2011)
Coseismic and postseismic slip of the 2011 magnitude-9 Tohoku-Oki
earthquake. Nature 475:373–376. https://​doi.​org/​10.​1038/​natur​e10227
Rivalta E, Segall P (2008) Magma compressibility and the missing source forsome dike intrusions. Geophys Res Lett 35:L04306. https://​doi.​org/​10.​1029/​2007G​L0325​21
Roman DC, LaFemina PC, Bussard R, Stephens K, Wauthier C, Higgins M, Feineman M, Arellano S, de Moor JM, Avard G, Cruz MM, Burton M, Varnam
M, Saballos A, Ibarra M, Strauch W, Tenorio V (2019) Mechanisms ofunrest and eruption at persistently restless volcanoes: Insights from the2015 eruption of Telica Volcano, Nicaragua. Geochem Geophys Geosys
20:4162–4183. https://​doi.​org/​10.​1029/​2019G​C0084​50
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
Sato H, Takahashi H, Yamamoto E, Fukuo N, Uehara M, Terasawa Y (1980)
Development of the crustal tilt observation method using borehole-typetiltmeters. Jisin 33:343–368. https://​doi.​org/​10.​4294/​zisin​1948.​33.3_​343
Sato E, Fukui K, Shimbori T (2018) Aso volcano eruption on October 8, 2016,
observed by weather radars. Earth Planets Space 70:105. https://​doi.​org/​
10.​1186/​s40623-​018-​0879-4

Page 17 of 17

Stix J (2018) Understanding fast and slow unrest at volcanoes and implications
for eruption forecasting. Front Earth Sci 6:56. https://​doi.​org/​10.​3389/​
feart.​2018.​00056
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/​s40623-​018-​0855-z
Syarifuddin M, Oishi S, Hapsari RI, Shiokawa J, Mawandha HG, Iguchi M (2019)
Estimating the volcanic ash fall rate from the Mount Sinabung eruption
on February 19, 2018 using weather radar. J Disaster Res 14:135–150.
https://​doi.​org/​10.​20965/​jdr.​2019.​p0135
Tanada T, Ueda H, Nagai M, Ukawa M (2017) NIED’s V-net, the fundamental
volcano observation network in Japan. J Disaster Res 12:926–931. https://​doi.​org/​10.​20965/​jdr.​2017.​p0926
Terada A (2018) Kusatsu-Shirane volcano as a site of phreatic eruptions. Jour
Geol Soc Japan 124:251–270. https://​doi.​org/​10.​5575/​geosoc.​2017.​0060
Terada A, Hashimoto T (2017) Variety and sustainability of volcanic lakes:
response to subaqueous thermal activity predicted by a numericalmodel. J Geophys Res Solid Earth 122:6108–6130. https://​doi.​org/​10.​
1002/​2017J​B0143​87
Terada A, Sudo Y (2012) Thermal activity within the western-slope geothermalzone of Aso volcano, Japan: Development of a new thermal area. Geothermics 42:56–64. https://​doi.​org/​10.​1016/j.​geoth​ermics.​2012.​01.​003
Tseng KH, Ogawa Y, Nurhasan TSB, Ujihara N, Honkura Y, Terada A, Usui Y,
Kanda W (2020) Anatomy of active volcanic edifice at the KusatsuShirane volcano, Japan, by magnetotellurics: hydrothermal implications
for volcanic unrests. Earth Planet Space 72:161. https://​doi.​org/​10.​1186/​
s40623-​020-​01283-2
Ueda H, Kozono T, Fujita E, Kohno Y, Nagai M, Miyagi Y, Tanada T (2013) Crustal
deformation associated with the 2011 Shinmoe-dake eruption as
observed by tiltmeters and GPS. Earth Planet Space 65:4. https://​doi.​org/​
10.​5047/​eps.​2013.​03.​001
Wagner W, Pruss A (2002) The IAPWS formulation 1995 for the thermodynamic
properties of ordinary water substance for general and scientific use. J
Phys Chem Ref Data 31:387–535. https://​doi.​org/​10.​1063/1.​14618​29
Wessel P, Smith WHF (1998) New, improved version of generic mapping tools
released. EOS Trans AGU 79:579. https://​doi.​org/​10.​1029/​98EO0​0426
Yaguchi M, Ohba T, Numanami N, Kawaguchi R (2019) Constituent mineral
and water-soluble components of volcanic ash from the 2018 eruptionof Mt. Motoshirane of Kusatsu-Shirane Volcano Japan. J Disaster Res
14:991–995. https://​doi.​org/​10.​20965/​jdr.​2019.​p0991
Yamada T, Kurokawa AK, Terada A, Kanda W, Ueda H, Aoyama H, OhkuraT, Ogawa Y, Tanada T (2021) Locating hydrothermal fluid injection of
the 2018 phreatic eruption at Kusatsu-Shirane volcano with volcanic
tremor amplitude. Earth Planet Space 73:14. https://​doi.​org/​10.​1186/​
s40623-​020-​01349-1
Yamamoto M, Aoyama H, Tsutsui T, Terada A, Kanada W, Ogawa Y (2018) Temporary seismic observation at Kusatsu-Shirane volcano, Japan. Abstractsof Japan Geoscience Union Meeting 2018 SVC41-P10.
Yamaoka K, Geshi N, Hashimoto T, Ingebritsen SE, Oikwa T (2016) Special issue“the phreatic eruption of Mt Ontake volcano in 2014.” Earth Planets Space
68:175. https://​doi.​org/​10.​1186/​s40623-​016-​0548-4
Yoshida K, Hasegawa A, Okada T, Iinuma T, Ito Y, Asano Y (2012) Stress before
and after the 2011 great Tohoku-oki earthquake and induced earthquakes in inland areas of eastern Japan. Geophys Res Lett 39:L03302.
https://​doi.​org/​10.​1029/​2011G​L0497​29
Zobin VM, Bretón M, Ramírez JJ, Santiago H (2020) Transition from passive
to pre-extrusion degassing of ascending andesitic magma before thelava dome-building eruption as derived from the seismic signals and tiltchanges: Volcán de Colima, México, August–September 2004. J Volcanol
Geotherm Res 401:106971. https://​doi.​org/​10.​1016/j.​jvolg​eores.​2020.​
106971

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