Correlation Between GNSS‐TEC and Eruption Magnitude Supports the Use of Ionospheric Sensing to Complement Volcanic Hazard Assessment

Despite the global threat posed by large‐scale eruptions to communities, to the climate, and to the consequent impacts on the world economy, many active volcanoes still lack of adequate ground‐based instrumentation. Satellite‐based remote sensing has been used to complement volcano monitoring and risk assessment for volcanic ash, but this technique is often limited by weather conditions. In this work, we explore the ionospheric total electron content (TEC) perturbations measured by GNSS to provide additional information and complement conventional monitoring systems. To this end, we measure the GNSS TEC perturbation associated with the acoustic‐gravity waves generated by 22 volcanic explosions. We introduce a new metric—the Ionospheric Volcanic Power Index (IVPI)—to quantify the energy transferred to the ionosphere by volcanic explosions. We evaluate the IVPI against several well‐established metrics from seismic and infrasonic volcano monitoring as well as satellite remote sensing. Our results show that the IVPI successfully correlates with the Volcanic Explosivity Index (VEI) for events larger than VEI 2. Moreover, the IVPI shows strong correlation with both the acoustic source power and the ash plume height, from which depends the style of volcanic activity. Moderate correlation between IVPI and peak ground velocity (PGV) requires further study in order to evaluate the role of different parameters (seismic magnitude, attenuation, style of faulting, crustal structure, etc.). Our results suggest that ionospheric monitoring by GNSS TEC can help to characterize volcanic eruptions, opening new exciting avenues for continuous volcano monitoring and warning systems by remote sensing.

A common remote sensing technique for volcanic ash detection is the use of thermal infrared imagers on satellites (e.g., MODIS, AVHRR, GOES, and SEVIRI), which record: the abundance of fine-grained volcanic ash; the temperature contrast between the volcanic cloud and the underlying surface; and the relative abundance of water within the volcanic cloud. However, this technique is far from straightforward (Brenot et al., 2014), and is limited by ice clouds in cold regions, and by the high water-vapor quantity in the tropics (Pavolonis et al., 2006), as well as by ice-coating phenomenon affecting the ash particles (Taisne et al., 2019). Another method to approximate the position of volcanic ash plumes is based on the measurement of gases emitted by volcanoes, especially sulfur dioxide (SO 2 ), that is assumed to follow a similar atmospheric dispersion pattern than the ash. SO 2 is easily tracked due to its low background concentrations in the atmosphere and its absorption bands in different distinct spectral regions (Carn et al., 2016). However, this technique only yields useful information when the eruption is accompanied by a detectable amount of SO 2 .
In recent years, infrasound has become a common technique for volcano monitoring. The International Monitoring System (IMS) infrasound network operated by the Comprehensive Nuclear-Test-Ban Treaty Organization was originally designed to detect any violations of the treaty , but it has also demonstrated its capability to detect and locate geophysical signals generated by natural hazard events (e.g., Hedlin et al., 2012). Infrasound sensors have been used to detect volcanic explosions for which the signal propagation over large distances is controlled by atmospheric winds and temperature (e.g., Dabrowa et al., 2011;Fee et al., 2011;. Large explosive events, such as the eruptions of Kasatochi (2008) and Okmok (2008) volcanoes in Alaska, USA, have been detected by infrasound sensors located as far away as 4,400 and 5,200 km from the source, respectively . The 2005 eruption of Manam volcano, Papua New Guinea, was recorded in Madagascar, about 11,000 km from the vent (Dabrowa et al., 2011). Infrasound is a cloud cover independent remote monitoring technique for volcanic activity, which is especially important in regions like Southeast Asia which has near-constant cloud cover (Taisne et al., 2019). Nevertheless, the density of the infrasound networks is currently limited, and the detections rely on the geometry of stations relative to the source and atmospheric attenuation conditions. In addition, detections of signals from large volcanic eruptions can be hampered by the presence of microbarom signals, i.e., coherent acoustic background noise in the 0.1-1.0 Hz frequency band. Such signals are generated by the non-linear interaction of the ocean waves, bathymetry, and coastlines (De Carlo 2020; Ouden et al., 2020;Smets, 2018). Therefore, additional infrasound network deployments in volcanically active regions are necessary to fully exploit this technique.
In remote environments and where physical and meteorological conditions hinder the application of these techniques, new methodologies could complement the current set of tools available to quantify the power of volcanic events and potentially further mitigate the losses. A promising method is based on the observation that powerful blasts, including volcanic eruptions, strong earthquakes, or even nuclear explosions, generate acoustic-gravity waves that propagate upward in the atmosphere and perturb the ionospheric plasma density (Occhipinti, 2015, and reference therein). During the upward propagation the acoustic-gravity waves, which are strongly amplified by the decreasing density of the neutral atmosphere, cause the perturbations of the ionospheric plasma. The perturbations result in the variation of total electron content (TEC), which is measured by GNSS as the integral of electron density along the line-of-sight between the GNSS receiver and the satellite (see material and methods for more details). Perturbations in the ionospheric TEC have been observed for a number of large earthquakes (e.g., Cahyadi & Heki, 2013Calais & Minster, 1995;Mikumo & Watada, 2010;Occhipinti et al., 2013) and tsunami (e.g., Artru et al., 2005;Occhipinti et al., 2008Occhipinti et al., , 2011Occhipinti et al., , 2013Occhipinti et al., , 2006Rakoto et al., 2018;L. M. Rolland et al., 2010); the ionospheric signatures of those events have been detected (Occhipinti et al., 2013) and explored to estimate the tsunami risk (Manta et al., 2020), the seismic magnitude (Occhipinti et al., 2018), as well as the oceanic displacement during tsunami propagation off-shore (Rakoto et al., 2018). Recent studies have shown the possibility of detecting tsunami-driven gravity waves in real-time, opening toward the possibility to a future implementation of ionospheric TEC monitoring to the existing tsunami early warning systems (Savastano et al., 2017;Savastano et al., 2019). The analysis of TEC has also revealed perturbations after volcanic eruptions at Pinatubo Volcano (Cheng & Huang, 1992;Igarashi et al., 1994), Asama Volcano (Heki et al., 2006), Soufrière Hills Volcano (Dautermann et al., 2009), Calbuco Volcano (Shults et al., 2016), and Kelud Volcano (Nakashima et al., 2016). However, the relationship between the TEC perturbation induced by the volcanic acoustic gravity waves and the corresponding magnitude and explosiveness of the eruption has not yet been analyzed.
We have recently shown that the ratio between the maximum level of power spectral density of the TEC perturbation, produced by the tsunami genesis, and the mean background level (or Ionospheric Tsunami Power Index) quantifies the power of an earthquake, providing additional information on the volume of the maximum water displaced during the seismic rupture and proving the potential role of ionospheric monitoring in tsunami risk estimation (Manta et al., 2020). Following the same approach of Manta et al., (2020), in this work, we define the Ionospheric Volcanic Power Index (IVPI), evaluating the relationship between the TEC perturbation and the volcanic eruption characteristics. With this investigation, we are able to determine the feasibility of using ionospheric sounding techniques to infer the style of volcanic eruptions. We compare the TEC perturbations observed during 22 different eruptive events to different parameters commonly related to the magnitude of the eruption (reported VEI, plume height, reduced seismic amplitude, and acoustic power). We investigate events that occurred at 12 different volcanoes, located at different latitudes, to include the effects of the magnetic field on the TEC perturbations (Astafyeva, 2019;Occhipinti et al., 2008;L. M. Rolland et al., 2013). The results of this work provide insight into the feasibility of using ionospheric sounding techniques to complement current volcano monitoring systems. The potential to remotely gain additional information by GNSS-TEC observations can be especially valuable for isolated volcanoes where monitoring is limited or non-existent.  (Heki et al., 2006); Calbuco volcano 2015, Chile (Shults et al., 2016); Kelud volcano 2014, Indonesia, (Nakashima et al., 2016). To identify variations in the TEC data recorded from the stations operating at the time of each eruption, we compare the TEC signal of six quiet days prior to the eruptions to the signal recorded at its onset.

GNSS-TEC Data
The TEC is the total electron content measured by GNSS receivers, originally proposed by Mannucci et al. (1993) as a new method for monitoring the ionosphere. The TEC corresponds to the integral of the electron density along the line-of-sight between the receiver and the satellite measured in TEC Units (1 TECU = 10 16 el/m 2 ). Following the methodology described in Manta et al. (2020), we computed the TEC as: where L1 and L2 are the carrier phase, corrected for phase ambiguities using the program Ninja within the GPS-Inferred Positioning System and Orbit Analysis Simulation Software developed by the Jet Propulsion Laboratory (JPL), California Institute of Technology (Blewitt, 1990(Blewitt, , 2000Lichten & Border, 1987) and f 1 and f 2 are the high and low GPS frequencies, respectively.
We considered the area located around the height of the maximum of ionosphere ionization (the F2 layer, fixed at the altitude of 300 km), as the main contributor to the observed TEC variations. Therefore, the observed perturbations are visualized at the ionospheric pierce point (IPP), which represents the intersection of the line-of-sight between each GNSS satellite-receiver pair and the ionospheric layer of maximum ionization. The TEC time series for all the satellite-station pairs were filtered with bandpass finite impulse response Butterworth filter with a range of 2 mHz-10 mHz, to remove the contributions from daily ionospheric variabilities, satellite motions, and instrumental biases. Spectrograms are then calculated with 30min moving windows and 99% overlap.
We then introduced the Ionospheric Volcanic Power Index (IVPI), to quantify the intensity of the disturbances induced by the eruptions in the ionosphere. We highlight that Equation 2 is the same as the one described by Manta et al. (2020) with a different alias to imply the different area of application (volcanic explosion parameters estimation rather than tsunami risk estimation), and it is defined as: PSD MAX corresponds to the maximum level of power spectral density of the TEC recorded during the event, and MBL is the mean background level. MBL is defined as the average maximum level of PSD of the TEC recorded by the same satellite-station pair during a 2-h window, starting from the event onset time, of the 6 days preceding the event. The size of the time window has been chosen based on previous studies that show the presence of an important ionospheric response during the first 2 h of the eruptions (Dautermann et al., 2009;Heki et al., 2006;Nakashima et al., 2016;Shults et al., 2016). We averaged the maximum PSD calculated in the 6 days before the event by the total number of days to further limit the effect of perturbations unrelated to the event.
For each event, among all the available TEC time-series recorded by the different GNSS satellite-receiver pairs, we selected for further analysis those which meet the following criteria: (1) presence of continuous recording during the event; (2) IPPs distance from the source <1,000 km. If these criteria were met, we visually inspected the TEC time-series and chose those showing less background noise in the 6 days before the event, to limit the effect of meteorological perturbations.
Furthermore, we investigated the effects caused by different geometrical parameters on the amplitude of the TEC signal, such as: the distance between the volcano and the GNSS receiver; the elevation angle of the satellite-receiver raypath; the distance between the volcano and the IPP; as well as the incidence angle between the satellite-receiver line-of-sight and the gravito-acoustic wavefront propagation. To be able to compare the IVPI values measured for the different events analyzed in this work, we selected satellite-station pairs with similar geometric characteristics.
MANTA ET AL.

Explosiveness, Seismic, and Infrasonic Measures
We collected information from the literature on several parameters related to the intensity of the eruptions. In particular, we used, when available, the Volcanic Explosivity Index (VEI), reported plume height, seismic Peak Ground Velocity (PGV), and infrasonic acoustic power. VEI is determined by different factors, including the ejecta volume, plume height, eruption duration, and other qualitative observations (Newhall & Self, 1982). The events that we analyzed have VEI values in the range from 2 to 4 as reported by the Smithsonian Institution's Global Volcanism Program (GVP) (Global Volcanism Program, 2012).
Another parameter used to characterize the explosivity of eruptions is the plume height. It is measured by ground-based sensors, including radars (Donnadieu, 2012;Donnadieu et al., 2016;Schneider & Hoblitt, 2013), LiDAR (Marenco et al., 2011), lightning detectors (Cimarelli et al., 2016), satellite-based sensors and infrasound (Perttu Taisne et al., 2020;A. J. Prata & Grant, 2001). We used plume height values reported in the literature, obtained from direct measurements and modeling (see Table 1 for references). The plume height for the 2016 Sinabung event, Indonesia (SI03), is based on the report issued by the Volcanic Ash Advisory Center (VAAC) in Darwin.
Further, we collected and analyzed broadband seismic data available from GEOFON and IRIS data centers.
For each event, we analyzed the data coming from a single seismic station with distances ranging from 32 km to 330 km. Following Zobin (2012), we measured the amplitude of the PGV of the first arrival extracted from the vertical component of the sensors filtered between 0.5 and 8 Hz. PGV values were adjusted to the reference distance of 10 km from the source to take into account the attenuation caused by geometric spreading and the loss of energy due to the anelasticity of the medium (Aki & Richards, 2002). The amplitude of the waves-function at the epicentral distance, r, is expressed following Battaglia and Aki (2003): where A 0 is the source amplitude, f is the frequency, β is shear wave velocity, Q is the quality factor, and n is the exponential factor controlling the amplitude decay of the seismic signal. According to experimental data, n varies between about 0.3 and 3, depending on the type of seismic wave and distance range considered (Bormann et al., 2012). Considering the generally large distances of the recording stations from the sources analyzed in this study, we assume that the surface waves are the dominant features, consequently n = 0.5. We use f = 4.25 Hz which represents the mean value of the frequency range (0.5-8 Hz) to filter the data. We use the range of shear wave velocities between 2 and 3 km/s. We chose quality factors in the range of values obtained by Koyanagi et al. (1995) for Kilauea Volcano, and De Natale et al. (1987): For events measured at distance <70 km, we used Q = 50, which better represents the highly fractured and weakened volcanic edifice; and Q = 300 for the events recorded at distance greater than 70 km considering the presence of more competent rocks.
Therefore, the PGV measured at r was adjusted as: Finally, we analyzed the infrasound signal recorded at the single station for each eruptive event. For each event-station pair the source acoustic power E, was calculated following Dabrowa et al. (2011): where P is the acoustic pressure amplitude measured from the single sensor or calculated from the beamform obtained for the stations in array configuration. The source to sensor distance is denoted as r, ρ a is the atmospheric density, c a is the atmospheric sound speed, and t is the time window of the signal. This equation assumes energy attenuation through hemispherical spreading, and for the remote detections (r > 250 km) this is therefore the estimate and not the expected actual acoustic power of the event, as the complexity of atmospheric propagation has not been accounted for. Nevertheless, our estimate of acoustic power is used to compare events following Dabrowa et al. (2011). Acoustic power was calculated at the vent, and at the reduced distance of 1 km and 10 km from the vent. In the case of multiple stations, the station with the maximum acoustic power was used for comparison (Table S1).

Results and Discussion
We used the available GNSS datasets to determine whether ionospheric TEC could provide useful information to volcano observatories and to the VAACs for mitigation of volcanic ash hazards.

GNSS-TEC and IVPI
For each event, we measured the TEC for all the satellites orbiting close to the volcano and we analyzed TEC measurements in both time and frequency domains to better detect the perturbations.
In Figure 2, we report an explanatory example of the analysis of the Anak Krakatau flank collapse and eruption on the December 22, 2018 (KR01) which triggered a tsunami that caused widespread losses and damage to the nearby communities (Perttu, Caudron et al., 2020). We analyzed data from 23 GNSS stations (21 stations from SuGAr with 15-s temporal resolution shown in Figure 2a that satisfied the criteria for the analysis of IVPI, we used the ones that showed stronger TEC perturbations (PRN 3, 7, and 23, shown in Figure 2a) to visualize the oscillation pattern in the hodochrone (Figure 2b). Tsunami are known to produces atmospheric internal gravity waves that can reach the ionosphere approximately 40 min after the event (e.g., Occhipinti et al., 2013). For this reason, we applied a bandpass Butterworth filter between 1 and 10 mHz in order to preserve any low frequency signal that can be linked to gravity waves attributable to the tsunami triggered by the Anak Krakatau flank collapse. The time series and spectral analysis of the TEC signal recorded by satellite PRN3 with respect to station LNNG (Figures 2c and 2d) show the two different perturbations induced in the ionosphere, and more generally visible in the hodochrone (Figure 2b): the first one traveling at the horizontal speed of ∼1 km/s with origin time corresponding to the time of the flank collapse and a frequency content between 1 and 4 mHz; and the second one, propagating at the same speed (∼1 km/s) but with stronger in amplitude than the first one, appearing around 22:30 (local time) in conjunction with the second and stronger pulse and a frequency content below 2 mHz. The first perturbation coincides with the strong infrasound signal recorded throughout the region and is mostly related to the M w 5.1 earthquake (20:55 local time, NW-SE trending focal plane) recorded by the regional seismic network and triggered by flank collapse (GEOFON Program website, Perttu, Caudron et al., 2020b). According to Perttu, Caudron et al., (2020) the tsunami generated by the flank collapse reached the coast of Java, Indonesia at 21:30 (local time). Considering the estimated propagation velocity of 1 km/s and the first arrival of the second perturbation (22:30, local time), the stronger gravity component observed in the second pulse can be attributed to the presence of a sustained plume in the atmosphere, rather than with the genesis of the tsunami. Indeed, the rise of the eruptive plume, constituted by a mixture of hot gas and particles, can displace the atmospheric medium inducing the propagation of gravity-waves (De Angelis et al., 2011;Prata et al., 2020;Ripepe, De Angelis, Lacanna, & Voight, 2010).
We applied the same approach to analyze all the eruptions reported in Table 1. The analysis of the available GNSS data revealed 13 events out of 22 with clear signatures in the ionospheric TEC. For these events, Table 1 reports the measure of the IVPI index as well as other parameters linked to the size of the eruptions such as VEI, plume height, PVD, and acoustic power. Figure 3 shows an example of the spectral analysis of the TEC signal recorded during the week preceding the 2014 Kelud eruption, KL01 (Indonesia) used to calculate the IVPI. It is possible to see how, in comparison to the quiet days preceding the event, the eruption is marked by high power of the spectrum, with a frequency component ranging from 2 mHz up to 10 mHz. As reported by Nakashima et al. (2016), oscillations in the ionosphere are visible from 24:00 to 2:00 (local time).
Interestingly, it was possible to remotely detect the 2011 Cleveland volcano eruption, CL01 (Alaska, USA). In this region, many active volcanoes are not monitored by ground-based instrumentation, consequently the use of remote techniques is crucial. Regional-scale seismic and infrasound observations have been already reported in the literature (De Angelis et al., 2012), showing the potential of infrasound detection in this remote region to detect the volcanic explosion about 50 min after the beginning of eruptive activity. We have also shown that the ionospheric TEC is capable of detecting the event as early as 40 min after the eruption started (Figure 4b  Another interesting result is the detection of the 2008 Okmok eruption, Alaska (OK01). The observation of gravity waves induced by the rising plume was detected analyzing the ultra long-period (ULP) signal at OKFG seismic station (De Angelis et al., 2011).
For some of the events, the detection by ionospheric sounding has been limited due to different factors, resulting in low or negative IVPI values (see Table 1 and Figure 5). The strongest influences affecting the computation of TEC and consequent IVPI are: the observation geometry (e.g. the elevation angle of the satellite-receiver raypath, the distance of the GNSS receiver from the source, and the distance of the IPPs from the source), the wind direction, and the local geomagnetic field. The raypath geometry between the satellites and the GNSS receiver affects the quality of the detection as discussed by Cahyadi and Heki (2015), who observed that the strongest TEC signal is recorded when the IPP is located between the epicenter and the GNSS station so that the line-of-sight crosses mainly the positive part of the wavefront. Moreover, satellites at low elevation angles have a better detection capability (Occhipinti et al., 2013).    Table 1 for legend reference and terminology. IVPI, ionospheric volcanic power index; VEI, volcanic explosivity index.
where λ, φ, and h depict the geographic longitude, latitude, and terrestrial altitude respectively.
In some cases, the small number of GNSS stations drastically reduces the possibility to detect the event due to a lack of TEC observations in the area around the source. Indeed, the IPP should be inside a radius of 1,000 km in order to detect the propagation of ionospheric disturbances. Another parameter that may affect the detectably of volcanic ionospheric disturbances is the impact of horizontal atmospheric winds (Shults et al., 2016), which, however, has not being quantitatively investigated so far.
The complex dynamics of the eruption, in addition to an unfavorable observation configuration (limited network, network geometry, wind direction, and magnetic field), can further reduce the possibility of detecting ionospheric perturbations related to volcanic eruptions. The Tolbachik volcano eruption in 2012 is a clear example: the TEC analysis does not reveal the peak of intensity related to the onset of the event. The eruption at Tolbachik Volcano was less energetic at the onset, the November 27, 2012, and gained intensity between late 28th and early 29th November as showed by the detection at the IS44 array (Albert et al., 2015).

Relationship Between IVPI and VEI
The VEI, proposed by Newhall and Self (1982), is a practical and widely used metric for categorizing the scale of explosive eruptions and accounts for both eruption magnitude (volume) and intensity (plume height). For the event discussed in this paper, we used values of VEI as reported by the Smithsonian Institution's GVP. For some of the events analyzed, the VEI was not reported because the specific eruption is considered part of a sequence of explosions and is therefore excluded from Figure 5a. The relationship between this metric and the IVPI is shown in Figure 5a, where it is possible to notice a strong positive correlation (Pearson's r = 0.95; p < 0.001). Interestingly, the 3 events (AS01, AS02, and SI01) displaying a VEI = 2 show a negative IVPI corresponding to the detection of a weaker TEC perturbation compared to the average daily variability normally present in the ionosphere. The three eruptions are very different in terms of PGV ( Figure 5c) and all have small plume with maximum height of 5 km (Figure 5b).
Further, the analysis showed that all the most powerful events considered (VEI = 4) produced strong TEC perturbations apart from the Tolbachik Volcano. While this might be surprising at first glance, the low IVPI can be explained by considering the fact that the event was less energetic at the onset, as mentioned in the previous section. Indeed, also seismic data recorded at 17:15 (local time), shown in Figure 5c, suggests that the event was not energetic enough. This is also confirmed by infrasound data, revealing why the acoustic signals did not propagate at large distances (Fee & Matoza, 2013).
In optimal conditions, with favorable geometry of the GNSS network available and a sufficient number of satellite-receiver pairs, the IVPI metric correlates with the VEI: if IVPI < 4.5, it is more likely that the event has VEI = 2, while a IVPI between 4.5 and 10 might suggest a VEI = 3 event, and finally VEI = 4 when IVPI > 10. A more precise threshold estimate can be validated with additional observations, which will be the subject of future study. Figure 5b shows the relationship between the IVPI and the plume height (Pearson's r = 0.81, p < 0.001), where two groups of events can be distinguished: those with small associated plumes (height < 6 km), clustered at the left of the plot, and the remaining 10 events producing greater plumes and presenting a nicer linear correlation. The presence of these two different trends can be attributed to the fact that small explosive eruptions, accompanied by smaller plumes, generate acoustic-gravity waves characterized by higher frequency content. These are strongly attenuated before reaching ionospheric heights, causing only very small or no perturbation in the GNSS-TEC (Sutherland & Bass, 2004). Indeed, the low values of IVPI shown by the first group of events indicate that acoustic-gravity waves generated by explosions with small plume height are trapped in the lower atmosphere and do not reach the ionosphere. On the contrary, large explosive eruptions, generally accompanied by higher plumes, generate waves comprising a lower frequency range that are less affected by the atmospheric absorption of sound (Dabrowa et al., 2011). Therefore, for larger events, given the dependency of the VEI on the plume height, the IVPI can provide useful information on the explosive nature of the eruptions. However, among the second linear group, we observe that despite the considerable altitude reached by the ash emitted by the Etna eruptions (between 10 and 15 km of altitude), these events caused limited perturbations ( Figure 6) and two of these events display a negative IVPI.

Relationship Between IVPI and Plume Height
Mt. Etna volcano is one of the major degassing volcanoes in the world (Allard et al., 1991), and the events in this study are part of a rapid sequence of summit eruptive episodes observed in December 2015. The high-energy events, or paroxysms, were characterized by Strombolian activity followed by high lava fountaining and abundant tephra emission (Corsaro et al., 2017). Electromagnetic perturbations were also expected in the upper atmosphere, especially for the ET01 event, which was accompanied by lightning activity, linked to electrostatic discharges in the ash plume as the eruption columns can carry significant electrical charge (M. R. James et al., 1998;M. James et al., 2000;Mather & Harrison, 2006). In addition, this event happened at night, during which the optimal conditions for the detection of ionospheric disturbances are found (Rozhnoi et al., 2014). However, these events were all accompanied by low seismic activity, with the ET01 event characterized by lowest PGV (Figure 5c), suggesting that this event did not emit strong acoustic-gravity waves detectable in the ionosphere as compared to ET02. Moreover, the reason why the large plumes seen in the Etna events produced a small ionospheric signal can be linked to a different infrasound production mechanism compared to the other events. Previous studies have suggested that the relationship between plume height and acoustic energy, which leads to perturbation in the TEC signal, is dependent on different factors according to the eruptive style. Indeed, the Plinian and sub-Plinian eruptions are accompanied by greater acoustic energy linked not only to the volume of emitted material such as in the case of the Etna events, but also to the gas-thrust region of the volcanic plume.
From Figure 5b we can also notice that when considering each volcano separately, the IVPI positively correlates with the ash column heights. This behavior is observed at Mt Etna, Calbuco and Merapi: in the case of Calbuco, the two events analyzed were characterized by a plume height of 17 km (CA01) and 21 km (CA02) and we measured a IVPI of 17.3 and 21.6, respectively; for Merapi, the two plumes reached an altitude of 4.5 km (ME01) and 17 km (ME02), and the corresponding IVPI was 4.7 and 12 (see Table 1 A specific relationship between IVPI and ash plume seems to characterize the volcanic activity and further studies are warranted to identify the individual correlation between eruption plumes at different volcanic systems and to unveil the physical and chemical properties controlling the amplitude of ionospheric disturbances.

Relationship between IVPI and Peak Ground Velocity (PGV)
We evaluated whether the energy transferred by the events to the ionosphere, quantified using the IVPI, is proportional to the PGV measured for the vertical component of the broad-band seismic station, which is linked to the power of the eruptions (Nishimura et al., 2012). The relationship between IVPI and the log value of PGV is shown in Figure 5c. The regression analysis showed a moderate correlation between the two metrics in linear scale, as confirmed by Pearson's r = 0.60 (p < 0.004).
An inspection of the events that did not cause greater ionospheric perturbations than the disturbances related to normal fluctuations in TEC (IVPI ≤ 0), showed an unexpected trend for the three events recorded at Sinabung, an Indonesian stratovolcano awakened in 2010 after 1,200 years of quiescence. SI01 refers to one of the first activities recorded after the initial small phreatic eruption of August 28, 2010 that marked the beginning of his active state. Among the events analyzed at this volcano, SI01 has the greatest PGV, recorded during an ash explosion that followed the volcano-tectonic earthquake swarm , but the lowest IVPI. SI02 is part of the second episode of small explosions that began on September 15, 2013 which, however, is associated with the highest IVPI and the smallest PGV. The difference between the Sinabung events is that SI01 was mainly accompanied by deep volcano-tectonic earthquakes, that were absent during the SI02 event. The SI02 event was characterized by only shallow volcano-tectonic earthquakes . The greater IVPI of SI02 can be related to the evidence that TEC is especially influenced by shallow earthquakes, which lead to disturbances in the sensitive F2 layer of the ionosphere (Shah & Jin, 2015). It has also been shown that seismo-ionospheric disturbances related to a shallow hypocenter cause greater disturbances as compared to deep hypocenter earthquakes (Kon et al., 2011;Lin, 2013;Liu et al., 2006).
Another surprising result is that similar PGVs are accompanied by a large variation of IVPI. This is especially visible when comparing KL01, KR01, and SI03, which have an IVPI of 18.4, 4.4, and −3.4, respectively, despite a range of PGV between 0.7 and 1.2 mm/s. The relationship between the two metrics seems more complex compared to the other measurements discussed, as seismic data are affected by many parameters. The nature and amplitude of the ground motions at a certain distance from the epicenter depends on the tectonic regime, magnitude, magnitude-dependent attenuation, and differences in crustal structure, to cite a few. Further, a more rapid attenuation of amplitudes with distance is observed at volcanoes displaying shallow crustal seismicity (McVerry et al., 2006). For these reasons, a more precise relationship between seismic information recorded during volcanic eruptions and the IVPI could be established by taking into account more variables such as earthquakes magnitude, depth, duration, and wave profiles.

Relationship Between IVPI and Infrasound
Infrasound detection of volcanic activity is dependent not only on the magnitude of the acoustic source, and therefore the type of activity, but also on the atmospheric conditions and relative location between the source and infrasound sensors. Volcanic infrasound can be recorded with local networks (<15 km), regional stations or arrays (15-250 km), or remote arrays (>250 km) (Fee & Matoza, 2013 (Caudron et al., 2015).
Generally, eruptions that produce higher plumes also result in higher acoustic power and lower frequency infrasound (Dabrowa et al., 2011;. Figure 5d shows the correlation between IVPI and the acoustic power recorded at the source. With the publicly available infrasound data, we were able to analyze only 10 over 22 events. The data show proportionality between the maximum acoustic power recorded at the source, and the intensity of the ionospheric perturbation expressed in IVPI values. As discussed in the previous section, while small explosive events produce acoustic waves characterized by high frequency content (0.5-4 Hz, Fee, Graces et al., 2010, Tungurahua), large explosive events produce a broader range of frequencies. In particular, as the total acoustic energy increases, the lower limit of detectable frequencies decreases (Matoza et. at., 2019). Moreover, lower frequencies of infrasound (e.g. Okmok Volcano, 0.01-0.5 Hz) are observed in presence of large ash plumes typically associated with Plinian eruptions ).
According to Sutherland and Bass (2004), the attenuation of acoustic energy in the atmosphere varies with the height and is frequency dependent. The energy decreases approximately with frequency to the power of 2 (Fee & Matoza, 2013). The relationship between IVPI and the acoustic power showed in Figure The events with the largest IVPIs such as Calbuco (CA01 and CA02), Kelud (KL01), and Okmok (OK01) are also associated with the strongest infrasonic signal and are all Plinian and sub-Plinian eruptions (VEI = 4).
On the other side, the events with VEI < 3 clustered on the left lower side of the graph, and three of them display negative values of IVPI (AS02, ET01, and TO01).

Conclusions
In this work, we evaluated the possibility of exploiting remote sensing techniques based on ionospheric seismology to complement volcanic surveillance. To this end, we investigated the relationship between the IVPI, a new metric that captures the magnitude of the ionospheric perturbations, to several parameters that are commonly used to quantify the volcanic explosive energy. We found that our metric strongly correlates with the VEI and the ash plume height, allowing a quick risk assessment especially for those volcanoes where a well-developed ground-based network is not yet available, and other techniques cannot be applied (e.g. lack of instruments or clouds coverage). In the optimal situation, the acoustic-gravity waves generated by a volcanic explosion can be detected after 8 min from the event, allowing a quick estimation of the height of the plume. A deeper examination is required to understand the link between the Peak Ground Velocity (PGV) and the IVPI. The relationship between acoustic power measured from the well-established infrasound technique and IVPI supports the idea that the new index can complement volcano monitoring where infrasound stations are not available or regions where the wind noise limits their use. The results confirm the potential of the proposed tool to be applied to complement existing monitoring techniques. The proposed method showed an efficient applicability in remote areas such as the Aleutian volcanic arc, Alaska, and in tropical regions, such as the Sunda volcanic arc, Indonesia, where harsh weather conditions and almost persistent cloud coverage can limit the applicability of the existing satellite-based techniques. We will carry out further examination of the relationship between IVPI and other parameters used to quantify explosive energy, such as mass eruption rate and infrasound, to support the applicability of ionospheric monitoring to remotely access the magnitude of large explosive events.

Data Availability Statement
The datasets generated during and/or analyzed during the current study are available in the DR-NTU repository, https://doi.org/10.21979/N9/GIP5TJ