Introduction

Volcanic glow is a phenomenon in which the plume above the crater glows red during the nighttime due to scattering of red light radiating from the crater floor. Since the presence of volcanic glow needs the presence of a high-temperature source, such as high magma levels in the conduit (cf. Harris and Ripepe 2007), its brightness (intensity) and spectra contain information regarding the size and temperature of this source, which is often hidden from view by the crater rim. We here examine volcanic glow at Sakurajima volcano (Japan) and relate variation in glow to the source of Vulcanian explosions.

Temporal changes in brightness of volcanic glow have been correlated with rates of degassing and shallow system processes. For example, Harris and Ripepe (2007) reported temporal variations in vent temperature and gas flux at Stromboli (Italy), where strong activity was associated with glow, weak with no or weak glow. Patrick et al. (2010) quantified vent glow at Halemaʻumaʻu (Kīlauea, Hawaiʻi) using brightness value in color images, and suggested that the temporal fluctuation of glow reflected changes in vent geometry and exposure of the high-temperature lava surface. Patrick et al. (2010) also reported that glow variation was correlated with ground deformation and seismicity. Johnson et al. (2014) used the red intensity of the glow as a proxy for the gas flux from the vent at Santiaguito (Guatemala) and found that glow intensity correlated with ground inflation-deflation cycles, suggesting that deformation, and hence glow, cycles were controlled by gas accumulation and venting beneath the vent. Volcanic glow sometimes also occurs as a precursor to eruption, indicating high levels of high-temperature magma in the shallow conduit (Iguchi et al. 2017; Yamamoto et al. 2017). Takeo et al. (2013), for example, reported the appearance of volcanic glow before Vulcanian eruptions at Shinmoe-dake (Japan) and suggested that it corresponded to leakage of volcanic gas from the shallow part of the conduit.

Here, we report temporal changes in volcanic glow prior to Vulcanian eruptions at the Showa crater of Sakurajima volcano, Japan. The onset of Vulcanian explosions has been understood as resulting from the failure of a lava plug formed by degassing and crystallization of magma in the uppermost part of the conduit (Stix et al. 1997; Clarke et al. 2007; Wright et al. 2007; Giachetti et al. 2010; Miwa et al. 2013; Bain et al. 2019; Gaunt et al. 2020). Following the disruption of the plug, overpressurized gas-rich magma (hereinafter, referred to as a gas pocket, following Iguchi et al. [2008]) is released, and the underlying magma is fragmented by sudden decompression and evacuated as ash plumes (Self et al. 1979; Ishihara 1985; Druitt et al. 2002; Formenti et al. 2003; Ohminato et al. 2006; Iguchi et al. 2008; Burgisser et al. 2011; Yokoo et al. 2013; Clarke et al. 2015). The triggering processes for Vulcanian explosions are not yet fully understood and there are two different models: ‘bottom-up’ and ‘top-down’ (Iezzi et al. 2020). At Sakurajima, isotropic expansion and cylindrical contraction are considered to have occurred at a depth of approximately 2 km in the conduit (i.e., the depth of explosion earthquakes). Upward propagation of a pressure wave or movement of fluids then triggers the expansion of the gas pocket and subsequent disruption of the plug. This is the bottom-up model (Ishihara 1985; Tameguri et al. 2002). However, the link between the explosion earthquake at depth and expansion of the gas pocket has not been clarified. Thus, instead, the top-down model simply considers that the plug is disrupted when the overpressure that has accumulated at the top of the conduit exceeds the fracture strength of the plug, and a decompression wave then propagates downward in the conduit, inducing fragmentation of magma (Druitt et al. 2002; Formenti et al. 2003; Clarke et al. 2015). Although this model considers instantaneous disruption of the plug, the actual disruption process occurs over a finite time and is more complicated (Iguchi et al. 2008; Yokoo et al. 2009).

We analyzed 90 Vulcanian eruption videos captured by a high-sensitivity camera installed at the foot of Sakurajima volcano with a line of sight to Showa crater (Aizawa et al. 2016; Muramatsu et al. 2018). Volcanic glow was quantified by the red intensity (R-value), and the green-to-red intensity ratio (GR ratio) was used to examine the spectral change of the glow. Infrasound data obtained at the camera location were also compared with the temporal change of the glow and allow us to infer triggering processes for Vulcanian eruptions.

Data

Eruption videos

Eruption videos used for the analysis were acquired using a high-sensitivity charge-coupled device (CCD) camera (Hitachi KP-DE500, Hitachi Kokusai Electric) installed at Kurokami Station (KURM) of Kyoto University, 3.5 km east of Showa crater (Fig. 1a). KURM has a line of sight to Showa crater (Fig. 1b), which is located east of the summit crater of Minamidake (Minamidake crater; Fig. 1a). The location of KURM thus allows visual observation of the eruptions (Johnson and Miller 2014; Aizawa et al. 2016) and infrasound analysis with a relatively unobstructed propagation path (Kim et al. 2015; Ishii et al. 2020). GPS time-synchronized videos were captured at normal speed (30 fps) and were continuously recorded on a hard-disk drive (Aizawa et al. 2016; Muramatsu et al. 2018). The videos used here were acquired between December 2011 and September 2015, when Showa crater was very active with more than 1,000 Vulcanian eruptions per year. This enables us to compare a large number of eruption events where we selected eruption events based on the following criteria: (1) the event occurred during the nighttime, (2) observable glow existed before the eruption, and (3) the crater was not obscured by cloud. After visual inspection of 3,650 eruptions from December 2011 to September 2015, we chose 90 eruption videos as suitable for the analysis (Supplementary Table S1).

Fig. 1
figure 1

(a) Location of Sakurajima volcano and observation sites. The contour interval is 200 m. The map inset at the upper left shows the location of Sakurajima on Kyushu Island (SW Japan). There are two active craters: Showa crater (red circle) and Minamidake crater (black circle). (b) Snapshot from an eruption video (Event 10; Table S1). The white rectangle (20 × 80 pixels) indicates the region of interest (ROI). The white dashed line represents the crater rim. (c) Typical waveform for infrasound accompanying Vulcanian eruptions at Showa crater (Event 28; Table S1). The inset is an enlarged view of the infrasound preceding phase (gray-shaded part)

Full size image

Infrasound data

Infrasound data were obtained at KURM (Fig. 1a). The propagation distance from the vent to KURM is approximately 3.58 km (Supplementary Fig. S1). The data were recorded by an ACO Type 7144 microphone from December 2011 to May 2013, and by a Hakusan SI104 microphone from June 2013 to September 2015. Both microphones have flat responses above 0.1 Hz. Data were logged at a sampling frequency of 200 Hz by a 24-bit resolution logger. The observed infrasound waveforms typically have an impulsive compression phase, and a subsequent rarefaction phase and long-period coda (Fig. 1c). Waveforms are similar to those observed for Vulcanian eruptions from Fuego (Johnson et al. 2004; Marchetti et al. 2009), Lokon-Empung (Yamada et al. 2016), Mount Cleveland (Iezzi et al. 2020), Popocatépetl (Matoza et al. 2019), Suwanosejima (Yokoo and Iguchi 2010), Tungurahua (Ruiz et al. 2006; Anderson et al. 2018), Volcán de Colima (Arámbula-Mendoza et al. 2020), and Shinmoe-dake (Aniya et al. 2022). Vulcanian eruptions occasionally show a modest precursory pressure increase just before the main impulsive phase, which is called an infrasound preceding phase (Fig. 1c; Yokoo et al. 2009), and has been reported in Vulcanian eruptions from some volcanoes (Yokoo and Iguchi 2010; Yamada et al. 2016; Anderson et al. 2018; Iezzi et al. 2020). In our dataset, an infrasound preceding phase was found in 64 of 90 eruptions. All of the infrasound waveforms are shown in Supplementary Fig. S2.

Methods

Quantification of volcanic glow

Volcanic glow observed by our camera can be interpreted as visible red light from a red-hot source on the crater floor scattered by a steam plume above the crater. Following Yokoo et al. (2013), we use the intensity of the red color gun (R-value, 0 to 255, 8 bit) to quantify the intensity of the glow, and the green-to-red intensity ratio (GR ratio) to investigate relative temperature changes. The effect of external light (moonlight) was removed by subtracting the linear trends of the RGB color channels from the original trends (Witsil and Johnson 2020), and a median filter with a kernel of five was applied to reduce noise. The glow shows some spatial variation due to variations in the shape and location of the steam plume (Supplementary Videos S1 through S4). To reduce this effect, we calculated the average R- and G-values for a region of interest (ROI) covering the glow region in the image (white rectangle in Fig. 1b) for each frame (i.e., approximately every 0.03 s). The resulting time series was smoothed using a moving average with a 0.3-s time window to reduce temporal variation and to obtain the GR ratio. We applied a backward linear weighted moving average, which uses data from the past 0.3 s and puts more weight on recent data and less on past data. The weighted moving average has a shorter time delay than the non-weighted (simple) moving average. The calculated R-values and GR ratios were plotted along the time axis to examine temporal changes in the glow intensity and relative temperature, respectively.

The GR ratio method is based on color-ratio pyrometry using a color digital camera (Fu et al. 2010; Densmore et al. 2011; Kuhn et al. 2011; Sankaranarayanan et al. 2021). This method involves optical pyrometry to estimate the radiation temperature by taking the signal intensity ratio of multi-color channels (Grum and Becherer 1979; DeWitt and Nutter 1988; Levendis et al. 1992). The intensity of radiation emitted from a surface with emissivity \({\varepsilon }_{\lambda }\) is described by Planck’s law:

$$\begin{array}{c}I(\lambda ,T)={\varepsilon }_{\lambda }{I}_{\mathrm{b}}\left(\lambda ,T\right)=\frac{2{\varepsilon }_{\lambda }h{c}^{2}}{{\lambda }^{5}\left({e}^{hc/\lambda {k}_{\mathrm{B}}T}-1\right)},\end{array}$$
(1)

where \({I}_{\mathrm{b}}\) is the spectral radiance [W m−2 sr−1 nm−1] emitted by a blackbody per unit area, solid angle (steradian, sr), time, and wavelength. In addition, \(h\) is the Planck constant, \({k}_{\mathrm{B}}\) is the Boltzmann constant, \(c\) is the speed of light, \(T\) is the temperature of the radiation surface, and \(\lambda\) is wavelength. Signal intensity (\({S}_{\mathrm{R}}\), \({S}_{\mathrm{G}}\), \({S}_{\mathrm{B}}\)) measured by the CCD camera with an RGB waveband measurement system can now be represented as (Fu et al. 2010)

$$\begin{array}{c}{S}_{\mathrm{R}}=\uppsi {\int }_{\Delta \lambda }{F}_{\mathrm{R}}\left(\lambda \right)\cdot {\varepsilon }_{\lambda }\cdot {I}_{\mathrm{b}}\left(\lambda ,T\right)d\lambda , \\ {S}_{\mathrm{G}}=\uppsi {\int }_{\Delta \lambda }{F}_{\mathrm{G}}\left(\lambda \right)\cdot {\varepsilon }_{\lambda }\cdot {I}_{\mathrm{b}}\left(\lambda ,T\right)d\lambda ,\\ {S}_{\mathrm{B}}=\uppsi {\int }_{\Delta \lambda }{F}_{\mathrm{B}}\left(\lambda \right)\cdot {\varepsilon }_{\lambda }\cdot {I}_{\mathrm{b}}\left(\lambda ,T\right)d\lambda ,\end{array}$$
(2)

where \(\uppsi\) is the instrument constant, [\({F}_{\mathrm{R}}\left(\lambda \right)\), \({F}_{\mathrm{G}}\left(\lambda \right)\), \({F}_{\mathrm{B}}\left(\lambda \right)\)] are the spectral response functions for the camera (i.e., the spectral sensitivity of the color filter), and \(\Delta \lambda\) is the bandwidth of the color filter. The intensity ratio of \({S}_{\mathrm{G}}\) to \({S}_{\mathrm{R}}\) (GR ratio) is thus,

$$\begin{array}{c}\frac{S_{\mathrm G}}{S_{\mathrm R}}=\frac{\int_{\Delta\lambda}F_{\mathrm G}\left(\lambda\right)\;\cdot\;\varepsilon_\lambda\;\cdot\;I_{\mathrm b}\left(\lambda,T\right)d\lambda}{\int_{\Delta\lambda}F_{\mathrm R}\left(\lambda\right)\;\cdot\;\varepsilon_\lambda\;\cdot\;I_{\mathrm b}\left(\lambda,T\right)d\lambda}=k_{\mathrm c}\frac{I_{\mathrm b,\lambda_{\mathrm G}}(T)}{I_{\mathrm b,\lambda_{\mathrm R}}(T)}=f\left(T\right),\end{array}$$
(3)

where \({k}_{\mathrm{c}}\) is the calibration coefficient for the camera. Here, we assumed that emissivity (\({\varepsilon }_{\lambda }\)) is constant (gray-body assumption) over the visible region so that it was excluded from the integral. We examined the thermal response of the camera, and confirmed that the GR ratio is a monotonically increasing function of temperature in the range of 700 to 1,000 °C (Supplementary Figs. S3a and S3b). Although we obtain a relation between the GR ratio and temperature (Fig. S3b), estimation of the true temperature involves difficulties derived from the model assumptions and atmospheric effects (Harris 2013). We also examined the time response of the camera to incident light and confirmed that the camera responds instantaneously to changes in light intensity without significant time delay (Supplementary Figs. S3c and S3d). The time response of the camera is sufficiently fast that it should not be a problem in the analysis (cf. Harris 2013).

Infrasound processing

As a first step, to reduce long-period noise from wind and instrumental drift, we subtracted linear trends found in the 3 s before the onset of the explosion (Anderson et al. 2018). No filter was applied because our waveforms have a relatively high signal-to-noise ratio. To allow comparison with the glow time series, the infrasound waveforms were shifted to their origin time at the vent (reverse-time correction) by subtracting the propagation time (\({t}_{\mathrm{prop}}\)) from the arrival time at KURM. Propagation time was obtained from the propagation distance (\({d}_{\mathrm{prop}}\)) (approximately 3,580 m) and the effective sound speed (\({a}_{\mathrm{eff}}\)):

$$\begin{array}{c}{t}_{\mathrm{prop}}={d}_{\mathrm{prop}}/{a}_{\mathrm{eff}},\end{array}$$
(4)

and by assuming linear propagation of the acoustic wave. The propagation distance was estimated by assuming the ray path to be the shortest route between the vent and the observation site (Minakami et al. 1970; Ruiz et al. 2006; Yokoo et al. 2009) (Fig. S1). The crater topography slightly changed during the observation period, but this does not affect the estimation of the propagation distance. The effective sound speed is calculated as

$$\begin{array}{c}{a}_{\mathrm{eff}}=a+\varvec{v}\cdot\varvec{n}=\sqrt{\gamma {{R}_{\mathrm{g}}T}_{\mathrm{a}}}+\varvec{v}\cdot\varvec{n},\end{array}$$
(5)

where \(a\) is the atmospheric sound speed, \(\gamma\) and \({R}_{\mathrm{g}}\) are the heat capacity ratio (1.4) and the gas constant (287 J/kg/K) for air, respectively, \({T}_{\mathrm{a}}\) is the atmospheric temperature in Kelvin, \({\varvec{v}}\) is the wind velocity vector, and \({\varvec{n}}\) is the unit vector in the propagation direction from the vent to KURM (80° from the north). The atmospheric temperature, and the wind velocity and direction, were obtained from radiosonde data from the Kagoshima Meteorological Office (KMO) (~ 11.4 km away from the crater; Fig. 1a) at the closest time to the eruption. We used the average values between those at the height of the crater rim (760 m a.s.l.) and KURM (60 m a.s.l.) (Fig. S1). Additionally, we assessed the uncertainty of the estimated effective sound speed (\({a}_{\mathrm{eff}}\)) and the propagation time (\({t}_{\mathrm{prop}}\)) by comparing with estimations from two other weather stations: Kihoku (KHK, ~ 18.2 km west of the crater) and Makinohara (MKN, ~ 19.3 km to the west) (Fig. 1a). Although there is some variation, the estimated values are generally similar between KMO, KHK and MKN (Supplementary Fig. S4). The variations in estimated sound speed and propagation time are up to ~ 7.5 m/s and ~ 0.22 s, respectively. We thus consider that the timing uncertainty of the reverse-time corrected infrasound is up to 0.22 s. We use the propagation time estimated from the data at KMO for the reverse-time correction because it is the closest site to Sakurajima and incorporates a vertical gradient of atmospheric temperature and wind velocity.

Although we assumed linear propagation for the estimation of the propagation time, supersonic wave propagation and nonlinear effects could occur near the source (Ishihara 1985; Yokoo and Ishihara 2007; Medici et al. 2014). However, the shock wave associated with Vulcanian explosions at Showa crater decays rapidly and changes to a sonic wave in the proximity of the source (Medici et al. 2014). Maher et al. (2020) also suggested that nonlinear effects are not significant for acoustic signals from Showa crater at propagation distances greater than 2 km. Iezzi et al. (2020) reported that nonlinear propagation has no quantifiable impact on the infrasound arrival time based on a numerical simulation of Vulcanian explosions from Mount Cleveland with a source-to-receiver distance of 3.8 km. We thus assume that nonlinear propagation does not affect the estimation of the propagation time at KURM.

We classified the events into those that had an infrasound preceding phase (N = 64, 71% of all events) and others (N = 26). The waveforms, and average of each type, are shown in Supplementary Fig. S5. In addition, we estimated the duration (\({\tau }_{\mathrm{p}}\)) and the amplitude (\({\Delta P}_{\mathrm{p}}\)) of the infrasound preceding phase (Fig. 1c). These parameters are summarized in Supplementary Table S1.

Results

Figure 2 is an example of temporal changes in the R-value, GR ratio, and infrasound during a single explosion (Event 5; Table S1; Video S1). Glow intensity (R-value) rapidly increases approximately 1 s before the onset of the explosion (thick red arrow in Fig. 2d), which synchronizes well with the infrasound preceding phase (gray arrow in Fig. 2d), even considering the timing uncertainty of the reverse-time correction (gray bar in Fig. 2d). This short-term glow increase preceding the explosion was found in 28 of 90 events (31%). Examples are shown in Supplementary Fig. S6 and Videos S2 through S4. We refer to these events as type-A events (Fig. 3a). We also found 36 events that had no short-term glow increase, but did exhibit an infrasound preceding phase (type B; Fig. 3b). The other 26 events had neither an infrasound preceding phase nor a short-term glow increase (type C; Fig. 3c). All events with a short-term glow increase had an infrasound preceding phase. Our event classification and associated parameters are summarized in Supplementary Table S1.

Fig. 2
figure 2

Temporal change in volcanic glow and infrasound for Event 5 (Table S1; Video S1). The top panels (a–c) are snapshots from the eruption video at the times indicated by the vertical dashed lines in (d). (d) Temporal change in R-value (red line), GR ratio (navy-blue line), and infrasound (gray line). The R-value shows a long-term gradual increase (thin red arrow) before the onset of the eruption. Approximately 1 s before the explosion, the R-value and the GR ratio rapidly increase (thick red and navy-blue arrows), synchronized with the infrasound preceding phase (thick gray arrow). The gray bar represents the timing uncertainty of the reverse-time correction of infrasound. All times are Japanese Standard Time (JST: UTC + 9 h)

Full size image
Fig. 3
figure 3

Classification of events. The left-hand panels give examples for (a) type-A, (b) type-B, and (c) type-C events. The red and gray lines mark temporal changes in the glow and infrasound, respectively. The vertical dashed lines and label ‘Ex’ represent the onset of the explosion. The gray bars represent the timing uncertainty of the reverse-time correction of infrasound. The right-hand panels show (d) the relationship between \({\tau }_{\mathrm{p}}\) and \({\Delta P}_{\mathrm{p}}\), and box-and-whisker plots of the values of (e) \({\tau }_{\mathrm{p}}\) and (f) \({\Delta P}_{\mathrm{p}}\), respectively

Full size image

Figures 3d through 3f show the relationship between the duration (\({\tau }_{\mathrm{p}}\)) and the amplitude (\({\Delta P}_{\mathrm{p}}\)) of the infrasound preceding phase. We see no clear relationship between \({\tau }_{\mathrm{p}}\) and \({\Delta P}_{\mathrm{p}}\) (Fig. 3d). Duration tends to be longer when accompanying a short-term glow increase (Fig. 3e), but no such relation could be observed for \({\Delta P}_{\mathrm{p}}\) (Fig. 3f).

The GR ratio is higher when accompanied by a short-term glow increase (navy-blue arrow in Figs. 2d and S6). Figure 4 shows the average temporal changes in the R-value (red), the GR ratio (navy), and infrasound (gray) during type-A (Fig. 4a), type-B (Fig. 4b), and type-C events (Fig. 4c). Plots of the R-value and the GR ratio for each event are given in Supplementary Figs. S7 and S8, respectively. The GR ratio increases coincided with a short-term glow increase and an infrasound preceding phase in type-A events (navy-blue arrow in Fig. 4a). In contrast, there is no noticeable increase in type-B events (Fig. 4b).

Fig. 4
figure 4

Temporal change in R-value (red lines), GR ratio (navy-blue lines), and infrasound (gray lines) for (a) type-A, (b) type-B, and (c) type-C events. The solid lines are the averages of all events for each type. The underlying shaded areas behind the R-value and the GR ratio represent the variations around the average (sample standard deviation; \(\pm 1\sigma\)) as the uncertainty of the average line. The R-value is plotted as the relative change from 3 s before the onset of the explosion (vertical dashed lines with ‘Ex’ label). For the GR ratio, the baseline is subtracted and the graphs are then normalized by their maxima before taking the average

Full size image

We observe a long-term increase in R-value before eruption onsets (thin red arrow in Fig. 2d; Fig. 5). Although a gradual weakening beginning approximately 6 min before, and subsequent disappearance of glow, eruptions at Showa crater have been reported (Yokoo et al. 2013), such a gradual decrease was not found in our dataset. Instead, we find a gradual increase before eruption. We examined temporal changes for 6 min prior to all 90 eruption onsets and found that long-term glow increases for almost all events (N = 73, 81%). The duration of the long-term glow increase (Fig. 5a) ranges between 3 and 330 s with a median of 32 s (Fig. 5b). The GR ratio does not show a significant increase during the long-term glow increase (Fig. 5a).

Fig. 5
figure 5

(a) Example of long-term (LT) glow increase (Event 9). The red and navy-blue lines represent the R-value and the GR ratio, respectively. The duration of the long-term glow increase was measured as the preceding time from the onset of the explosion (indicated by a label ‘Ex’). (b) Histogram of duration of long-term glow increase

Full size image

In summary, we found two types of temporal change in volcanic glow at Showa crater: a short-term rapid change approximately 1 s before the onset of Vulcanian eruptions (short-term glow increase) and a gradual change several minutes to several seconds prior to the onset of an eruption (long-term glow increase).

Discussion

Before Vulcanian eruptions, the upper conduit may be sealed by a solidified lava plug (Ishihara 1985; Wright et al. 2007; Clarke et al. 2015). We assume that red-hot regions at the top of this plug cause the observed glow. This red-hot region could be due to exposure of high-temperature material due to fracturing or cracking of the plug. Such a situation has been observed at Sakurajima (Ishihara 1985), Shinmoe-dake (Ichihara et al. 2013), and Santiaguito volcano (Sahetapy-Engel and Harris 2009; Johnson et al. 2014; Hornby et al. 2019). This process can be modeled in terms of thermal radiation from the high-temperature lava surface (cf. Sahetapy-Engel and Harris 2009). The spectral radiant flux [W nm−1] emitted by a surface with emissivity \({\varepsilon }_{\lambda }\) and area \(A\) can be written as

$$\begin{array}{c}E\left(\lambda,T\right)=\int_A\int_\Omega{\varepsilon_\lambda\cdot I}_{\mathrm b}\cos\theta\;dAd\Omega=A\varepsilon_\lambda\pi I_{\mathrm b}\left(\lambda,T\right),\end{array}$$
(6)

where \(\Omega\) is the solid angle and \(\theta\) is the radiation angle measured from the vertical direction. Here, we assumed isotropic hemispherical radiation (Lambertian radiation; Brewster 1992; Harris 2013). In Eq. (6), \({\varepsilon }_{\lambda }\) and \(A\) are the emissivity of the lava surface and the area of the red-hot region, respectively. Following Eq. (6), we thus consider that the glow intensity is a function of the surface temperature and the area of the red-hot region. The increase in temperature and/or area, which may be caused by an increase in the number of fractures in the plug, causes intensification of the glow.

As a secondary process, the red light emitted by the red-hot region illuminates the steam plume which are subsequently scattered to the camera. The scattering intensity depends on the particle size, number density, and spatial distribution of the scatterer (Ishimaru 1978). Detailed discussions of the scattering process are presented in Appendix A. Although the particle size, number density, and homogeneity of the scatterer could not be constrained from our data, we assume that an increase in the volume of the steam plume also contributes to the increase in glow intensity. An increase in gas emission contributed to the intensification of the glow at Stromboli (Harris and Ripepe 2007) and Santiaguito (Johnson et al. 2014).

In summary, the temporal change of the glow is caused by the change in the surface temperature (\(T\)) and the area (\(A\)) of the red-hot region (fractures on the plug), and the volume of the steam plume (the scatterer). Based on this model, we can consider plausible mechanisms to explain the short- and long-term glow increases.

Short-term glow increase

The short-term glow increase coincides with an increase in the GR ratio (Figs. 2d, 4a, and S6), indicating that the short-term glow increase is caused by a temperature increase. This suggests exposure of high-temperature sources. We consider that opening of tensile cracks causes the temperature increase and thus the short-term glow increase. Short-term glow increase always coincides with an infrasound preceding phase (type A; Fig. 4a). Yokoo et al. (2009) suggested that the infrasound preceding phase is caused by swelling of the plug just before explosion, with swelling being induced by expansion of the gas pocket. Our observation suggests that swelling of the plug occasionally occurs with the opening of tensile cracks in the plug. On the other hand, 36 eruptions did not show a short-term glow increase despite exhibiting an infrasound preceding phase (these being the type-B events; Fig. 4b). In addition, 26 eruptions had neither an infrasound preceding phase nor a short-term glow increase (type C; Fig. 4c). Type-B and -C events imply that tensile cracks were not opened at the crater floor just before these explosions. This indicates that swelling of the plug does not necessarily result in the opening of tensile cracks. This may be due to rapid disruption of the plug, so that cracks do not have time to open. If swelling does not occur, the plug might be immediately disrupted without tensile cracks being opened.

The duration of the infrasound preceding phase (\({\tau }_{\mathrm{p}}\)) for type-A events tends to be longer than that for type-B events (Fig. 3e). Yokoo et al. (2009) suggested that disruption of the plug (explosion) occurs when the expansion rate for the gas pocket exceeds the fracture strength of the plug. Therefore, \({\tau }_{\mathrm{p}}\) may be related to the swelling time of the plug until the fracture strength is exceeded. If tensile cracks open during swelling of the plug (type-A event), then overpressure may leak from the gas pocket and the expansion rate would become lower. This may increase the time for the expansion rate to exceed the fracture strength of the plug, resulting in a longer \({\tau }_{\mathrm{p}}\). When tensile cracks do not open, the expansion rate may become higher and \({\tau }_{\mathrm{p}}\) may become shorter (type B). If the expansion is more rapid, the plug would be immediately disrupted without observable swelling (type C).

It has also been suggested that a time interval between the eruption and the previous eruption (i.e., repose time) relates to variation in Vulcanian eruption activity (e.g., Watt et al. 2007; Bain et al. 2019; Iezzi et al. 2020). We examined relationships between the repose time (Table S1) and the types of event, and found that the repose time for type-C events tends to be longer than for type-A and -B events (Fig. 6). In Strombolian systems, the existence of a high viscosity layer (HVL) in the uppermost part of the magma column is thought to modulate eruption style (Gurioli et al. 2014; Del Bello et al. 2015; Gaudin et al. 2017; Oppenheimer et al. 2020; Simons et al. 2020). Considering an analogy of the Strombolian HVL and Vulcanian plug, a longer repose time would promote the development of plug due to increasing viscosity and yield stress (Leduc et al. 2015; Capponi et al. 2017; Kurokawa et al. 2022), resulting in a thick and strong plug (Simons et al. 2020; Thivet et al. 2021). This regime might be related to type-C events in our observations. On the other hand, a shorter repose time can lead to a thin and weak plug. A thin HVL (plug) formed between explosions may involve tensile cracking and/or deformation with a gas release (Gurioli et al. 2014; Leduc et al. 2015; Suckale et al. 2016; Gaudin et al. 2017; Thivet et al. 2021). This regime would correspond to type-A and -B events in our observation. Alternatively, it has also been shown that sealing of a lava dome through the formation of tuffisites can affect the rate of Vulcanian explosions (Kendrick et al. 2016). If a longer repose time allows greater sealing of fractures, the strength of the plug may increase and possibly result in type-C events.

Fig. 6
figure 6

Relationship between the repose time and the type of events. The left-hand panel shows the relationship between the duration of the infrasound preceding time (\({\tau }_{\mathrm{p}}\)) and the repose time in minutes. The symbol color represents types of events. The vertical axis has a logarithmic scale. The right-hand panel shows box-and-whisker plots of the values of the repose time

Full size image

The amplitude of the infrasound preceding phase (\({\Delta P}_{\mathrm{p}}\)) does not exhibit a clear relation with the type of event, or to \({\tau }_{\mathrm{p}}\) (Figs. 3d and f). The amplitude may be related to vertical displacement of the plug due to swelling (Yokoo et al. 2009). This may imply that the plug displacement is affected by other parameters, such as overpressure of the gas pocket (Yokoo and Ishihara 2007), vent radius (Muramatsu et al. 2018), explosion depth scaled by energy (Bowman et al. 2014), and variations in gas flux (Harris and Ripepe 2007). We also examined relationships between \({\tau }_{\mathrm{p}}\), \({\Delta P}_{\mathrm{p}}\) and the maximum amplitude of infrasound, but could not find any clear relation (Supplementary Fig. S9).

Long-term glow increase

The long-term glow increase does not accompany a significant increase in the GR ratio (Fig. 5a), unlike the short-term glow increase, suggesting that it cannot be explained by a temperature increase. Like the short-term increase, the long-term increase in glow can be explained by an increase in the area of the red-hot region and/or the volume of the steam plume. Indeed, since the long-term glow increase lasts for up to several minutes (Fig. 5b), it may be affected by the fluctuation of the steam plume as well as the state of the plug.

To assess the contribution of the steam plume to the long-term glow increase, we extracted the contour of the glow region in each image (e.g., Figs. 7a through 7f). The glow in the image can be regarded as an anomaly in the R-value against the dark nighttime background across R is approximately zero. Thus, the glow region can be automatically detected by setting a threshold R-value, which was determined from the 95th percentile in the histogram of the R-values in the first frame of the video for each event (Supplementary Fig. S10). This image was converted to a binary array using the R-value threshold to extract the contour of the glow using the OpenCV functions of Python. Summing the R-value over the extracted glow region (i.e., R-value × glow area) was also used as a proxy for scattering intensity (\(\Sigma R\)).

Fig. 7
figure 7

Additional analysis of long-term (LT) glow increase. (a–f) Snapshots of spatial development of glow (Event 9; Table S1). The white lines in each snapshot are the contours of the glow. (g) Temporal changes in calculated parameters. The green, brown, red, and navy-blue lines indicate the apparent glow area, \(\Sigma R\), the R-value, and the GR ratio, respectively. These traces are normalized by the value at the onset of explosion (indicated by a solid line labeled ‘Ex’) and are shifted vertically for visibility. Vertical dashed lines correspond to the times of the snapshots (a–f)

Full size image

The results for Event 9 (Table S1) are displayed in Fig. 7g as an example. The apparent glow area and \(\Sigma R\) increase from approximately 30 s before the onset of the explosion, this being the long-term glow (R-value) increase, and suggesting an increase in the area of the glow and increase in scattering. The GR ratio, however, does not show a significant change during the long-term glow increase, indicating the absence of a significant temperature increase. Plots for all events (N = 73), and the average plot, are given in Supplementary Fig. S11. All show an increase in the R-value without an increase in the GR ratio. Based on these results, we suggest that the increase in fractures on the plug and subsequent gas emission (increase in the volume of the steam plume) both contribute to the long-term glow increase. This process may correspond to gas leakage from the gas pocket before the onset of Vulcanian eruptions, indicating decompression of the conduit (Iguchi et al. 2008; Yokoo et al. 2013; Takeo et al. 2013; Inza et al. 2014). The duration of gas leakage at Showa crater was estimated from ground deformation to be approximately 0 to 5 min (Yokoo et al. 2013), which is consistent with the duration of the long-term glow increase (Fig. 5b). Our observations thus support the occurrence and potential significance of gas leakage prior to Vulcanian eruptions.

Finally, we consider the role of long-term gas leakage in initiation of Vulcanian eruptions. The duration of the long-term glow increase ranges from 3 to 330 s (Fig. 5b), suggesting that decompression of the conduit due to gas leakage continues for a few seconds to hundreds of seconds before the onset of the explosion. This means that decompression proceeds gradually and there is a time lag between the onset of decompression and explosion. It has been shown that the fragmentation behavior of magma depends on the decompression speed and can occur with a significant time lag between the onset of decompression and fragmentation (Kameda et al. 2008, 2013). Slow decompression experiments on magma analogs have suggested that during decompression, volatile-rich and poorly crystalline magma become gradually supersaturated and eventually expand explosively when the pressure drop reaches the supersaturation threshold (Rivalta et al. 2013). Before Vulcanian eruptions at Showa crater, low-crystallinity and volatile-rich magma may exist beneath the plug (Miwa et al. 2013). We suggest that gradual decompression causes progressive supersaturation of magma beneath the plug, and explosive expansion occurs with a time delay to trigger expansion of the gas pocket and disruption of the plug. The longer duration of gas leakage (decompression) may result in a larger pressure drop, which would cause supersaturation to deeper parts of the conduit. This may partially explain the variability of the depth of explosion earthquakes at Sakurajima and other volcanoes characterized by Vulcanian activity (Iguchi et al. 2008).

Conclusion

We present a conceptual model to explain the temporal change in volcanic glow at Showa crater (Fig. 8). Analysis of the glow associated with 90 Vulcanian eruptions suggests that long-term (several seconds to minutes) gas leakage through fractures in the degassed magma plugging the upper conduit is a key process for the initiation of Vulcanian eruptions. This process may increase the volume of the steam plume above the crater and cause the long-term (3–330 s before the explosion) glow increase (Fig. 8). We suggest that gradual decompression of the conduit due to gas leakage induces progressive supersaturation of magma beneath the plug, resulting in explosive expansion of magma residing in the upper portion of the conduit, but with a time delay to the triggered explosion.

Fig. 8
figure 8

Conceptual model to explain observed temporal changes in volcanic glow. Increased fracturing of the lava plug and subsequent long-term gas leakage causes an increase in the amount of steam generated, leading to a change in long-term (LT) glow. Gas leakage induces slow decompression and supersaturation of volatile-rich magma beneath the plug. Just before the explosion, swelling of the plug can occur to cause opening of tensile cracks, which causes a short-term (ST) glow increase. The difference in the expansion rate of the gas pocket and the mechanical strength of the plug controls the disruption manner, where the gray-shaded part (labeled ‘Ex’) indicates the onset of explosion

Full size image

Of all events considered, only 28 (31%) were accompanied by rapid opening of tensile cracks due to swelling of the plug, and which coincided with an infrasound ‘preceding phase.’ This process may cause a significant temperature increase on the crater floor and glow increase within 1 s of the explosion (Fig. 8). However, the remaining 62 eruptions were not preceded by opening of tensile cracks, even when heralded by an infrasound preceding phase, indicating that opening of tensile cracks is sporadic and occurs only when the plug swells immediately before an explosion. These differences can be interpreted as being due to a difference in the expansion rate of the gas pocket below the plug, the gas flux, and/or the plug rheology and mechanical state (Fig. 8). We thus argue that volcanic glow observations can be an essential component in constraining the processes leading to Vulcanian explosions, informing on shallow system event-triggering processes.

Appendix A. Scattering process for volcanic glow

Here, we consider the scattering process during the occurrence of volcanic glow. The incident light from the red-hot region is scattered by the steam plume, which we assume consists of water particles. The particle radius likely ranges from approximately 5 to 20 μm, which are typical values for cloud and fog droplets (van de Hulst 1981). The wavelength of red light is approximately 625 to 750 nm, and may be smaller than the water particle size. In this case, scattering can be described using Mie theory (van de Hulst 1981). If the path length of the incident light (\({l}_{\mathrm{p}}\)) is longer than the mean free path for a photon (MFP), then the process is dominated by multiple scattering. The MFP is the inverse of the scattering coefficient (\({K}_{\mathrm{s}}\)): \(\mathrm{MFP}=1/{K}_{\mathrm{s}}=1/n\sigma\), where \(n\) is the particle number density, and \(\sigma\) is the scattering cross-section. The typical number density for clouds and fog is 108 m−3 (Ishimaru 1978), yielding an MFP on the order of 10 to 100 m. Considering the distance between the floor of the Showa crater and the typical height of the plume above the crater rim, the path length for light (\({l}_{\mathrm{p}}\)) will be greater than approximately 200 m (Fig. S1), and thus the optical depth defined by \({l}_{\mathrm{p}}/\mathrm{MFP}\) falls between 2 and 20. This corresponds to an optically thick condition (Brewster 1992), and multiple scattering will be dominant. In this case, scattering of the radiance, neglecting absorption, is described by the radiative transfer equation:

$$\begin{array}{c}\frac{dI\left({\varvec{r}},s\right)}{ds}=-{K}_{\mathrm{s}}I\left({\varvec{r}},s\right)+\frac{{K}_{\mathrm{s}}}{4\pi }{\int }_{4\pi }p\left(s,{s}^{^{\prime}}\right)I\left({\varvec{r}},{s}^{^{\prime}}\right)d{\omega }^{^{\prime}}+q\left({\varvec{r}},s\right),\end{array}$$
(7)

where \(I\left({\varvec{r}},s\right)\) is the radiant intensity along the direction \(s\), \(p\left(s,{s}^{^{\prime}}\right)\) is a phase function that defines the probability of a photon moving in direction \(s\) and being scattered in direction \({s}^{^{\prime}}\), and \(q\left({\varvec{r}},s\right)\) is the source term (Ishimaru 1978; Brewster 1992). The source term \(q\left({\varvec{r}},s\right)\) is proportional to the radiant flux from the light source (Eq. (6)). Unfortunately, we cannot quantify this process because the particle size, number density, and homogeneity of the plume cannot be constrained. In addition, if the steam plume contains larger particles, such as volcanic ash, the problem will become more complicated. However, we can conclude that scattering intensity depends on the scattering coefficient (\({K}_{\mathrm{s}}\)) and the spatial distribution of the steam plume, which will be controlled by the amount of steam and gas emitted through fractures in the plug.