Introduction

In 2018, a phreatic eruption occurred on the Kagamiike-kita and Kagamiike pyroclastic cone (PC) of the Mt. Motoshirane Pyroclastic Cones Group (MPCG) on the southern part of Kusatsu–Shirane volcano, Japan (Ogawa 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 been found (Yaguchi et al. 2019).

Until several minutes before the occurrence of the eruption, no significant precursors (e.g., increased seismicity, seismic velocity change, unusual ground deformation and unusual thermal activity) (Barberi et al. 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 chemistry of Yugama Crater Lake (YCL) and nearby fumarolic gas 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 occurred at 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 relationship between the hydrothermal systems beneath the MPCG 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 and 2), were in operation at Kusatsu–Shirane volcano. Such a tiltmeter network with high temporal resolution and dense 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 a mass flux of water (Briggs 1969; Kagiyama 1981; Terada and Sudo 2012; Narita et al. 2019), although surveillance cameras and photographs taken from the ground did not constrain the height of the ash cloud at the climax of the eruption.

Fig. 1
figure 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 (both operated by the Japan Meteorological Agency). b Topographic map of Kusatsu–Shirane volcano. Solid lines are topographic contours at intervals of 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 the ground 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

Full size image
Fig. 2
figure 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 photograph

Full size image

In this paper, we present a simple model of inflation–deflation during the 2018 MPCG eruption based on the tilt 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 underground structure provides a constraint of a fluid pathway beneath Kusatsu–Shirane volcano. On the basis of the underground structure, the observations during the eruption in 2018 at the MPCG and the progression of the unrest of the SPCG since 2014, we propose a schematic model of the 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; Terada and 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, as a 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 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 continuous seismic activity has been detected around the APC and MPCG (Mori et al. 2006). On 23 January 2018, a phreatic 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 the northern 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 height of the volcanic cloud was estimated to be 5500 m asl by weather radar (Meteorological Research Institute 2018), corresponding to 3400 m in height relative to the altitude 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) showed a whitish plume and black tephra deposit around the Kagamiike-kita PC, but no ash cloud (Fig. 2c). During the eruption, the MC, the west crater (WC), and other smaller craters were formed on the Kagamiike-kita 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 Research Institute for Earth Science and Disaster Resilience (NIED) operate force-balanced-type tiltmeters (Sato et 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, and 5 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 rotation of the sensor. Tilt records are sampled at 20 Hz (NIED) or 200 Hz (KSVO). In this study, acausal low-pass filters with cut-off frequencies of 0.034 Hz and 0.021 Hz were applied 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 the 2018 eruption occurred but did record the 2011 event.

Table 1 Details of the instruments in the tiltmeter network
Full size table

Ground deformation detected by the tiltmeter network

We detected a simple and significant variation in the tilt 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 of tilt 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 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).

Fig. 3
figure 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 filters with 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 the inflation and deflation phases as phase 1 and phase 2, respectively. The time of eruption onset was estimated by Yamada et al. (2021)

Full size image
Fig. 4
figure 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

Full size image

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 obscured by 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 the deformation source causing the tilt changes. Although the tilt vectors exhibit a simple radial pattern, the tilt change at N.KSHV was twice as large as that recorded at N.KSYV, even though the two sites are located at a similar 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 propose that a subvertical tensile crack exists beneath the Kagamiike-kita PC.

We modeled a finite rectangular tensile crack in an elastic and homogeneous half-space medium as proposed by 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; dip angle (D) relative to the ground surface; and dislocation (U) of the crack. Our model assumes Lamé’s constant \(\uplambda\) = \(\upmu\) and L/W = 2. We searched for the best-fit combination of seven parameters (X, Y, Z, L, A, D, and U) that produced the smallest misfit (root mean square sum of the difference between the observed and calculated values of each 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 searched from 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. The strike of the crack (A) and dip angle (D) were searched from 70° to 160° and 15° to 165°, respectively, at intervals of 1°.

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 the crack 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. Uplift of the northern side is consistent with the tilt changes recorded around YCL, while uplift on the southern side can 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 this is half of the horizontal length), and 0.35 m, respectively, which resulted in a volumetric expansion of ~ 5.1 × 105 m3.

Fig. 5
figure 5

Comparison of observed (red arrows) and modeled (blue arrows) tilt vectors of phase 1 (Fig. 3). 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 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

Full size image

To assess the uncertainty of our estimation of optimal crack location, we computed the spatial distribution of misfit values in horizontal and vertical cross-sections at 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”, 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”. 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.

Fig. 6
figure 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

Full size image

To assess the degree of uncertainty arising from the assumption 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, which corresponds to the optimal value obtained assuming L/W = 2. Additional file 1: Figure S1 compares the observed and modeled tilt vectors in phase 1, assuming L/W = 4, 1, and 0.5. In the cases given here, the center of 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 assumption of 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, and searched 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 VNet during the 2018 eruption was ~ 1.5 × 105 m3.

Fig. 7
figure 7

Comparison of observed (red arrows) and modeled (blue arrows) tilt vectors of phase 2 (Fig. 3). 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 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–W cross-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

Full size image

Vapor release from the main crater

Volcanic ash cloud during the 2018 eruption

Although vigorous tephra and gas emission occurred at 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 the eruption 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 cloud echo of the phreatic eruption at Aso volcano on 8 October 2016. The echo obtained from weather radar with a wavelength 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 height by the echoes (12,000 m) and duration of ash emission (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 an ash cloud at a range of 0–9 km from the MC at 10:05–10:10 (Meteorological Research Institute 2018). The track of 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 of the volcanic plume with distance from the MC (Fig. 7), based on reflectivity images from the Nagano radar (blue circle 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.

Simple plume model

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

$$h\, = \,Cx^{\frac{2}{3}} ,$$
(1)

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 by least-squares fitting of the observed the center of the ash height 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:

$$C\, = \,3.8\, \times \,10^{{ - 5}} u^{{ - 1}} \dot{q} ^{{\frac{1}{3}}} ,$$
(2)

where u is horizontal wind speed (m/s) and \(\dot{q}\) is heat discharge rate (W).

We calculate plume height (h) as a function of distance from the MC (x) assuming a heat discharge rate (\(\dot{q}\)) in the range of 1–500 GW under wind velocities (u) of 20 and 30 km/s, respectively (Fig. 8). Equations (1) and (2) are consistent with the empirical observation that a plume with strong buoyancy ascends rapidly, even in the case of strong horizontal wind (Briggs 1969).

Fig. 8
figure 8

Tracks of the ash cloud. Red bars represent ash height detected 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 Research Institute 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 \(\dot{q}\) in GW (109 W), as defined in Eq. (2). The source of the plume is fixed to the location of the MC

Full size image

If a steady-state plume is assumed for a duration of \(\Delta t\), the total heat energy emission Qp (J) can be calculated by

$$Q_{p} \, = \,\dot{q}\Delta t,$$
(3)

where Qp is derived from the heat of vapor Qv and volcanic ash Qt:

$$Q_{p} \, = \,Q_{v} + \,Q_{t} .$$
(4)

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

$$Q_{v} \, = \,{\text{H}}_{{\text{v}}} M_{v} .$$
(5)

We used an Hv value of 2,675 J/kg for vapor at 100 °C with a pressure of 1 \(\times\) 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

$$Q_{t} \, = \,{\text{C}}_{{\text{t}}} TM_{t} ,$$
(6)

where Ct (J/kg/K) is the heat capacity averaged among the tephras, T (K) is the temperature of the plume, and Mt (kg) is the total weight of ejecta. If Hv, Ct, and T are known or given, then \(\dot{q}\) can be estimated from Eq. (2), and if \(\Delta t\) is assumed, then Mv can be calculated using the above 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 least-squares 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 \(\dot{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 a large 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, the total 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, the mass 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 associated with 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 details of 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 high-temperature 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 a conceptual model of the hydrothermal system beneath Kusatsu–Shirane volcano.

Conceptual model of the hydrothermal system at Kusatsu–Shirane volcano

A self-sealed zone beneath volcanoes, formed by the precipitation 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 inferred from 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 the area where meteoric-derived hydrothermal fluids circulate through brittle crust at hydrostatic pressures (Fig. 9). Such hydrothermal systems have been recognized based on data including core and water samples from deep wells in geothermal fields, such as at Kakkonda, Japan (Doi et al. 1998).

Fig. 9
figure 9

Schematic diagram of the hydrothermal system associated with the 2018 MPCG eruption. HF clusters represent hypocenter distributions of high 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)

Full size image

The bottom of a hypocenter distribution can delineate the BD transition zone, because the hypocenter distribution of volcanic–tectonic earthquakes can mark the boundary of quasistatic behavior (Fournier 1999; Ingebritsen and Manning 2010; Castaldo et al. 2019). At Kusatsu–Shirane volcano, most hypocenters are located shallower 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 have detected 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 episodically breaches, then sudden injection of supercritical hydrothermal fluid at a lithostatic pressure to the region of hydrostatic 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 associated with 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 with the 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 and broadband seismometer recorded rapid inflation of 1 × 106 m3 at a depth of 1000 m, which was likely caused by water boiling (Maeda et al. 2017). The hydrothermal fluid associated with the 2014 Ontake eruption was derived from 3–6 km beneath the ground surface (Kato et 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 specific enthalpy due to its closer proximity to magma, leading more explosive eruption.

Deflation during phase 2

Tilt changes recorded during phase 2 indicate deflation of VDef = 3.6 × 105 m3 over a period of 10–20 min (Fig. 7). Assuming a density of 55.5 kg/m3, which is representative of the vapor phase at the center of the crack at a hydrostatic pressure of 10 MPa (Wagner and Pruss 2002), the corresponding mass is estimated to be 2 × 107 kg, which is 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 the compressibility 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 mass of 2 × 107 kg calculated above is a minimum estimate.

By comparing the mass Mv emitted as a plume and the minimum deflation mass estimated from the deflation volume VDef, we infer some of mass associated with the deflation was dispersed underground (Fig. 9). We suggest that underground water acted as a buffer against the impact of the sudden intrusion of hydrothermal fluid, which absorbed the ascending hydrothermal fluid in the pre-existing crack. Gradual subsidence around the Kagamiike PC after the 2018 MPCG eruption (Himenatsu et al. 2020) was likely caused by a relaxation of groundwater.

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 and unusual 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 the volume changes during the inflation–deflation event in 2011 using the same crack as for the 2018 MPCG eruption. As a result, we obtained an inflation of 5.8 × 104 m3 followed 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 that an enhanced flux of hot gases released from the magma heated a hydrothermal system at 2000 m depth, which caused the vaporization of 8.8 × 104 m3 (1.75 × 106 kg) of water, followed by discharge of vapor from the crack.

Fig. 10
figure 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 uplift and subsidence, respectively. Crack parameters such as X, Y, Z, L, A, and D are the same as for Fig. 5

Full size image

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 the heat 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 of a temperature difference of 100 °C and a heat capacity of 1000 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 the compressibility of the thermal fluid filling the crack and surrounding rocks is unknown (Rivalta and Segall 2008), the deflation volume VDef provides a minimum estimate. Assuming the thermal fluid is at boiling temperature at a pressure of 10 MPa, which corresponds to the center of the crack at 1 km depth under hydrostatic pressure, the density 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 is comparable 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 terms of 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 have occurred 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 this region 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 propose that the pre-existing crack has repeatedly acted as a fluid pathway from the self-sealed BD transition zone to the surface for many years. Indeed, similar rapid inflation (5.8 × 104 m3) followed by deflation (1.4 × 104 m3) without a phreatic eruption was observed on 27 May 2011 (Fig. 10). The uncompensated deflation volume can be dispersed or absorbed by surrounding low-enthalpy underground water, which acts as a buffer and suppresses explosions. We infer that the inflation of hydrothermal fluid in the crack in 2011 was too small to breach the surface rocks above the crack.

Yamada et al. (2021) reported that the source location of 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 vector of KSS before 10:01 (Fig. 4b) suggests a pressure source close to the tremor source near APC. We hypothesize that two cracks, beneath the eastern side of APC and beneath the Kagamiike-kita PC, were opened at the same time, but an eruption of APC failed to eventuate.

Aligned small craters associated with past small phreatic eruptions are common at Kusatsu–Shirane volcano (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 a microearthquake 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 SO4 concentration of YCL, and an increase in the CO2/H2S ratio of nearby fumaroles (SPCG unrest in 2014) (Ohba et al. 2019b). Moreover, 3 months after the 2018 MPCG eruption, a shallow inflation and a marked increase in seismicity occurred at SPCG from April 2018 (MPCG-SPCG 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. These findings suggest that the supply of magmatic fluid from the sill-like source to shallow depths enhanced the volcanic activity of Kusatsu–Shirane volcano (Munekane 2021).

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 the C2 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 the surface (Fig. 5); therefore, we propose that the hydrothermal fluid that was supplied from the sill-like source triggered the 2018 MPCG eruption, which highlights the classification 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 to a 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 the fluid pathway to the SPCG was sealed during the period of quiescence after SPCG unrest in 2014, leading to the eruption at the MPCG in 2018 in response to the increase in 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 the sill-like source inflation continued, indicating that the fluid 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 and monitoring 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 a phreatic 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 with the 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 estimated to 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 mass in the volcanic plume to be 0.1–2.0 × 107 kg, which is smaller 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 the Kagamiike PC after the 2018 MPCG eruption inferred from SAR was likely caused by a relaxation of intruded groundwater at shallow depths. This is the second time that rapid inflation followed by deflation has occurred at 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 the ejecta mass of the 2018 phreatic eruption is small, the total heat associated with the ground deformation is estimated to be ≥ 1014 J, which is comparable with or exceeds 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, precursors of 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 subvertical crack between the reservoir and ground surface, causing the 2018 MPCG eruption.