Investigation of Ionospheric Response to June 2009 Sarychev Peak Volcano Eruption

Global Navigation Satellite Systems have been extensively used to investigate the ionosphere response to various natural and man-made phenomena for the last three decades. However, ionospheric reaction to volcano eruptions is still insufficiently studied and understood. In this work we analyzed the ionospheric response to the 11–16 June 2009 VEI class 4 Sarychev Peak volcano eruption by using surrounding Russian and Japanese GPS networks. Prominent covolcanictotal electron content (TEC)ionospheric disturbances (CVIDs) with amplitudes and periods ranged between 0.03–0.15 TECU and 2.5–4.5 min were discovered for the three eruptive events occurred at 18:51 UT, 14 June; at 01:15 and 09:18 UT, 15 June 2009. The estimates of apparent CVIDs velocities vary within 700–1000 m/s in the far-field zone (300–900 km to the southwest from the volcano) and 1300–1800 m/s in close proximity toSarychev Peak. The characteristics of the observed TEC variations allow us to attribute them to acoustic mode. The south-southwestward direction is preferred for CVIDs propagation. We concluded that the ionospheric response to a volcano eruption is mainly determined by a ratio between explosion strength and background ionization level. Some evidence of secondary (F2-layer) CVIDs’ source eccentric location were obtained.


Introduction
The ionosphere is the layer of charged particles extending within 60-1000 km above the Earth's surface. Global Navigation Satellite Systems (GNSS) have been extensively used to investigate it for the last three decades. The main advantages of GNSS-sounding compared to other techniques are the high measurement accuracy of total electron content (TEC) relative variations [1] as well as a huge and quickly growing number of globally distributed GNSS-sites. According to numerous studies the ionosphere is sensitive to different man-made and natural phenomena such as rocket launches, explosions [2] (including nuclear tests [3] and asteroid atmospheric blasts [4]), earthquakes [5], tsunamis [6], tropical storms [7], volcano eruptions [8], etc. Contrary to other phenomena the ionosphere response to volcano eruptions is still insufficiently studied and understood mostly due to a limited number of investigated volcanic events and the absence of dense GNSS-arrays in the vicinity of many active volcanoes.
The June 1991 Mt. Pinatubo (15.12°N, 120.33°E, Philippines) eruption was the first volcanic event in which covolcanic ionospheric disturbances (CVIDs) were explored by using the first generation of GNSS-system−NNSS [9]. The CVIDs generated by the September 2004 Asama volcano (36.24°N, 138.31°E, Honshu Island, Japan) explosive eruption were detected and analyzed by using Japanese GPS Earth Observation Network (GEONET) [8]. The CVIDs caused by the July 2003 Soufrière Hill volcano (16.43°N, 62.11°W Montserrat, Lesser Antilles) and February 2014 Kelud volcano (07.72°S, 112.31°E East Java, Indonesia) eruptions were investigated by surrounding GNSS-networks in publications [10][11] and [12].The CVIDs invoked by the April 2015 eruptions of Calbuco volcano (41.33°S, 72.61°W southern Chile) were studied in reference [13]. The volcanic explosivity index (VEI) of the eruptions discussed above vary from 2 to 4. The amplitudes of the discovered CVIDs range from 0.03 to 0.50 TECU (1 TECU = 10 16 electrons/m 2 ). Their velocity estimates vary from a few hundred to 1200 m/s that fits well the gravity and acoustic mode of travelling ionospheric disturbances. Almost all of the described CVIDs arise within 10-20 min after a volcanic explosion and/or ash ejection. They last from a few tenths of minutes up to 2.5 h after the source event. Most CVIDs were detected at distances rangingfrom 100-200 to 700-800 km from a volcano. Some arguments for the relationship between the strength of eruption and relative amplitude of the detected CVIDs were provided in [13].
Mt. Pinatubo, Kelud and Asama volcanoes belong to the so-called "Ring of Fire"−tectonically active rim of the Pacific Ocean comprising 328 active volcanoes. Its northwestern part (the Commander Islands, the Kamchatka Peninsula and the Kuril Islands) is highly eruptive but it has never been investigated for CVIDs. By this work, we start to fill this gap by studying the ionospheric response to the June 2009 Sarychev Peak volcano (48.10°N, 153.20°E, Matua Island, the Kuril Islands Arc) explosive eruption using observation data of Russian and Japanese GNSS-networks.

June 2009 Sarychev Peak Volcano Eruption
Sarychev Peak stratovolcano (elevation 1446 m above the sea level) is situated on Matua Island. It belongs to the Kuril Islands Arc and is located as far as 400 km to the southwest from the Kamchatka Peninsula tip (Russian Federation) and 750 km to the northeast from Hokkaido Island (Japan) (Figure 1). Strong explosive eruption of Sarychev Peak volcano occurred during 11-16 June 2009 after 33 years of inactivity. Its scale is characterized by the VEI = 4 (https://volcano.si.edu/database/search_eruption_results.cfm (accessed on 25 December 2020).). The active phase of the eruption consisted of two series of explosions separated by a 13 h pause [14]. More than 20explosions of different strength and duration resulted in eruptive columns of 3-21 km height. Ash plumes propagated within a few hundred up to 800 km from the volcano. The eruption chronology was reconstructed primarily based on satellite data including observations from the International Space Station [14,15] and infrasonic measurements provided by the International Monitoring System infrasound network [16,17]. No seismic records of the eruption were obtained by the nearest seismic arrays due to remote location of Matua Island. The main characteristics of the largest eruptive events are summarized in Table 1. Starting epochs of individual explosions were determined by infrasonic data with accuracy of about 5 min [16]. According to the publication [14] the highest ash plume reached 21 km altitude above the sea level. It was ejected by the 13 June explosion. However, the work [17] showed that series of the subsequent explosions generated notably larger infrasonic signals (see events 7-10 in Table 1). The a priorianalysis of daily traveltime diagrams formed for each day between 09-19 June 2009 and subsequent inspection of the filtered vertical TEC-series (see the next section for the explanations) also showed that the group of explosions occurred at 18:51 UT on 14 June at 01:15 and 09:18 UT on 15 June produced the most intense CVIDs comparing other eruptive events. Thus in the next sections we will focus only on these events.  Table 1. The main characteristics of the largest eruptive events with explosion plume altitudes ≥10 km above the sea level. The start time and duration of each paroxysm are given according to infrasonic observations [17]. The eruptive column height was estimated by satellite data and analysis at SVERT − Sakhalin Volcanic Eruption Response Team (http://www.imgg.ru/en/svert (accessed on 25 December 2020)). Amplitude of infrasonic signals is given for the nearest global IMS infrasonic network site located in Petropavlovsk-Kamchatsky [17]. The background total electron content(TEC) values were calculated by using IRI2016 online reference ionosphere model (https://ccmc.gsfc.nasa.gov/modelweb/models/iri2016_vitmo.php (accessed on 25 December 2020)).

Data and Methods
Sarychev Peak eruption occurred under the quiet geomagnetic conditions (index Kp≤2), absence of solar flares and at a low level of the local seismic activity.
The data of four GNSS-arrays were used to study the ionospheric response to June 2009 Sarychev Peak eruption ( Figure 1): IGS-network (International GNSS Service, http://www.igs.org/, accessed on 25 December 2020); KurilNet (Kuril GPS Network) deployed under collaboration of a few Russian and USA research institutions [18]; KamNet (Kamchatka GPS Network) maintaining by the Kamchatka Branch of Geophysical Survey of the Russian Academy of Sciences [19] and GEONET (GPS Earth Observation Network) operating by Geospatial Authority of Japan [20]. We used 31 GEONET sites almost evenly distributed through the eastern part of Hokkaido Island. The most distant GPS station was located of about 900 km whereas the closest site was situated only 5.7 km from active vent of the volcano. In total, observations of 48 GNSS stations were used in this work. All observation sites were equipped by dual-frequency receivers and tracked only GPS satellites with 30 sec data sampling rate.
The main theoretical aspects of GNSS application to ionosphere sounding are widely known and can be found in a number of publications (see, for example, papers [21][22][23] and references in them). In this work we adopted a classical approach to TEC data analysis, assuming that all TEC variations occur in a thin F2-layer of maximum ionization [24] (the shell model [1]). The intersection of the line of sight (LOS) between a GPS satellite and a GPS-receiver with this layer is called an ionospheric piercing point (IPP). The IPP sequences were used to analyze spatiotemporal evolution of CVIDs propagation, to calculate distances for the traveltime diagrams and to estimate apparent CVIDs propagation velocities.Thus knowledge of the F2-layer altitude (hmax) is important for the subsequent TEC analysis.
The value of the F2-layer altitude demonstrates diurnal, seasonal and long-term variations caused by different sources (the Earth rotation, solar activity, different atmospheric and geophysical processes, etc. [25]). The F2-layer height changes within 250-400 km depending on region and geophysical conditions (latitude, observation local time, season, level of geomagnetic activity, etc.). It is problematic to obtain accurate estimate of instantaneous hmax value in a region under study. Usually some mean hmax value is used. The hmax=300 km is often adopted (see, for example, papers [8,[26][27]). The better approach is to use data of ionozonde(s) located as close as possible to an observation GNSS-network. In this work we accessed data of the Petropavlovsk-Kamchatsky ionosonde (52.97°N 158.24°E, http://www.wdcb.ru/stp/data/ionosphere/46501_petropavlovsk_kamchatsky, accessed on 25 December 2020) located about 630 km to the northeast from the volcano. O-mode (f = 0.834f0F2) reflection altitudes (hpF2-height) obtained by the Petropavlovsk-Kamchatsky ionozonde were used for the hmax approximation. The hpF2-height values exhibit notable diurnal changes. However, the ionozonde data do not provide estimates of hpF2-height errors. That is why we fixed hmax value equal to 280 km basedon linear approximation of daily hpF2-height estimates obtained by the ionosonde during 10-17 June 2009 (see Figure S1).
For each day within 09-19 June 2009 the TEC-SUITE Ver. 0.7.6 utility (http://www.gnss-lab.org/tec-suite.html, accessed on 25 December 2020) was used for reconstructionof slant relative TEC variations from the raw GPS data stored in RINEX format. The slant TEC series were converted to vertical TEC variations by using Klobuchar's mapping function [28]: where VTEC−the vertical TEC value, STEC−the slant TEC value obtained using GNSS-observations, Rz−the Earth's radius, hmax−the altitude of the maximum ionization F2-layer, θs−the GNSS-satellite elevation angle which can be easily calculated using the broadcast ephemeris of the GNSS-satellite and the GNSS-receiver position. The equation takes into consideration the Earth's sphericity.
For the further analysis we performed a band-pass filtering of the TEC series by applyinga moving average with a time window of 1-5min (3-17mHz). According to a number of our tests this interval provides the most optimal signal to noise ratio. The similar frequency band was used to study CVIDs in other works (e.g., [13]). The brief description of the moving average method and its application to TEC data smoothing and filtering is given in File S1.
According to special studies (see, for example, [1,29]) and a great number of researchworks investigating ionospheric response to different natural and man-made phenomena (see, for example, references in publication [22]) errors of differential TEC estimates obtained by using dual-frequency phase GNSS-measurements do not exceed 0.01-0.02 TECU. As can be seen from the following illustrations (see the next section) the typical noise level of our filtered TEC series do not exceed 0.02 TECU.
The first CVIDs were met by the satellite G20 at 19:01:30 UT in 215 km to the southwest from Sarychev Peak, i.e., 10.5 min after the beginning of the eructation (Figures 2 and 3a). The single train of quasi-periodic oscillations with the maximum amplitude of 0.09 TECU, with the period of 3 min and duration of about 20 min was detected by the group of satellite G20 IPPs located southward from the volcano ( Figure  2). As can be seen from     Oscillations of much longer duration were revealed by the satellite G23 to the southwest from Sarychev Peak. They lasted of about 74.5 min from 19:05 to 20:19.5 UT (Figures 2 and 3b). The LOS G23−SHIK passed approximately orthogonal to the CVIDs front. It crossed the area right over the volcano and provided the detailed disturbance profile within 480-74 km out of the volcano. The discovered disturbances are characterized by the oscillation period of 3 min and maximum amplitude of 0.07 TECU which is slightly less than the amplitude of the CVIDs detected by the satellite G20. The appropriate TEC variations are visible up to 650 km out of the volcano in the traveltime diagram ( Figure 4). However, their propagation speed cannot be determined from the diagram. The most probable reason of this problem could be associated with a large space length of the wave packet. TEC disturbance maps also show prominent ionospheric variations at 300-600 km out of the volcano (Figure 5d (Figures 2 and 3c). The amplitudes of these variations are significantly less contrary to disturbances observed in opposite directions. The LOS G04−KHAC met CVIDs with period of about 3 min at 19:02 UT as far as 170 km to the north from Sarychev Peak. The period of disturbances and their detection distance are close to analogous parameters of the first CVIDs captured by the satellite G20 in opposite direction with respect to the volcano. However, their amplitude is notably less and waveform shape is different (compare Figure 3a,c).
TEC series of four other satellites G02, 04, 13 and 25 also exhibit CVIDs signatures between ~20:02−20:45 UT at 600-900 km to the southwest from the volcano. Their amplitudes only barely exceed background TEC variations (the last box in Figure 3c). However, they can be clearly recognized in TEC disturbance maps (see Figure 5f). Their apparent propagation velocities are estimated as large as ~700-800 m/s. The general view of the CVIDs spreading can be obtained from Animation S1.

4.2.01:15. UT, 15 June Explosive Event
According to infrasonic data ( Table 1) the next explosive event occurred 6.4 h after the previous explosion. Satellites G09, 15, 18, 26-28 showed notable variations in filtered TEC time series obtained within more than 3 h after the eructation.
The first TEC oscillations most probably correlated with the volcanic activity were detected by LOS G15−KETC at 01:19:30 UT about 51 km to the south of the volcano, i.e., 4.5 min after the beginning of the eruption episode (Figures 6 and 7a). The period and maximum amplitude of these variations were equal to ~2.5 min and 0.06 TECU accordingly. The compact wave packets with maximum amplitude of 0.06-0.09 TECU and period of 2.5-3 min were detected by LOSs G15−MATC and G26−MATC at 75 and 228 km to the east and northwest from the volcano. They were discovered in 8.5 and 13.5 min after the first CVIDs detection (Figure 7a). The TEC pulse with amplitude of about 0.05-0.06 TECU was observed by LOS G09−PARM at 01:27:30 UT 100 km to the north of Sarychev Peak.   Figure 7b. The LOS G28−SHIK passed over the volcano from the southwest to northeast and captured intense TEC variations between 01:27-02:01 UT starting from 146 km ahead and ending by 105 km behind the volcano ( Figure 6). The duration of the most intense TEC oscillations was equal to 34 min which is close to the appropriate duration of infrasonic signals triggered by this eruptive event [17]. The maximum TEC perturbation amplitudes were observed at 70-80 km to the southwest from the volcano within 01:38:30-01:40 UT and quickly faded in accordance with northeastward movement of satellite G28. The appropriate ionospheric variations can be clearly seen in the traveltime diagram ( Figure  8). Covolcanic TEC oscillations with longer period of 4-4.5 min and amplitudes of about 0.03-0.04 TECU were detected by the satellites G09 and G27 after 03:00-03:30 UT near Kunashir and Iturup Islands (Figures 6 and 7c).
Analysis of TEC disturbance maps formed for the time interval 01:15-04:30 UT showed that the intense CVIDs were observed between 01:27-02:02 UT to the southwest and west within ~200 km of the volcano (Figure 9a-d). They propagated with apparent velocities of about 1000 m/s. Clear CVIDs became visible within 02:18-02:50 UT at 700-900 km southwestward from Sarychev Peak (Figure 9e,f). The last portion of relatively long-period CVIDs were recognized between 03:27-04:20 UT about 600-700 km in the same direction away the volcano (Figure 9g,h). The appropriate propagation velocities of CVIDs visible in the far-field zone ranged within 800-900 m/s. The movie illustrating CVIDs propagation is given by Animation S2.

4.3.09:18. UT, 15 June Explosive Event
The last eruptive event initiating significant CVIDs occurred 8 h after the preceded explosion. Only three satellites G21, 24 and 29 detected TEC variations which could be associated with the volcano activity. The LOSs of the satellites G24 and G29 scanned the region within 300 km to the southwest, west and northeast from Sarychev Peak ( Figure  10). The first CVIDs were detected by the LOS G29−SHIK in 10 min after the explosion at ~90 km to the south from the volcano. The LOS G24−SHIK caught the CVIDs at almost the same time and direction (Figure 11a,b). The maximum amplitudes of the recorded disturbances ranged between 0.06-0.09 TECU. The amplitudes of TEC variations observed by the satellite G24 were notably higher than signals detected by satellite G29. Period of the most intense oscillations was equal to 2-2.5 min whereas the following much lower amplitude variations were characterized by periods of about 3-3.5 min.  The satellite G21 IPPs tracks cover the region to the southwest of the volcano as far as 250-550 km from it. The chain consisting of two 3.5-4.0 min period wave packets was observed. The oscillations lasted of about 56 and 20 min correspondingly. Their amplitudes were notably lower comparing the CVIDs amplitudes detected by satellites G24 and G29 (Figure 11c). The traveltime diagrams also confirmed excitation of prominent TEC disturbances within 0-600 km from the volcano (Figure 12). TEC disturbance maps demonstrate appearance of the first CVIDs to the south and southwest from the volcano at 09:29:30 UT as far as 100 km from it (Figure 13a). The most intense variationsstarted at 09:33 UT and lasted until 09:59 UT (Figure 13b-d). They propagated with apparent velocities of about 1200-1300 m/s southwestward and westward from Sarychev Peak. At 10:23 UT the CVIDs reached the group of IPPs located at 400-500 km to the southwest from Sarychev Peak (Figure 13e,f). Their estimated propagation velocities ranged within 1000-1200 m/s approximately. The disturbances were being observed until ~11:15 UT. It should be noted that the intense background TEC variations were detected around eastern Hokkaido and its vicinities within 9:00-11:30 UT. These variations obviously were not correlated with CVIDs. However, they probably masked covolcanic perturbations which could be detected at 800-900 km from the volcano as they were observed for the previously considered eruption episodes.
The movie illustrating CVIDs propagation is given by Animation S3.

Discussion
The 11-16June 2009 eruption of Sarychev Peak volcano provided us a lot of data for CVIDs analysis due to its relatively large scale, long duration, significant number of explosive episodes during the eructation and availability of the GPS networks around the volcano. Analysis of filtered relative vertical TEC variations for the total eruption period showed that a number of paroxysms triggered detectable CVIDs but only three events formed significant perturbations with maximum amplitudes up to 0.09-0.15 TECU. It is notably less comparedwith other class 4 eruptions (see, for example, [12,13]). According to publication [13], the most probable explanation of this phenomenon could be related to differences in background ionospheric TEC conditions. Indeed, the calculated background TEC values listed in Table 1 Table 2in [13]). Comparison of strength of explosions and the background ionization level which is given in Table 1 allows us to conclude that their combination determines ionospheric response to a volcanic eruption. For example, the high explosion strength of 16:43 UT, 15 June paroxysm is neutralized by the low background TEC level. However, the 08:13 UT, 12 June and 09:27 UT, 13 June high ionization level was compensated by the small size of the explosions (see the appropriate infrasonic signal amplitudes in Table 1). Other factors such as LOSs configuration with respect to a volcano, the IPPs' remoteness from it also may play an important role in apparent TEC variation formation.
The observed waveforms, oscillation periods (2.5-4.5 min), estimated propagation velocities (700-1800 m/s) point out at acoustic nature of the detected CVIDs. Spectral analysis of the detected TEC signals leads to the same conclusion ( Figure 14). It shows that disturbances of a longer period are observed in the far-field zone compared with the volcano's vicinities. Clear azimuthal inhomogeneity of CVIDs propagation was observed for all the considered eruption episodes. In agreement with previous studies [22], the southwestward spread of ionospheric disturbances prevailed among other directions. In this azimuthal sector the prominent covolcanic TEC oscillations were registered up to 900 km from the volcano. However, in northern quadrants we didnot detect CVIDs further than ~200-300 km from Sarychev Peak. Figures 2-3, 6-7 and 10-11 demonstrate fast decay of disturbance amplitudes to the west from the Kuril Islands Arc. Figure S2 shows an example of such amplitude attenuation. Most probably, this effect could be associated with structure of the magnetic field (Inclination = 60.8°, Declination = −7.2° in the volcano vicinities according to EMM2009-2019 model, https://www.ngdc.noaa.gov/geomag/calculators/magcalc.shtml#igrfwmm, accessed on 25 December 2020) [22]. However, some directivity of CVIDs emission and other possible sources including ionospheric winds also cannot be excluded from consideration.
The observed southwestward propagation of CVIDs in close proximity of the volcano on 15 June during 09:40−09:45 UT raises a question about the location of a source of these perturbations (see Figure S3). Most of the previous studies assumed that the CVIDs source islocated straight over a volcano. However, our observations contradict this assumption and fit better to eccentric position of a virtual ionospheric source which should be placed to the northeast from Sarychev Peak. This assumption may also explain a registration of intense positive phase signals to the northeast from the volcano. It does not match modeling results given in [22] and [30]. These studies stated that only negative TEC signal part can be observed to the north from a source in the northern hemisphere. A possible physical explanation of an eccentric location of the ionospheric CVIDs source may be associated with the upper atmosphere horizontal winds [13]. These flows may notably shift a virtual ionospheric source position at an altitude of the F2-layer. However, further investigations should be undertaken in order to verify this suggestion.

Conclusions
In this work, we analyzed ionospheric response to the 11-16June 2009 Sarychev Peak volcano eruption. Three explosive episodes occurred at 18:51 UT on 14 June (Event 1), at 01:15 UT on 15 June (Event 2) and at 09:18 UT (Event 3) were considered. Prominent covolcanic quasi-periodic TEC variations with amplitudes and periods ranged between 0.03-0.15 TECU and 2.5-4.5 min were discovered. For two former volcanic events the notable TEC disturbances were observed more than 3 h after the explosions. Propagation speed estimates of detected CVIDs varied from 700-1000 m/s in the far-field zone (300-900 km from the volcano) to 1300-1800 m/s in close proximity to the volcano. The characteristics of the observed TEC variations obtained allow us to attribute them to the acoustic mode signals.
Duration of the most intense TEC oscillations triggered by the eruptive events seems to be correlated with the duration of infrasonic signals (see Table 1): Event 1-76 min, Event 2-34 min, Event 3-26 min. TEC signals amplitude does not directly depend on explosion strength and eruptive column height. It is constrained by a complex of factors including background ionization level, wave front−LOS relative configuration and other reasons.
The southwestward propagation of the CVIDs prevails among other directions. In this direction the prominent covolcanic TEC disturbances could be detected up to 900 km away from Sarychev Peak whereas in north azimuth quadrants they were discriminated only within 300 km.
In this study, we obtained some evidence of possible eccentric position of ionospheric CVIDs source (secondary F2-layer source). The limited observation data amount remains an open question for future investigations.
Quick ionospheric response to strong volcanic explosions (in this study: Event 1−10.5 min, Event 2−4.5 min, Event 3−10 min) opens promising prospects for development of volcanic eruption early warning systems based on real-time TEC data provided by GNSS networks. The equatorward CVIDs propagation in northern hemisphere, configuration of GNSS orbital segments should be taken into account to design optimal LOS and IPPs coverage around a monitored volcano. Additional investigations should be undertaken in this direction in the near future.