SUMMARY
Continuous Plinian eruptions often excite atmospheric modes of ∼3.7 and ∼4.4 mHz, which are observed as harmonic oscillations of ionospheric total electron content (TEC) by global navigation satellite system (GNSS) receivers. Such TEC oscillations started shortly after the great eruption of the Hunga Tonga–Hunga Ha'apai (HTHH) submarine volcano at ∼4:14 UT, on 2022 January 15. Here I analyse GNSS data at stations within ∼4000 km from the volcano to study temporal and spatial distribution of such atmospheric modes. Strong ∼3.7 mHz TEC oscillations in near fields started shortly after the eruption onset and propagated outward with the sound speed from HTHH. Later such TEC oscillations became strong again with the amplitude peak at the distance ∼1400 km from HTHH. Such far field oscillations occurred also above New Zealand and the Solomon Islands, ∼3000 km from HTHH. Their amplitudes seem correlated with those of the 0S29 solid earth mode, suggesting that vertical surface vibrations underneath may play a role in maintaining the atmospheric mode. Onset of the far field TEC oscillations are synchronized with the local sunrise, possibly controlled by diurnal changes in the ionospheric electron density.
1 INTRODUCTION
Large volcanic eruptions often disturb ionosphere. They are observed as changes in ionospheric total electron contents (TECs) shortly after such eruptions using dual-frequency global navigation satellite system (GNSS) receivers deployed worldwide over the ground. Ionospheric disturbances associated with volcanic eruptions have two types. The first one is characterized by short N-shape disturbances, caused by sound wave pulses excited by Vulcanian-type explosions (Heki 2006; Cahyadi et al.2021).
The other type is long-lasting harmonic oscillations of TEC associated with atmospheric modes excited by continuous Plinian eruptions (e.g. Heki & Fujimoto 2022). Periodic surface pressure changes by strong atmospheric modes generate Rayleigh surface waves. They propagate globally and excite a few spheroidal modes of the solid earth free oscillation having frequencies close to the atmospheric modes. Seismometers recorded them as bichromatic ground vibrations having frequencies of ∼3.7 mHz (fundamental atmospheric mode) and ∼4.4 mHz (first overtone), after the 1991 Pinatubo eruptions (Kanamori & Mori 1992; Widmer & Zurn 1992).
Since GNSS-TEC observations started in 1990 s, there have been reports of harmonic TEC oscillations by strong volcanic eruptions, for example the 2014 eruption of Kelud volcano, Java Island (Nakashima et al. 2016), the 2015 eruption of the Calbuco volcano, Chile (Shults et al. 2016) and the 2021 Fukutoku-Okanoba eruption, Japan (Heki & Fujimoto 2022). However, none of these eruptions were strong enough to excite solid earth free oscillations.
The eruption of submarine volcano Hunga Tonga–Hunga Ha'apai (HTHH) at 04:14 UT, 2022 January 15, was the largest eruption since significant numbers of ground GNSS stations were deployed worldwide. Atmospheric disturbances of various kinds have been recorded by barometers, seismometers and by GNSS stations in various frequency bands (Matoza et al.2022). Ringler et al. (2023) reported that this eruption excited the solid earth modes. The intensity of the 0S29 (3.7 mHz) recorded by worldwide deployed seismometers was 11 times as strong as the Pinatubo case. They also showed that barometers near HTHH recorded ∼3.7 mHz oscillations of the surface pressure.
There have been numerous reports of ionospheric disturbances caused by the 2022 HTHH eruption. They report, for example multiple initial impulsive disturbances caused by explosions (Astafyeva et al.2022), worldwide propagation of ionospheric disturbances caused by Lamb waves (Themens et al.2022; Chum et al.2023) and detailed properties of its propagation (Ajith et al.2022; Kong et al.2022; Chen et al.2023; Li et al.2023). The ionospheric signatures of multiple passages of Lamb wave are recorded by dense GNSS networks in Japan (Heki 2022), and Indonesia (Muafiry et al.2022). Associated phenomena include simultaneous emergences of anomalies in both hemispheres connected by geomagnetic fields (Lin et al.2022; Shinbori et al.2022), and large-scale electron depletions (He et al.2022). Tonegawa & Fukao (2023) found ∼3.6 mHz oscillations of ocean floor and surface pressure around the Japan Trench, due possibly to the resonance between the ocean wave and the atmosphere.
In this paper, I study GNSS-TEC data of atmospheric resonant oscillation excited by this eruption, observed within 4000 km from HTHH. Muafiry et al. (2023) analysed data from the dense network in New Zealand and found atmospheric mode became active in New Zealand ∼14 hr after the eruption. Despite large number of past studies, such harmonic oscillations of TEC have not been sufficiently addressed so far. Here I expand Muafiry et al. (2023) using >10 GNSS stations deployed in the southern Pacific Ocean, to clarify the spatial and temporal distribution of the TEC oscillations, and their link to solid earth free oscillations.
2 ATMOSPHERIC MODES BY THE HTHH ERUPTION OBSERVED BY GNSS-TEC
2.1 GNSS data
I downloaded raw GNSS data of stations in the Pacific Basin from NASA Crustal Dynamics Data Information System (CDDIS; Fig. 1a) and converted the carrier phase data into time-series of slant TEC (STEC), number of electrons along the line-of-sight (LOS) connecting the ground stations and GNSS satellites (see Heki 2021 for details of the technique). The GNSS data include those from Global Positioning System (GPS), Russian GLONASS, European Galileo, Chinese Beidou Satellite System (BDS) and the Japanese Quasi-Zenith Satellite System (QZSS). BDS and QZSS include four geostationary satellites, and QZSS include three satellites having quasi-zenith orbits. These satellites enable us to observe long continuous TEC records that are quite useful for ionospheric volcanology (Heki 2022; Heki & Fujimoto 2022).
(a) Twelve GNSS stations (large red dots) around HTHH (star) analysed to study atmospheric modes. Data from two stations with smaller symbols (FTNA and TONG) were not suitable for the analyses in this study. Dark concentric circles denote distance from HTHH with the unit of wavelength of 0S29 solid earth mode. Atmospheric modes are recorded by barometers at MSVF and AFI (Ringler et al.2023). (b) LOS from a station close to the volcano (black lines) penetrates both positive (red) and negative (blue) electron density anomalies and is not suitable for observing them as TEC changes. On the other hand, LOS from a station ∼800 km to the north (e.g. SAMO; white line) is suitable for the purpose. Modified from fig. 6 of Heki & Fujimoto (2022). Geomagnetic inclination (B) is assumed as −35°. The electron density anomalies are meaningful only in a relative sense.
The data are available at TONG, the nearest station to HTHH. However, they are not used for analysing atmospheric modes because the LOS from a station close to the volcano penetrate both the positive and negative electron density anomalies in ionosphere and is not suitable for observing the disturbances as TEC (Fig. 1b). We also excluded the data from FTNA because of frequent data loss. The STEC data were converted to vertical TEC (VTEC) by removing biases so that the VTEC align with those calculated from global ionospheric map (GIM) models downloaded from CDDIS.
2.2 TEC oscillations at SAMO
Fig. 2(a) shows the VTEC time-series at SAMO (Samoa), ∼800 km north–northeast of HTHH, over the entire eruption day in UT. In addition to strong disturbances starting tens of minutes after the eruption, VTEC plots and their close-ups clearly show the occurrence of harmonic oscillations with a frequency ∼3.7 mHz (period: ∼270 s). We calculated the amplitude of this oscillation and show them in Fig. 2(b) with colours. Fig. 3 demonstrates the method to infer the oscillation amplitudes for synthesized data. I first perform wavelet transformation using the Gabor mother function having 3.7 mHz frequency (Fig. 3 middle). Then I obtain the smooth curve for its time-variable amplitude by initially doing the wavelet transformation twice with ±π/8 radian phase shifts, and then by adding their squares and calculating their square roots (Fig. 3, bottom).
(a) VTEC on 2022 January 15 at SAMO, observed using GPS, Galileo, GLONASS and BDS. Inset shows the close-up of 6–10 UT for GPS and Galileo. Satellites showing strong oscillation signals are drawn in red (GPS) and blue (Galileo). 1 TEC unit (TECU) corresponds to 1016 electrons m−2. (b) Strength of 3.7 mHz (period: 270 s) harmonic oscillation component as a function of time and distance from the volcano. The vertical dashed line indicates the eruption time (4:14 UT).
An example of synthetic 3.7 mHz oscillation signals with an amplitude of 1 TECU (top panel), and its wavelet transformation (middle panel) and the amplitude of the mode (bottom panel) recovered by combining two wavelet-transformed values with a time lag of 1/4 of the oscillation period. I used the Gabor mother function for wavelet transformation (red curve).
We notice that SAMO showed an activation of the oscillation shortly after the onset of the eruption ∼4:14 UT. It seems to extend as far as ∼1000 km in distance from HTHH. Here we call them the ‘near field’ activity. This activity continued until ∼11 UT. These near field TEC oscillation signals propagate outward from HTHH by ∼0.8 km s−1, the sound velocity in the F region of the ionosphere (see Section 3.2).
This is consistent with past cases, for example the 2014 Kelud and the 2021 Fukutoku-Okanoba eruptions. It should be noted that this is not the horizontal propagation of the atmospheric mode itself. The fundamental atmospheric mode is essentially the acoustic wave trapped between the surface and the mesopause. Although its horizontal propagation is quite slow, atmospheric waves leaking across the mesopause may propagate as acoustic waves in the ionospheric F region by ∼0.8 km s–1 outward. The near field TEC oscillations in Fig. 2 would have been generated in this way.
After this near field TEC oscillation ended, a weaker oscillation seems to reactivate at ∼18 UT at distances exceeding 1000 km. This is a part of the far field TEC oscillation to be discussed in the next section.
Barometers also record atmospheric mode as 3.7 mHz oscillation of air pressure at MSVF, Fiji, and possibly at AFI, Samoa (Ringler et al.2023). They started ∼40 min after the eruption onset. Tiwari et al. (2024) report that the smaller scale eruptions of HTHH on January 13 and 14 also caused similar atmospheric modes from GNSS-TEC observations at SAMO. They were not strong enough to excite the solid earth mode. In this paper, I concentrate on the data on January 15.
2.3 TEC oscillations at TUVA
Fig. 4 indicates the VTEC time-series at TUVA (Tuval), ∼1500 km north–northwest of HTHH, over the eruption day. The TEC oscillation during the first few hours can be seen also at TUVA with a few satellites having the ground projection of the ionospheric piercing points (IPP) of LOS ∼1000 km from HTHH. This would be the similar near field activity to those observed at SAMO. Fig. 4 shows that this near field activity decayed by ∼11 UT. After that, the oscillation became stronger twice, the first 12–14 UT and the second 18–22 UT. They have peaks at distances of ∼1400 km from HTHH. It is interesting to note that they are weaker at closer distances. Here I call them ‘far field’ activities because of its isolated strength peaks in space.
(a) VTEC on 2022 January 15 at TUVA, observed using GPS, Galileo, GLONASS, BDS and QZSS. Inset shows the close-up of 17–21 UT for BDS and QZSS. Three satellites showing strong oscillation signals are drawn in green (BDS) and black (QZSS). (b) Strength of 3.7 mHz (period: 270 s) harmonic oscillation component as a function of time and distance from the volcano. Far field atmospheric mode signals are clear for distances 1200–1600 km around 12–14 and 18–22 UT.
2.4 Frequency components
Fig. 5 compares the close-up of detrended time-series of the TEC oscillations, and their frequency spectra obtained by the Blackman–Tukey method for the near field (Figs 5a and b 6–10 UT, at SAMO) and the far field (Figs 5c and d, 18–22 UT, at TUVA) activities. For the former, we used Galileo satellites with conventional orbits with half-day orbital periods (E03, E15 and E08) to obtain data with IPP close to HTHH, say within ∼500 km in horizontal distance. For the latter, we used BDS and QZSS satellites having geostationary (C01, C04) and quasi-zenith (J02) orbits.
Comparison of the frequency components for the TEC oscillations shortly after the eruption observed at SAMO using three Galileo satellites (a), and those reactivated later observed at TUVA using three BDS and QZSS satellites (c). In both cases the 3.7 mHz component is prominent while the 4.4 mHz component cannot be identified (b, d). Dashed lines in (b) and (d) indicate the second and third harmonics of the 3.7 mHz oscillation.
In both the near field and far field activities, the 3.7 mHz peak is clear for the three satellites. However, it is not clear if the 4.4 mHz peak, the first overtone, exists. This is consistent with the seismometer records showing much smaller amplitudes of 4.4 mHz (Ringler et al.2023). This makes a strong contrast to the 1991 Pinatubo eruption. Both 3.7 and 4.4 mHz ground vibrations were observed in that case, and the latter had a larger amplitude (Kanamori & Mori 1992; Widmer & Zürn 1992). Such a difference in the major excitation frequency might reflect the plume height of the eruption and different energy distributions of the two modes in altitudes (Watada & Yoshida 2023).
The frequency spectra for the far field TEC oscillation 18–22 UT (Fig. 5d) show sharper peak at 3.7 mHz than those for the near field oscillation (Fig. 5b). This would reflect three factors. The first is the longer arc enabled by geostationary and quasi-zenith orbits of the GNSS satellites. The second factor would be the slower movement of IPP of the LOS connecting the ground stations and these geostationary/quasi-zenith satellites. The third factor would be the absence of clear propagation of the far field TEC oscillation, as discussed later in Section 3.2. The second and the third factor reduces frequency shifts of the observed TEC oscillation due to the movements of LOS relative to wave fronts.
3 DISCUSSIONS
3.1 Activated region of the far field TEC oscillations
Fig. 6(a) shows geographic distribution of the intensity of the 3.7 mHz TEC oscillation (positions show the IPP of LOS at 300 km) during the period 17–21 UT. In addition to the two stations given in Fig. 2 (SAMO) and 4 (TUVA), I show GNSS stations shown in Fig. 1. There is a dense network in New Zealand recording the oscillations (Muafiry et al.2023), but I here let three stations represent that region. Fig. 6(a) shows that large amplitude oscillations do not occur everywhere, that is they occur within several isolated regions.
(a) Intensity of 3.7 mHz TEC oscillation during 17–21 UT observed at 12 GNSS stations using all available satellites. They are plotted on their IPP tracks assuming thin ionosphere at 300 km altitude. Strong atmospheric mode signals tend to occur around the intersections of concentric circles where 0S29 has maximum amplitudes (0.5λ, 1.0λ, 1.5λ, 2.0λ) and the belt of strong ionization (average position of EIA, shown as the blue curve). (b) compares the average amplitudes of the oscillations at different distances from HTHH using 100 km bins. Error bars are 10 times as large as the standard deviation of the averages of the data within individual bins.
Those around TUVA, and SOLO have relatively high amplitudes for their distances from HTHH. It is interesting to note that these two stations are located just beneath the equatorial ionospheric anomaly (EIA). It is two belts of strong ionization located to the north and south of the geomagnetic equator (around ±15° in geomagnetic latitude), and Fig. 6(a) gives the average position of the southern part (Immel et al.2006). Smaller but significant oscillations occur for limited satellites observed at PTVL, AUCK, CKIS and THTG. Some of these stations are nearly 3000 km apart from HTHH. Considering that surface barometers recorded significant atmospheric mode only within 1000 km of HTHH (Ringler et al.2023), the GNSS-TEC technique may have higher sensitivity to atmospheric modes than ground barometers.
The significant TEC oscillations 17–21 UT also seems to be distributed on concentric circles drawn with solid (broken) curves indicating the distances 1, 2 and 3 times (1.5, 2.5 and 3.5 times) as long as the wavelengths of the solid earth mode 0S29 (Fig. 6a). Along these rings the amplitudes of the vertical crustal oscillations appear to have maxima (phases are reversed between 1λ and 1.5λ rings). Fig. 6(b) suggests that amplitudes of the TEC oscillations tend to have peaks at the intersections of EIA with the 1λ ring, and possibly at the 1.5λ and 2λ rings.
The atmospheric mode signatures propagated slowly from HTHH would weaken as they propagate a long distance. Indeed, surface barometers record the atmospheric modes only within 1000 km from HTHH (Ringler et al.2023). Fig. 6(b) may indicate that the local resonance between the ground motion and the atmosphere is essential to keep the atmospheric mode above the detectable amplitude in far field.
IGS stations in eastern Australia (TOW2) do not show any 3.7 mHz oscillations. TEC oscillation signals are not seen, either, from GNSS stations in Africa and Europe close to the HTHH antipode, where the ground vibrations by 0S29 are relatively large. This indicates that the ground vibration alone does not drive the fundamental atmospheric mode but suggests that the atmospheric modes propagated from the volcano need to be fuelled by the ground vibration to maintain their amplitudes.
3.2 Propagation of the TEC oscillations
As found in past cases, harmonic TEC oscillation signals caused by atmospheric modes by Plinian volcanic eruptions are observed to propagate outward from volcanoes at the acoustic wave speed in the F region of the ionosphere (∼0.8 km s−1, e.g. Nakashima et al. 2016). Fig. 7(a) shows that this is the case for the 2022 January 15 HTHH eruption. The time–distance plots with all the available station–satellite pairs for the near field atmospheric modes in 7–9 UT clearly show that the 3.7 mHz signals propagate outward from HTHH with the prescribed speed. This would be the propagation of the acoustic wave in the ionospheric F region that leaked from the atmospheric mode trapped beneath the mesopause, rather than the propagation of the atmospheric mode itself.
Time-distance plots showing the 3.7 mHz oscillation components extracted using a Gabor wavelet (Fig. 3) using all the available station–satellite pairs. For the near field data (left-hand panel), the oscillation signals propagate outward from the volcano by the acoustic wave speed (∼0.8 km s–1). On the other hand, for the far field data (right-hand panel), such a propagation is not clear. We can see propagation towards the volcano by ∼4 km s–1 in one region and ∼0.8 km s–1 in another region.
Fig. 7(b) shows the same plot in 18.5–20.5 UT for the far field TEC oscillations. There the signal propagation signatures are not simple, that is the phases are delayed for pairs closer to HTHH as if the signals propagate ‘towards’ HTHH. The propagation speeds also seem non-uniform, that is it looks acoustic wave speed (∼0.8 km s–1) in one region and the Rayleigh wave speed (∼4 km s–1) in another region. Muafiry et al. (2023) studied the propagation of the TEC oscillation signals using the dense GNSS network in New Zealand. They report propagation of the oscillation signals toward HTHH by ∼0.5 km s–1. I do not have a model that can explain the complicated propagation.
This complex situation of propagation makes a sharp contrast between the far field and near field TEC oscillations. It is at least clear that the far field atmospheric oscillations are not caused by acoustic waves propagated from HTHH through ionosphere. Instead, they might have been caused by the acoustic mode occurring right there, possibly by the resonance between the atmospheric mode propagated from HTHH and the vertical ground oscillations of 0S29 underneath. Various propagation velocities as found in Fig. 7(b) would reflect that the ground motion has not become pure standing waves but is still a mixture of Rayleigh waves propagating in multiple directions. Also factors such as wind and turbulence would influence the velocity fields as seen in Fig. 7(b).
3.3 The start and end of the far field TEC oscillations
Fig. 8 compares the time-series of the far field TEC oscillations, two satellites from each GNSS station, observed at various longitudes. It suggests that at eastern stations such as CKIS, the oscillation started before 17 UT. On the other hand, at SOLO, the westernmost site, it starts shortly before 20 UT. The time lag of ∼3 hr is roughly equivalent to that of the sunrise at the two points. Indeed, the far field TEC oscillations appear to have started shortly after the VTEC curves show positive bending by solar radiation increase associated with sunrise. This is clearly seen in Fig. 4(inset).
Comparison of the start times of the far field 3.7 mHz TEC oscillations at stations with different longitudes. The stations in the eastern and western longitudes are displayed in the upper and lower parts of the figure, respectively. Short vertical dashed lines given for individual stations indicate sunrise time at the GNSS stations on the eruption day. The oscillation seems to start earlier and later in eastern and western stations, respectively.
Fig. 8 also shows that the two stations TUVA and SOLO, located beneath EIA, show larger amplitudes of the TEC oscillations than other stations. Such amplitudes would reflect three factors, that is the intensity of the oscillation of neutral atmosphere (factor #1) and the degree of ionization (factor #2). The incidence angle of the LOS to the wave front is also an important factor (factor #3). The onset of the oscillation triggered by sunrise and larger amplitudes beneath EIA may suggest the importance of factor #2.
The geometry factor (factor #3) would be similar for stations close to local zenith. In Fig. 6(a), strong oscillation signals are seen by satellites near zenith at stations like TUVA and SOLO while they are not seen at stations like PTVL and NRMD. This contrast suggests that LOS geometry (factor #3) is not a dominant factor.
On the other hand, the importance of factor #1 is clear in Fig. 4. There the TEC oscillation signals at TUVA are strong during 12–14 UT, local nighttime when the VTEC remains low without solar radiation. In that case, factor #1 would have been strong enough to compensate for unfavourable factor #2. Fig. 4 also indicates that far field TEC oscillations decayed at ∼14 UT, remained invisible for ∼4 hr, then it resumed at ∼18 UT. The disappearance at ∼14 UT might be due to factor #1 (declining atmospheric mode intensity), and the restart at ∼18 UT might be due to factor #2 (sunrise).
Fig. 4 also shows that the TEC oscillation at TUVA ends by ∼22 UT although VTEC continued increasing toward the noon. Gradual decay of the eruption intensity is clear in the diminishing amplitudes of the solid earth mode 0S29 in seismometers (Ringler et al.2023). The disappearance of the TEC oscillation ∼22 UT would be therefore due to factor #1.
Far field TEC oscillations at TUVA are not seen before ∼12 UT, ∼8 hr after the eruption onset. This large time lag might reflect the time for the atmospheric mode to travel ∼1000 km from HTHH, suggesting that the group velocity of the atmospheric mode would not exceed tens of metres per second. Another possibility is that it took hours for stable atmospheric mode to develop around TUVA. Fig. 4(b) shows sporadic short-lived excitations of the atmospheric mode during 8–12 UT. The atmospheric mode may have propagated faster, having travelled ∼1000 km in a few hours. However, it may have taken significant time for the mature atmospheric mode coherent over a large region to grow.
4 CONCLUSIONS
Here I summarize the findings on the atmospheric modes caused by the eruption of HTHH on 2022 January 15, observed using GNSS-TEC technique at stations in southern Pacific Ocean. The monochromatic TEC oscillation of ∼3.7 mHz, corresponding to the frequency of the fundamental atmospheric mode, started shortly after the eruption and continued nearly until the end of the day. I found it appropriate to classify them into the early near field activity (0–1000 km from HTHH) and the later far field activity (1000–3000 km from HTHH). The former is stronger near the volcano and propagated outward at the sound wave velocity in the F region of the ionosphere. On the other hand, the latter appears to have peak intensities where the solid earth mode 0S29 have maxima and did not show outward propagation from the volcano. The far field atmospheric modes may need to be fuelled from the resonant vertical surface vibration underneath. Their start times synchronize with local sunrise. This suggests that the signal strengths of TEC oscillations are governed by the background VTEC as well as the atmospheric mode intensity.
DATA AND RESOURCES
The GNSS raw data and global ionospheric maps are downloaded from NASA/CDDIS (cddis.nasa.gov).
Acknowledgement
I thank two anonymous reviewers for constructive comments. I also thank Kiwamu Nishida (ERI, Univ. Tokyo), Kazunori Yoshizawa (Hokkaido Univ.), Mala Bagiya (Indian Inst. Geomagnetism), Ihsan Naufal Muafiry (ITB) and Shuanggen Jin (SHAO) for discussions. KH is supported by Chinese Academy of Sciences, President's International Fellowship Initiative (grant number 2022VEA0014).
CONFLICT OF INTEREST
The author declares that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
AUTHOR CONTRIBUTIONS
K. Heki: conceptualization, data curation, formal analysis, methodology and writing.
