Comparison of volcanic explosions in Japan using impulsive ionospheric disturbances

Using the ionospheric total electron content (TEC) data from ground-based Global Navigation Satellite System (GNSS) receivers in Japan, we compared ionospheric responses to five explosive volcanic eruptions 2004–2015 of the Asama, Shin-Moe, Sakurajima, and Kuchinoerabu-jima volcanoes. The TEC records show N-shaped disturbances with a period ~ 80 s propagating outward with the acoustic wave speed in the F region of the ionosphere. The amplitudes of these TEC disturbances are a few percent of the background absolute vertical TEC. We propose to use such relative amplitudes as a new index for the intensity of volcanic explosions.


Introduction
The Earth's ionosphere ranges from ~ 60 to > 800 km in altitude and is characterized by large number of free electrons. Ionospheric conditions are controlled by solar radiation and often disturbed by geomagnetic activities. In addition to such disturbances caused by space weather, the ionosphere is disturbed by events below (Blanc 1985), such as earthquakes (Heki 2021), tsunami (Occhipinti et al. 2013), human-induced explosions (Kundu et al. 2021), and volcanic eruptions.
Ionospheric total electron content (TEC) can be easily measured by comparing the phases of two microwave carriers from global navigation satellite system (GNSS) satellites, such as Global Positioning System (GPS) (e.g., Hofmann-Wellenhof et al. 2008). Ground GNSS networks have been deployed to monitor crustal movements, and these networks were found useful to study ionospheric disturbances by volcanic eruptions. There are two types of ionospheric TEC responses to volcanic eruptions.
The first type is the long-lasting harmonic TEC oscillations ( Fig. 1). They are atmospheric modes excited by continuous acoustic waves generated typically by Plinian eruptions. The interference of upward and downward acoustic waves between the ground surface and the mesopause causes resonant oscillation of atmosphere. They have prescribed frequencies reflecting the vertical atmospheric structure (Tahira 1995). It was found after the 13 July 2003 eruption of the Soufrière Hills volcano, Montserrat, in the Lesser Antilles (Dautermann et al. 2009a, b).
This type of disturbance also occurred after the February 2014 eruption of the Kelud volcano, eastern Java Island, Indonesia (Nakashima et al. 2016). They reported that harmonic oscillations caused by atmospheric resonance excited by the Plinian eruption of the Kelud volcano lasted for ~ 2.5 h after the eruption started. Shults et al. (2016) found similar TEC oscillations after the 2015 April Plinian eruption of the Calbuco volcano, Chile. Cahyadi et al. (2020) also found such harmonic TEC oscillations lasting ~ 20 min following the 2010 November 5 eruption of the Merapi volcano, central Java Island. Although the onsets of these continuous eruptions are not always clear, the TEC oscillations emerge 20-30 min after the eruptions started. Cahyadi et al. (2020) also suggested that the TEC oscillation amplitudes relative to background TEC represent the mass eruption rate, and the products of such amplitudes and the duration provides a new index for the total amount of the ejecta. Recently, Shestakov et al. (2021) reported the TEC oscillations lasting for an hour during the 2009 eruption of the Sarychev Peak volcano, Kuril Islands, Russia.
The second type of disturbances occur 8-10 min after volcanic explosions by short pulses of acoustic waves propagating upward from the surface to the ionospheric F region (Fig. 1). They make short-term N-shaped impulsive TEC responses as Heki (2006) observed with the GPS-TEC method after the Vulcanian explosive eruption of the Asama volcano, Central Japan, on September 1, 2004. Despite many reports of the first type of disturbances, this second type of disturbances have not been reported since the Asama eruption. In this study, we report four new examples of this type and compare their amplitudes together with the 2004 Asama case.
Intensity of a volcanic explosion has been studied by atmospheric pressure changes associated with the airwave (infrasound) generated by the eruption (e.g., Matoza et al. 2019). However, geometric settings of such sensors relative to volcanoes are diverse, and amplitudes of such airwaves are difficult to serve as a universal index to describe the explosion intensity. Volcanic explosivity index (VEI) is used to describe the intensity of the eruptions (Newhall and Self 1982). However, this index is determined by the amount of ejecta and does not directly indicate the explosion intensities. In this study, we explore the possibility to use the amplitude of ionospheric disturbance that occur ~ 10 min after a large explosion as the new index. For this purpose, we compare ionospheric TEC responses to five recent explosive volcanic eruptions of four volcanoes in Japan 2004-2015 comparing the GNSS-TEC data from GEONET (GNSS Earth Observation Network), a continuous GNSS network in Japan.

GNSS data
We calculated TEC by multiplying a certain factor to the phase difference of the microwave signals in two frequencies, L1 (~ 1.5 GHz) and L2 (~ 1.2 GHz), from GNSS receivers in the Japanese GEONET (here we use only GPS satellites). TEC indicates number of electrons integrated along the line-of-sight (LoS) connecting the ground stations and GNSS satellites. We often represent the point of observation using the latitude and longitude of ionospheric piercing point (IPP). It represents the point of intersection of LoS with the hypothetical thin layer within the F region at the altitude of the highest electron density, assumed as ~ 300 km based on routine observations of ionosonde (https:// wdc. nict. go. jp). We plot their surface projections, often called sub-ionospheric points (SIP), on the map.
We downloaded the raw GEONET GNSS data on the days of the eruptions from Geospatial Information Authority of Japan (https:// terras. gsi. go. jp). We use the phase differences between L1 and L2 microwave carrier phases and converted them to TEC. Such slant TEC (STEC) values are often converted to vertical TEC (VTEC) after removing inter-frequency biases in GNSS receivers and satellites. Here, however, we use STEC throughout the study to capture TEC signatures from LoS penetrating the wavefront with shallow angles (Fig. 1). We use VTEC calculated from Global Ionospheric Map (GIM; Mannucci et al. 1998) only to normalize the amplitudes of the STEC changes with the background VTEC. The details of the GNSS-TEC technique to study lithospheric phenomena are available in the book chapter by Heki (2021).
Interaction of the electron movement with geomagnetic fields allows us to observe such TEC disturbances from stations located to the south/north of volcanoes in mid-latitude regions of northern/southern hemispheres (Heki and Ping 2005;Heki 2006;Rolland et al. 2013;Kundu et al. 2021). Figure 2 shows an example of the STEC change time series before and after the 2011 February 11 eruption of the Shin-Moe volcano in Kyushu, SW Japan, observed at the station 0729 located in a small island to the south of Kyushu. We cannot take advantage of the dense network, because the volcano is located near the southern edge of SW Japan, and the TEC signatures are visible only at the southern side. We isolated the short-term signals by fitting the polynomials (with degree 7) to STEC time series and showing residuals from such reference curves. The degree of the polynomials is tuned so that the long-period fluctuations are effectively removed. The short-period TEC signatures made by volcanic eruptions are insensitive to the selection of the polynomial degree as demonstrated in Additional file 1: Fig. S1.
Small pulses occurring ~ 10 min after the eruption (Satellites 4, 10, 13) are caused by acoustic waves propagating from the volcano to the ionosphere. Incessant small fluctuations of TEC observed by Satellites 7,8,12,17,and 26 are natural variabilities of TEC intrinsic for low elevation satellites. The Dst index time series shows that geomagnetic activity is low on 2011 Feb. 11 and on the other eruption days (Additional file 1: Fig. S2). Additional file 1: Fig. S3 compares the TEC data in Fig. 2 with those on the previous (2011 Feb. 10) and the next (2011 Feb. 12) days of the eruption.

Five volcanic explosions of four volcanoes in Japan
We selected five recent explosive volcanic eruptions in Japan with clear TEC disturbance signals. The Asama volcano, central Japan, started eruptive activity at 11:02 UT on September 1, 2004, with a Vulcanian explosion associated with strong airwaves (Nakada et al. 2005). For this eruption, there have been reports of ionospheric disturbance using GNSS-TEC by Heki (2006) and by HF-Doppler measurements by Chonan et al. (2018). Plume height was unknown due to cloudy weather, and they detected the atmospheric pressure change exceeding 205 Pa at a sensor located ~ 8 km to the south (Yokoo et al. 2005).

Fig. 1
Ionospheric disturbance caused by continuous (Type 1 left) and explosive (Type 2 right) volcanic eruptions can be detected by differential ionospheric delays of microwave signals of two carrier frequencies (L1 and L2) from GNSS satellites. Strong continuous eruptions sometimes excite atmospheric modes and long-term oscillatory disturbances in ionosphere. For explosive eruptions, we often find short-term impulsive disturbances in ionosphere 8-10 min after eruptions, the time required for acoustic waves to reach the ionospheric F region. The acoustic wave makes electron density anomalies (pairs of positive and negative anomalies as shown with red and blue colors in the figure) on the southern side of the volcano (for northern hemisphere cases) due to interaction with geomagnetic field (blue arrow) VEI of this eruption is reported as 2 according to Global Volcanism Program (2013) (same source for VEIs of the other eruptions).
Sakurajima is an active stratovolcano in the Kagoshima Prefecture, Kyushu, and is one of the most active volcanoes in SW Japan. Since 2009, several hundreds of explosions occur every year in the volcano. There were 125 eruptions in 2009 October, in which 101 were explosive. The explosion of Minamidake, Sakurajima, at 07:45 UT on October 3, 2009, was one of the strongest eruptions in its activity since 2009. Plume reached the height ~ 3000 m above the caldera rim, and the atmospheric pressure change exceeding 294.5 Pa was detected at a sensor ~ 5 km southeastward. At another observatory, ~ 11 km to the west of the vent, pressure change of 74 Pa was observed (JMA 2010).
Shin-Moe Volcano is also located in the Kagoshima Prefecture, Kyushu. We studied two VEI 2 explosive eruptions in 2011. In the first eruption (Jan. 31 22:54 UT), the plume reached the height of ~ 2000 m above the caldera rim, and the atmospheric pressure change exceeding 458.5 Pa was observed at a sensor ~ 2.6 km southwestward. In the second eruption (Feb. 11, 02:36 UT), the plume reached the height of ~ 2500 m above the caldera rim, and the atmospheric pressure change exceeding 244.3 Pa was observed at the same sensor (JMA 2013).
The last volcanic eruption was Kuchinoerabu-jima volcano, located at a tiny island Kuchinoerabu-jima ~ 100 km to the south of Kyushu. A VEI 3 eruption occurred on 29 May 2015 (00:59 UT). The plume height was ~ 9000 m, and pyroclastic flow reached the ocean. An atmospheric pressure of 62.2 Pa was observed at a sensor located ~ 2.3 km northeastward (JMA 2015). Nakashima (2018) studied ionospheric disturbances caused by this eruption using 1 Hz high-rate GNSS data. Figure 3 compares the impulsive TEC changes with periods of 1-2 min, ~ 10 min after the explosion. They are all STEC and the time series show the residual from the best-fit polynomials with degrees 7-9. We also see faint harmonic oscillations (similar to Type 1 disturbance in Fig. 1) sometimes follow the eruptions, e.g., the 2015 Kuchinoerabu-jima eruption. Here, we focus on the N-shaped TEC disturbances.

Comparison of ionospheric disturbances by the five volcanic explosions
Strong disturbances can be seen only from GNSS stations located to the south of the volcano due to the interaction with geomagnetic fields (Heki 2006). Therefore, we could not fully take advantage of the dense GNSS network, because the Sakurajima, Shin-Moe and Kuchinoerabu-jima volcanoes are all located in southern Kyushu and GNSS stations are sparse to their south.
Following Heki (2006), we try to adjust a simple function made of a set of positive and negative pulses, to the disturbances observed by five eruptions (Fig. 4). This function has a maximum and minimum at t = − σ and t = σ, respectively (drawn as a smooth curve in red in Fig. 4). The two parameters a and σ representing the amplitude and period, respectively, were tuned to minimize the root-mean-squares (rms) of differences between the synthesized and the observed disturbances. We also move this function in time to find the optimal arrival time of the disturbance with the time step of 2.5 s. This procedure minimizes the influence from the temporally sparse data (30 s sampling). The values σ = 19.5 resulted in good fits to the majority of time series, which corresponds to 78 s (1.3 min) as the approximate period (i.e., 4 × σ) of the disturbance. This is consistent with the period inferred form the 1 Hz data for the Kuchinoerabu-jima eruption (Nakashima 2018). From the adjusted values of a, we obtained the peak-to-peak amplitude, i.e., f(− σ) − f(σ), as summarized in Fig. 5. Here we focus on the amplitudes of the  Fig. 3 for each eruption. The period of the TEC fluctuation (Heki 2006) is fixed to 1.3 min, and adjusted the amplitudes and time lags to minimize the difference between the model function and the observed TEC observed disturbances, and we do not discuss their propagation velocity (~ 0.8 km/s). The same calculation has been done for all the three examples from each eruption and the average peak-to-peak amplitudes are compared in Fig. 5.
Volcanic eruptions excite acoustic waves in neutral atmosphere layer and ionospheric electrons move together with such neutral atmospheric molecules. Naturally, the strength of the TEC disturbances is largely influenced by the electron density in the F region of the ionosphere. As the index for the explosion intensity, it will be reasonable to normalize the amplitudes of STEC changes with background electron densities in the F region. Here we used background VTEC to normalize such amplitudes and express the TEC amplitudes relative to them (dark blue squares in Fig. 5). VTEC values at the time and location of eruptions are obtained from GIM.
Another important factor is the geometric condition of line-of-sight with the acoustic wavefront. Additional file 1: Fig. S4 shows the propagation of an N-shaped acoustic pulse in the standard atmosphere (Kundu et al. 2021). Additional file 1: Fig. S5 compares the synthesized STEC waveforms at various distances from the volcano with various elevation angles of satellites. The amplitude is not so sensitive to the distance from the volcano for the distance range studied here. However, different elevation angles cause different incidence angles of line-of-sight with the wavefronts, and change the disturbance amplitudes up to a few tens of percents. We corrected for such geometric factors for the five eruption cases by converting to an arbitrary standard geometry and plotted the disturbance amplitudes after correction with light blue squares in Fig. 5.

Discussion
To compare intensities of volcanic explosions, we often use VEI. They are either 3 (2015 Kuchinoerabu) or 2 (other 4 eruptions) for those studied here. VEI does not have finer scales and is not useful to compare intensities of explosive volcanic eruptions of this class. Intensities of volcanic explosions can be studied also by measuring amplitudes of airwaves (atmospheric pressure changes). However, different distance of the ground sensors from the volcanoes and different topographic and vegetation conditions makes it difficult to compare such intensities for different volcanoes.
In Fig. 5, we compare atmospheric pressure changes by the airwaves for the January 31 and February 11 explosions of the Shin-Moe volcano detected using the same sensor at the YNN station (JMA 2013). These two eruptions show similar amplitudes of STEC changes. However, the background VTEC at the time of the February eruption was more than twice as strong as those in the January eruption. Hence, relative amplitude of the January eruption becomes more than twice as large as the  Blanc (1985) showing the attenuation of airwaves in the Earth's atmosphere. The frequency 12.8 mHz corresponds to the higher end of the atmospheric bandpass filter February eruption. After correcting for geometric factors, the ratio becomes slightly less than twofold. This is in agreement with the difference of the pressure changes for these two eruptions (458.5 Pa for the January 31 eruption, and 244.3 Pa for the February 11 eruption).
This suggests the validity of using the relative amplitudes of the ionospheric STEC changes as the new index to describe the intensities (explosion energies) of volcanic explosions for different volcanoes. Its benefit is that we do not rely on the deployment of infrasound sensors, i.e., we can use this index whenever permanent GNSS networks are available on the southern/northern side of the volcano in northern/southern hemisphere.
Its drawback is that this index can be used only for strong volcanic explosions occurring when number of ionospheric electrons are sufficient (e.g., during daytimes). In fact, there were two explosions of the Shin-Moe volcano (Feb. 1 20:25 UT, and Feb. 13 20:07 UT) with stronger airwaves than the February 11 02:36 UT eruption. However, we cannot find ionospheric disturbances for these explosions because of small background VTEC early in the morning (Additional file 1: Fig. S6).
We also looked for such TEC signatures outside Japan. However, we failed to add more cases due to the lack of GNSS stations in appropriate places or to the insufficient intensities of the explosions. The only exception is the TEC signatures made by the human-induced explosion in 2020 August in Lebanon (Kundu et al. 2021).
As shown in Fig. 3, TEC changes by the five different volcanic explosions have similar periods of ~ 1.3 min. Such a uniformity suggests its origin in the atmospheric structure rather than characteristics of the volcanic eruptions. In Fig. 6, we compare this period with the diagram of atmospheric attenuation of acoustic waves with various periods at different altitudes (Blanc 1985). There, 1.3 min corresponds to the shortest period of the airwaves that can reach the altitude of ionospheric F region (~ 300 km) without large attenuations, i.e., 12.8 mHz corresponds to the high end frequency of the atmospheric bandpass filter. Infrasound records observed at ground sensors associated with explosive volcanic eruptions have stronger powers in periods much shorter than 1.3 min (e.g., Matoza et al. 2019). However, only those with periods 1.3-4.0 min can reach the ionospheric F region without strong attenuation. Because the original spectrum had larger powers for higher frequencies, we would have detected the 12.8 mHz component as the TEC changes at the F region altitude.