Global Propagation of Ionospheric Disturbances Associated With the 2022 Tonga Volcanic Eruption

In this study, we use measurements from over 4,735 globally distributed Global Navigation Satellite System receivers to track the progression of traveling ionospheric disturbances (TIDs) associated with the 15 January 2022 Hunga Tonga‐Hunga Ha'apai submarine volcanic eruption. We identify two distinct Large Scale traveling ionospheric disturbances (LSTIDs) and several subsequent Medium Scale traveling ionospheric disturbances (MSTIDs) that propagate radially outward from the eruption site. Within 3,000 km of epicenter, LSTIDs of >1,600 km wavelengths are initially observed propagating at speeds of ∼950 and ∼555 ms−1, before substantial slowing to ∼600 and ∼390 ms−1, respectively. MSTIDs with speeds of 200–400 ms−1 are observed for 6 hrs following eruption, the first of which comprises the dominant global ionospheric response and coincides with the atmospheric surface pressure disturbance associated with the eruption. These are the first results demonstrating the global impact of the Tonga eruption on the ionospheric state.

• Highly directional Large Scale traveling ionospheric disturbances (LSTIDs) are generated from the Tonga Eruption at initial speeds up to 950 m/s and propagate potentially as far as 20,000 km • Outside the local area of the eruption, Medium Scale traveling ionospheric disturbances (MSTIDs) coincident with the surface pressure wave form the dominant ionospheric response • Secondary MSTIDs continue to be produced in the vicinity of the eruption for up to 6 hr following the event

Supporting Information:
Supporting Information may be found in the online version of this article.
Due to coverage limitations, most observational studies rarely examine the global impact of such events. We here use the global coverage GNSS observations to track ionospheric structures generated by the Hunga Tonga-Hunga Ha'apai (herein simply referred to as Tonga) that began at 04:15 UT 15 January 2022, centered at (20.546°S, 175.390°W) (https://earthquake.usgs.gov/earthquakes/eventpage/pt22015050/origin/detail). We will further examine and characterize the global ionospheric response to this eruption.

Data
GNSS TEC measurements have long represented a vital data set in the study of the global ionospheric response to a wealth of phenomena. Using measurements of carrier phase and code pseudorange on a pair of UHF signals, the integrated electron density of the ionosphere can be found (D. R. Themens et al., 2013). In this study, we are interested solely in the relative variations in TEC, which are typically calculated by applying a high-pass filter to un-leveled, un-bias-corrected slant TEC (sTEC) (Komjathy et al., 2012;Zhang et al., 2019). This approach, however, can lead to an elevation dependence in the amplitude of variations due to the slant propagation path. To mitigate this potential artifact, we instead fully calibrate our sTEC measurements as in D. R. Themens et al. (2013) by applying satellite bias estimates from the Center for Orbit Determination (CODE) (ftp://ftp.aiub.unibe.ch/) and then applying the Ma and Maruyama (2003) receiver bias estimation method to get a first-order estimate of the absolute TEC. By calibrating our sTEC, we can project the sTEC to vertical TEC (vTEC) using the standard thin ionospheric shell model with a shell height of 400 km without also imposing a significant projection function-dependent trend to the vTEC (D. R. Themens et al., 2015).
This vTEC is then detrended by removing the 30-min-wide, boxcar smoothed vTEC from the unsmoothed vTEC, which is consistent with what has been used previously in studies of the global propagation of Medium Scale traveling ionospheric disturbances (MSTIDs) and Large-Scale traveling ionospheric disturbances (LSTIDs) in Zhang et al. (2019); however, as noted in Coster et al. (2017), the length of the chosen detrending window can have an impact on the characteristics of the residual waves. Use a detrending window that is too wide and you risk introducing substantial trends from quiescent ionospheric variability and masking smaller scale structures behind stronger large-scale variability. Use a detrending window that is too narrow and you risk removing parts of the desired signal. In our case, we sit on the boundary of these two scenarios, where we have both LSTIDs and MSTIDs of interest. LSTIDs generally have propagation speeds between 400 and 1,000 ms −1 , wavelengths greater than 1000 km, and periods of 30-180 min, while MSTIDs tend to propagate at speeds of 250-1,000 ms −1 , have wavelengths of several hundred kilometers, and have periods of 15-60 min (Francis, 1975;Hunsucker, 1982;Zhang et al., 2019). Using the 30-min window will allow us to easily identify the MSTIDs, but may artificially suppress the observed amplitude of LSTID structures.
To further ensure the reliability of the data, any intervals with missing data, such as partial intervals at the start of arcs of lock or across cycle slip events are omitted to ensure the consistency of the smoothing interval and mitigate the impact of some cycle slips. Accounting for differential code biases (DCBs) may be an unnecessary step as the detrending process would remove slow moving trends imposed by erroneous DCB projections, but we still apply DCB corrections as a precautionary measure. To limit the amount of horizontal smearing by oblique ray paths, a lower elevation limit of 40° is used.
An example of the detrended vTEC at several stations in New Zealand (A), Australia (B), Japan (C), Eastern Canada (D), South Africa (E), and Northern Europe (F) is presented in Figure 1. To track the propagation of these waves globally, we have gathered data from 4,735 distributed GNSS receivers from dozens of local and global receiver networks. These stations correspond to all public data that were available within 2 days of the eruption. A map of the stations used is presented in Figure 1g, with the sample stations from Figures 1a-1f marked with crosses.
The aforementioned processing is applied to each receiver, the ionospheric pierce points (IPPs) of the TEC measurements are calculated using the rapid precise GPS orbits provided by the Crustal Dynamics Data Information System (CDDIS) (https://cddis.nasa.gov/archive/gnss/products/), and the resulting geolocated vTEC perturbations are mapped for the entire day of 15 January 2022. The resulting vTEC measurements are then plotted  on maps according to their IPP location for every minute of the day. Example videos of vTEC anomalies are provided in the supplementary material in Supporting Information S1.
To compare the propagation of the ionospheric waves to those at the surface, a barotropic version of the TIGAR (Transient Inertia Gravity and Rossby wave dynamics) model has been run at triangular truncation T170 horizontal resolution (i.e., 512 zonal points and 256 meridional points). TIGAR solves primitive equations on the sphere using the Hough harmonics (Vasylkevych & Žagar, 2021) thereby providing the time evolution of Rossby and inertia-gravity waves. The atmospheric response to the Tonga eruption is simulated using a homogeneous background (i.e., isothermal and with no winds). The eruption is represented by a Gaussian perturbation at the location of the eruption superimposed on the barotropic atmosphere with a mean depth of 10,114 km. The simulated horizontally-propagating gravity waves represent the Lamb wave in the hydrostatic, isothermal atmosphere. The mean depth and the parameters of the source were chosen to fit the observed horizontal structure of the eruption and the observed amplitude and timing of surface pressure perturbations.
One method to look at the global propagation of ionospheric disturbances induced by a known source is to organize the vTEC anomalies according to their distance from the source, as in Chen et al. (2017) to characterize earthquake-driven ionospheric TIDs. By binning all vTEC anomalies with latitudes between 55°S and 55°N in 50 km bins with respect to radial distance from the eruption location, we can isolate coherent structures associated with the event. We have focused here on regions below 55° since the high degree of other activity at the poles significantly complicates the identification of waves. This approach is applied in Figure 2a, where we see several strong vTEC anomaly structures fanning out from the eruption epicenter (y-axis) and from the time of eruption onset (x-axis). Note that the color scale has been saturated at 0.2 TECU to make wave features easier for the reader to see. Wave amplitudes should instead be inferred from individual station measurements to avoid attenuation/smearing from data binning. Also plotted in Figure 2a is a superposition of the radial propagation of the TIGAR-modeled Lamb wave pressure disturbance originating from the eruption. This pressure wave propagates at 315 ms −1 and is coherent with the generation and propagation of primary MSTIDs from the event.
To assess the speed of the disturbances propagating from the eruption event, we use two main methods. The first method is to use triangulation of the irregularities from neighboring GNSS raypaths. If the TEC signature of a TID is detected along at least three GPS satellite-receiver raypaths in a localized region, the relative detection time and corresponding IPP locations can be used to estimate the propagation velocity. This method, demonstrated in Figures 2b and 2c, assumes the disturbance has a uniform structure perpendicular to the direction of propagation and requires the assumption of an IPP altitude (400 km), if the altitude of the disturbance is unknown. This TEC triangulation technique has been previously applied in estimating the velocities of ionospheric structures associated with magnetospheric substorms (Watson et al., 2011) and ultralow frequency magnetic field pulsations (Watson et al., 2015).
The second method for characterizing the wave propagation speed is to trace the slope of coherent wave structures in plots of the distance from the epicenter versus. time, such as is presented in Figure 2a. One can then use the fitted slope for each independent structure as an estimate of the propagation speed. If a coherent wave structure exhibits a break in its slope, such as was seen in Figure 3a at 2,200-2,300 km before 06:00 UT, each section of the structure is treated separately. A demonstration of the result of both speed estimation methods is presented in Figures 3b and 3c for the regional domain around New Zealand. Using the derived speeds we can also determine an approximate wave generation time by extrapolating the fitted linear trends to zero range. This implied initiation time is used as the x-axis in Figure 3c, not to be confused with the times in Figures 3a and 3b. Similarly, the triangulated speeds for each wavefront can also be averaged and mapped backward to estimate their implied time of origin. Both approaches are done separately for ranges within and beyond 2,300 km to highlight the change in phase speed of the waves with radial distance from the epicenter.

Results and Discussion
In Figures 1a-1f, we see a number of waves arriving at sample GNSS receivers at roughly the time of the surface pressure anomaly, with the largest amplitudes generally coincident with the surface pressure wave; however, these plots are restrictive and can be challenging to associate directly with the wave fronts propagating from the Tonga eruption.
Looking at Figures 2a, 3a, and 3b, we can begin to characterize some of the wave structures emanating from the Tonga eruption. We first see an LSTID propagating radially at speeds of 950 ± 170 ms −1 , with a dominant wavelength >1,600 km and period of 48 ± 10 min, before slowing substantially to 600 ± 30 ms −1 after reaching ∼2,300 km range (Figure 3b). The trailing edge of this LSTID had a speed of 710 ± 115 ms −1 before the break at ∼2,300 km range, where it slows to 550 ± 15 ms −1 . A second LSTID with the same dominant  wavelength and period, then propagates from the eruption at 555 ± 45 ms −1 before similarly slowing substantially to 390 ± 15 ms −1 after ∼2500 km range. Wavelengths here have been determined by using a sliding least squares sinusoidal fit, similar to a wavelet transform, but for spatial data, and periods were determined using a wavelet transform at 2,500 km range. Examples of this analysis is provided in the supplementary material in Supporting Information S1.
Changes in propagation speed with radial distance can be clearly seen in Figure 3a, where the slopes of LSTID leading/trailing edges abruptly decrease. These breaks in propagation speed can be inferred more directly in Figure 3c by looking at the implied origin times of the waves. LSTIDs within 2,300 km of the eruption have origin times of ∼15 and ∼20 min after the eruption, respectively, while the portions of the LSTIDs beyond 2,300 km have implied origin times either coincident with the eruption or prior to it, depending on the velocity estimation method. A 15-min delay is consistent with the time it would take an acoustic wave to propagate into the thermosphere, while there is no indication that LSTIDs were generated at the origin prior to the eruption. The speed of the later portions of the propagation of these TIDs are consistent with horizontally propagating thermospheric gravity waves (e.g., the direct wave) (Mayr et al., 2013).
Following these two LSTIDs, smaller MSTIDs are seen propagating from the eruption source region for up to 6 hours. The first of these MSTIDs propagates at a speed of 337 ± 17 ms −1 with an initiation time of 4:32 UT and corresponds almost perfectly with the propagation of a surface pressure wave associated with the eruption (Figure 2a). This first MSTID and the MSTID after it are clearly visible out to ranges of at least 16,000 km and have dominant wavelengths of 300-350 km and periods of 10-18 min. Secondary MSTIDs with wavelengths ranging from 250 to 500 km and periods ranging from 10 to 25 min are also seen generated up to 6 hrs after the passing of the initial pressure disturbance. These secondary MSTIDs have speeds that range from 200 to 400 ms −1 , with most of the later MSTIDs having speeds toward the lower end of that range, consistent with the characteristics of ducted AGW waves reflected between the Earth and the lower thermosphere (Mayr et al., 2013).
The behavior of the ionospheric response to this event appears very similar to that of the 2011 Tohoku earthquake, but without the initial Rayleigh Wave (Tsugawa et al., 2011). Both events see initial LSTIDs with speeds near the ionospheric acoustic speed and subsequent MSTIDs with speeds near the surface acoustic speed.
One challenge with looking at the resulting disturbances in the manner shown in Figure 2a is that it assumes radial symmetry in the propagation of the TIDs. Many authors have noted preferential directional propagation of these structures associated with earthquakes (Inchin et al., 2020;Zettergren & Snively, 2019, for example). It is likely that highly directional wave fronts may not be apparent in representations like Figure 2a, which assume isotropy.
To get a clearer picture of the wave propagation without batching together all observed stations, in Figure 4 we generate analogous plots to Figure 3a but only for local regions with high GNSS station density. The boundaries of these regions are shown in Figure 1g. Using these plots, we can further examine the extent to which these waves propagate from the source.
In Figure 4, we see the abundance of waves of various sizes in the region directly adjacent to the eruption (New Zealand), and we can more clearly track the propagation of individual wave fronts. Interestingly, the LSTIDs do not appear to remain coherent in many other regions. Using the dotted propagation speed lines as reference, the first LSTID is seen clearly at New Zealand, Hawaii, and Northern Europe, with more inconclusive observation in Japan and Australia. Contrarily, the first LSTID is not apparent in South Africa and Eastern North America. Much of the same is seen for the second LSTID, but with a much clearer signature in Japan and a candidate signature in Eastern North America. Based on the geomagnetic location of the eruption and previous examinations of TIDs from similar sources, it is expected that waves will preferentially propagate along the magnetic field (Inchin et al., 2020;Zettergren & Snively, 2019). With the eruption being in the southern magnetic hemisphere, this would imply preferential northward propagation of these waves. Similarly, the declination of the magnetic field would imply preferential propagation to the North East and South West, aligned with the magnetic field. Examining Figure 4 again, one may note that the fast LSTID structures are barely discernible in Australia, suggesting substantial attenuation of the waves to the West. These modes are completely absent in South Africa, which is roughly along the same propagation trajectory from the source as Australia. This suggests a substantial disinclination of LSTID propagation to the West from the eruption. In contrast, the first and second LSTID signatures are the dominant structures seen at Hawaii to the North East, near the geomagnetic conjugate point of the eruption site. Interestingly, while the first LSTID front at New Zealand appeared to rapidly slow after ∼2,300 km range from the source, we see strong LSTID signatures at Hawaii ahead of the 700 ms −1 marker and very close to where we would expect a 900 ms −1 wave, consistent with the initial LSTID phase speed seen in New Zealand. It appears from this that, while the LSTID slowed substantially as it propagated over New Zealand no such reduction in speed occurred to the North East at Hawaii.
The appearance of the first and second LSTIDs in Northern Europe, which is at similar magnetic declination to Tonga, is suggestive of highly preferential propagation directly due North of the eruption site; however, a lack of stations in that direction over the Pacific, outside of the auroral oval, makes any conclusions regarding the appearance of the LSTIDs in Northern Europe challenging to assert. The strong presence of the LSTIDs at Hawaii is suggestive; however, more data is necessary. Further conclusions regarding the appearance of these LSTIDs in Northern Europe will thus have to await measurements from alternative sources, such as radio occultation. It should be further noted that geomagnetic activity was elevated during and prior to the day of the eruption, so it is possible for aurorally-driven equatorward propagating TIDs to also be present in these figures that could complicate analysis. We see this most clearly at New Zealand, Eastern North America, and South Africa in Figure 4, where we see equatorward propagating LSTIDs prior to the beginning of the eruption, and in Northern Europe, where we see an abundance of small and large scale structures that could, based on the auroral geometry, be aligned with the expected eruption wave propagation direction.
Subsequent MSTIDs originating from the eruption are seen at five of the seven test locations and arrive at times similar to the modeled atmospheric pressure disturbance. These waves are consistent with the speed of AGWs generated from the propagating surface perturbation, that is, Lamb wave. The MSTIDs propagating from this disturbance form the dominant global ionospheric wave response from this event, with clear signatures at nearly all sampled GNSS locations. In Northern Europe, MSTID wave signatures of eruption origin are far less obvious, likely due to interaction with waves of auroral origin and the sheer distance from the eruption location; furthermore, at Hawaii the MSTIDs are barely discernible. The lack of clear MSTID signatures at Hawaii may be a product of the lack of stations in that area; however, the observed amplitudes after the passage of the surface pressure anomaly are considerably weaker than the LSTID structures ahead of the surface pressure anomaly, which is only matched at New Zealand. Given that Hawaii is very near to the eruption geomagnetic conjugate point, it is possible that conjugate TID structures, generated by electric field perturbations associated with the TIDs at the eruption location, are interfering with TIDs propagating directly from the eruption (Zettergren & Snively, 2019). This, however, cannot explain why the LSTIDs do not suffer similar interference problems.

Conclusions
Perhaps the most striking feature of this event is the generation of an initial LSTID with a phase speed of >900 ms −1 . This wave is generated almost immediately after the event at nearly triple the surface acoustic speed and is close to the acoustic speed at the F-Region peak. This LSTID is shown to rapidly slow to propagation speeds near 700 ms −1 after propagating ∼2,300 km from the eruption location in the direction of New Zealand but did not exhibit such a reduction in speed at Hawaii, at even further ranges, and is perhaps a manifestation of a conjugate LSTID. Whilst the wave is a dominant structure in the vicinity of the eruption location, it appears to quickly dissipate depending on the direction of travel, where it is only barely discernible at Australia 4,000-6,000 km to the West of the eruption but forms the dominant wave response at Hawaii at similar ranges to the North East. Interestingly, the structure also seems to appear in Northern Europe between 12:00 and 13:00 UT, nearly directly due North from the eruption location. The initial LSTIDs generated from this event are likely acoustic waves generated from the initial eruption shock, while subsequent MSTIDs, propagating near the surface acoustic speed, correspond to AGWs generated by the subsequent atmospheric disturbance and persist to be generated up to 6 hrs after the eruption. While the shock-related LSTIDs demonstrate strong directionality and attenuate much faster than the MSTIDs, the MSTIDs exhibit no clear directional preference and form the dominant wave response seen globally except near the eruption conjugate point where these MSTIDs are much less dominant. Future work will need to combine modeling with additional measurements from radio occultation and of stratospheric/mesospheric winds to further elucidate the mechanisms at play.