2022 Tonga Volcanic Eruption Induced Global Propagation of Ionospheric Disturbances via Lamb Waves

The Tonga volcano eruption at 04:14:45 UT on 2022-01-15 released enormous amounts of energy into the atmosphere, triggering very significant geophysical variations not only in the immediate proximity of the epicenter but also globally across the whole atmosphere. This study provides a global picture of ionospheric disturbances over an extended period for at least 4 days. We find traveling ionospheric disturbances (TIDs) radially outbound and inbound along entire Great-Circle loci at primary speeds of ∼300–350 m/s (depending on the propagation direction) and 500–1,000 km horizontal wavelength for front shocks, going around the globe for three times, passing six times over the continental US in 100 h since the eruption. TIDs following the shock fronts developed for ∼8 h with 10–30 min predominant periods in near- and far- fields. TID global propagation is consistent with the effect of Lamb waves which travel at the speed of sound. Although these oscillations are often confined to the troposphere, Lamb wave energy is known to leak into the thermosphere through channels such as atmospheric resonance at acoustic and gravity wave frequencies, carrying substantial wave amplitudes at high altitudes. Prevailing Lamb waves have been reported in the literature as atmospheric responses to the gigantic Krakatoa eruption in 1883 and other geohazards. This study provides substantial first evidence of their long-duration imprints up in the global ionosphere. This study was enabled by ionospheric measurements from 5,000+ world-wide Global Navigation Satellite System (GNSS) ground receivers, demonstrating the broad implication of the ionosphere measurement as a sensitive detector for atmospheric waves and geophysical disturbances.


INTRODUCTION
The Tonga volcano eruption at 04:14:45 UT on 2022-01-15 was a huge geohazard event with farreaching effects, reportedly releasing 4-18 megatons (16-75 × 10 15 J) of thermal energy (Garvin, 2022) and causing a range of geophysical disturbances (Duncombe, 2022). Previous events and their effects in the charged upper atmosphere (e.g., Artru et al., 2005;Heki, 2006;Dautermann et al., 2009) are useful for comparison. The 1980 eruption of Mount Saint Helens was a VEI 5 (Volcanic Explosivity Index, see Newhall and Self, 1982) devastating disaster, comparable to the El Chichón eruption but less intense than the Pinatubo eruption at VEI 6. An estimated 24 megatons of energy release by this 1980 eruption (Kieffer, 1981) produced enormously impactful ionospheric disturbances at up to 9,000 km radius (Liu et al., 1982;Roberts et al., 1982). The reported ionospheric response to the Pinatubo eruption occurred at least 2,000-3,000 km distance across the Asian continent (Hao et al., 2006). Similar long distance effects occurred for the great Sumatra-Andaman earthquake (M 9.1 on Richter local magnitude scale) up to 5,000 km distance (Astafyeva and Afraimovich, 2006), and for the Tohuku earthquake (M 9.1) at up to 8,000-10,000 km distance in the US west coast (Crowley et al., 2016;Azeem et al., 2017).
Volcanic events can trigger severe disturbances that reach into the upper atmosphere above the epicenter, and in particular can produce periodic waves in both neutral and charged particles. A fundamental question for understanding the volcanic impact chain of response lies in characterization of the disturbance propagation mode in the upper atmosphere for given intensities of forcing and energy injection during the eruption. An eruption can excite both acoustic and infrasonic waves as compressional pressure waves, driving ionospheric plasma dynamics due to ion-neutral coupling. Tsunami waves are well known to be excited by the displacement of a large volume of water, and travel at a speed of~200 m/s for an ocean depth of~4,000 km (e.g., Astafyeva, 2019, and references therein). Ocean-atmosphere interaction via tsunami waves can induce atmospheric gravity waves which lead to ionospheric disturbances (e.g., Artru et al., 2005). In aggregate, these various volcano driven atmospheric wave modes are effective at causing ionospheric oscillations in the form of traveling ionospheric disturbances (TIDs) with periodicities spanning a few to 10s of min in the characteristic frequency domains of infrasonic, acoustic, tsunami, and gravity waves (e.g., Heki and Ping, 2005;Liu et al., 2006Liu et al., , 2010Hao et al., 2012;Zhao and Hao, 2015;Galvan et al., 2012;Zettergren and Snively, 2015;Chum et al., 2018;Astafyeva, 2019, and references therein). Ionospheric observations provide an effective and unique means of detecting these waves and some other oscillations, occurring in the entire atmosphere.
The extreme Tonga eruption provides an unprecedented scientific opportunity to gauge the global impact of this class of geohazard events on the whole atmosphere, and to improve our fundamental understanding of atmospheric wave characteristics during vertical and horizontal propagation. Themens et al. (2022) provided the first examination of a portion of the global extent of the ionospheric responses to the eruption, and reported some common TID modes as described earlier. Our study focuses on several important new features of eruption ionospheric effects. These include radially two-way (along full greatcircles) disturbance propagation in the global ionosphere for 4 days, and the fundamental roles of atmospheric Lamb waves that likely drove observed TIDs. These waves are recognized for the first time to cause a global impact over an extended period, well above their nominal dominant regime in the atmosphere.

METHOD AND DATA
We use GNSS total electron content (TEC) products from 5,000+ worldwide GNSS (GPS, GLONASS, and Beidou) receivers, generated (Rideout and Coster, 2006;Vierinen et al., 2016) and provided via the Madrigal distributed data system developed at the Massachusetts Institute of Technology's Haystack Observatory. In order to detect ionospheric responses associated with the Tonga eruption, we calculated differential TEC using an approach that effectively removes the background ionospheric "trend", as used in many previous TID studies (Zhang et al., , 2019aLyons et al., 2019;Sheng et al., 2020;Aa et al., 2021;England et al., 2021). Zhang et al. (2019a) provided more detailed discussions of this method. Differential TEC calculation of this nature is widely used for GNSS TEC based large and medium scale TID and ionospheric disturbance studies (Saito et al., 1998;Ding et al., 2007;Tsugawa et al., 2007;Komjathy et al., 2016;Zakharenkova et al., 2016;Azeem et al., 2017;Chou et al., 2018;Astafyeva, 2019).
The analysis uses individual receiver-satellite TEC data segments, subtracting a background TEC variation determined, in our technique, by a low-pass filtering procedure using a Savitzky-Golay low-pass filter (Savitzky and Golay, 1964). This residual is also called differential TEC (dTEC). We use a 30-min sliding window and a linear basis function for this particular study. To be completely free from impacts of the data edge associated with the use of a 30-min fixed length window, we removed data for the first and the last 15-min of each data segment. Finally, our analysis disregarded any data with satellite elevation < 15°. Final accuracy of this method ultimately derives from the accuracy of the GNSS phase measurement. Assuming that there is no loss of phase lock in the receiver, the error in differential TEC is less than 0.03 TEC units (Coster et al., 2012), as all satellite and receiver bias terms cancel out in a differential sense.

Global Extent of the Disturbances
The Tonga eruption provided an equivalent point source for observed atmospheric disturbances. We evaluated these disturbances based on the great-circle distance from the epicenter location (−20.5°N, −175.4°E) as identified by the US Geological Survey for the eruption induced magnitude M 5.8 earthquake origin (USGS, 2022). Figure 1 provides relevant geometry information and great-circle distance contours from the eruption location, as well as a great circle oriented at 26 and 206°azimuth from the epicenter. Superimposed is a background global map of GNSS TEC measurements at three post-eruption instances. Great circles assume 300 km height, characteristic of approximate ionospheric F region altitudes near the peak of the plasma population. The maximum great-circle distance is 21,000 km away from Tonga, located at (4.6°E, 20.4°N) near Sahara in North Africa. New Zealand was 2-4,000 km away, central US 12,000 km, South Africa 14,000 km, and Europe 18,000 km.
Upper atmospheric perturbations beyond 10,000 km have never been able to be examined before this eruption.
Both northern (|Azimuth| < 90°) and southern (|Azimuth| > 90°) great circles pass high latitudes. The great circle is presumed to be the shortest path along which disturbance energy and momentum in the neutral gas will flow radially from the epicenter. We note that, although global TEC is not evenly sampled by ground receivers with large gaps over the oceans, each observation is useful in distance-time analysis. Thus, in contrast to typical TEC studies, the distribution of disturbance propagation observations do not suffer severe gaps as demonstrated in the following distance-time figures. Themens et al. (2022) also presented these type of distance-time figures, and our analysis is similar except that we provide propagation estimates also based on azimuth bearing of great circles. The approach allows us to precisely locate propagation signatures and clearly identify inbound waves.
The distance-time variation of dTEC illustrated in Figure 2 indicate dramatic development of disturbance global propagation over a prolonged period. The southward propagation from Tonga to Africa sectors via the southern high-latitude region shows a defined envelope, as marked by fiducial arrows bounding enhanced disturbance (in dTEC) as a function of distance and UT. The width of the envelopes is~8 h in time with~350 m/s slopes. Results show that dTEC fluctuations reached the furthest distance at~20,000 km via the southern polar region. Northward propagation is predominately similar as indicated by envelope lines and their slopes, and also reached~20,000 km distance where it encountered the southern outbound propagation. Although dTEC signals became weak at several distances of 14,000 km and 16,000 km, corresponding to European sectors and midlatitudes, propagation signals reappeared beyond those distances perhaps due to wave modulation. In the following discussion, we examine detailed regional characteristics in near-field and far-field regions and provide further evidence of ionospheric perturbation arrivals.

Near-Field Ionospheric Disturbances
GNSS TEC measurements indicate immediate and vast near-field Tonga event atmospheric perturbations as demonstrated in FIGURE 1 | Geometry information of Tonga eruption impact distance (green lines) determined based on the great circle at 300 km height (white or yellow lines) that connects to the eruption region. Iso-distance lines up to 20,000 km are separated at 2,000 km interval. Great circles start at the Tonga epicenter for azimuth 26/206°. Background colors are differential TEC measured from ground-based receivers to GPS, GLONASS and Beidou navigation systems for the early stage of upper atmospheric responses at 0830 UT (A), 0620 UT (B), and 0920 UT (C). TID wave fronts are annotated by red and blue arrows in the three maps. Cyan lines are isogeomagnetic latitudes at the 15°interval.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org March 2022 | Volume 9 | Article 871275 Themens et al. (2022) and Figure 3. The earliest response was a clear positive dTEC occurring within 200 km of the epicenter almost instantaneously following the eruption at~04:15UT. This response, with~1 km/s radial propagation for the first 20 min, was an indication of supersonic infrasonic waves typically seen (as Rayleigh waves) during earthquake events (Astafyeva, 2019, and references therein). Immediately following, two enormous shocks occurred with dTEC magnitudes up to 3 TECu (1TECu = 10 16 electrons/m 2 ). Radial propagation initially occurred at 700 m/s speed, gradually slowing down to~450 m/s, and reached~5,000 km distance. The initial waves were clearly identifiable over the northern New Zealand area (~1,500 km away from the epicenter) as early as~0500 UT (e.g., Figure 3) and, specifically, at~06:20 UT with 2-D fronts ( Figure 1B). Subsequent waves were characterized by smaller amplitudes (0.1-0.2 TECu) at lower and relatively stable speeds of 360 m/s. These amplitudes were well within typical amplitudes for medium to large scale TIDs. These fluctuations had 10-30 min quasi-periodicity for at least 8 h (Figures 2,   3A-3C). The 2-D wave fronts exhibited horizontal wavelengths initially 500-1,000 km ( Figure 1B) and later 300 km ( Figure 1C).

Far-Field Ionospheric Responses
Beyond the near-field, outbound ionospheric responses propagated into different regions of the world. Within 195-315°azimuth bearing of the great circles, the disturbance propagation signals were evident between 13-18,000 km ( Figure 4A), particularly over South Africa, with large amplitudes and consistent~350 m/s propagation speeds. These disturbances lasted for at least 5 h with up tõ 10-30 min periods, arriving via southern great circles ( Figure 1A). These results were derived as dTEC averages in time and distance with 1 min bins size in time and 10 km bin size in distance. 2-D wavefronts during the shock arrival showed well organized disturbances with 500-1,000 km wavelengths ( Figure 4C), similar to Figure 1B in the nearfield.
FIGURE 2 | Distance-UT variation of dTEC for disturbance propagation southward (negative distance) and northward (positive distance) along the great circle paths at 300 km altitude on 15 January. White arrows provide envelope lines encompassing the ionospheric disturbances. The slopes of these lines are~350 m/s. Dashed lines with larger slopes (~700 m/s) follow the initial ionospheric shocks which terminated after 5,000-6,000 km. The two short red arrows mark the radial propagation in Japanese/US and European sectors.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org March 2022 | Volume 9 | Article 871275 The continental US (CONUS) has dense receiver networks, therefore only a narrow range of azimuths (55-58°) are taken into account, to minimize decoherence of the wave signature due to regional deviations in the group velocity of the wave fronts. Such deviations can be caused by e.g., prevailing wind velocity, atmospheric pressure, propagation direction (Taylor, 1932), and the ionospheric conditions. Figure 4B shows the first sign of disturbance arrival in the west coast (8-9,000 km distance) at 11-12:00 UT, and the earliest front departed off the east coast at 13,000 km by~16:00 UT. Throughout, propagation speed remained at~350 m/s. Some samples of 2-D wavefronts are shown in Figure 4E. dTEC enhancements were aligned with iso-distance lines, and separated zonally by~500 km spacing (wavelengths). Simultaneously, background TIDs were present, likely associated with geomagnetic disturbances during the recovery from a geospace storm with minimum Dst −94 nT at 23:00 UT on 2022-01-14 (according to WDC Kyoto, 2022b). The substorm FIGURE 3 | Near-field observations of initial and subsequent GNSS TEC fluctuations: the distance-time (regardless direction) variation within 5,000 km 6 h following the eruption (A); regional GNSS TEC fluctuations in New Zealand showing the evolution of fluctuation periodicities in space and time (B); near-field TIDs, the same as (A) but over 48 h (C) with red arrows marking the outbound~350 m/s wave propagation, and black arrows marking the potential returning waves at~350 m/s into Tonga after 15:00 UT on the following day 16 January. See Figure 4 for further indications of returning waves.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org March 2022 | Volume 9 | Article 871275 5 activity measured by the AE index (WDC Kyoto, 2022a) indicated an onset after 1100 UT and multiple hourly enhancements reaching~500 nT throughout the rest of the UT day. Thus the ionospheric short-term disturbance was characterized by hourly enhancements of TID amplitudes almost simultaneously in latitude and longitude over CONUS, well correlated to AE activity. Post-sunrise TIDs were evident in the US east coast, propagating eastward at~300 m/s with 20-30 min periods, similar to what were reported in Zhang et al. (2021). Background TIDs include also gravity waves in certain region potentially linked to the strong stratospheric polar vortex (Sato and Yoshiki, 2008;Becker and Vadas, 2018;Bossert et al., 2020). These background TIDs occurred on the distancetime plot as simultaneous enhancement bands, slant bands, and possibly other features. The arriving eruption-induced fluctuations segmented these pre-existing TIDs fronts into smaller structures elongated along iso-distance lines, and had sharply different characteristics that can be distinguished in the distance-time plot ( Figure 4B).

Wave Propagation Return to Tonga
Ionospheric fluctuations continued to propagate through the eruption antipode in North Africa toward Tonga on the next day. These returning TIDs were most evident over South Africa ( Figure 4D) where the disturbance phase was clearly toward shorter distance (light-blue arrows), starting at 03:00 UT at 15,000 km distance. The speed was~350 m/s. Figure 4E shows an example of wavefronts associated with the returning TIDs at 03:30 UT on 2022-01-16, close to local sunrise time. The timing of the returning TIDs is approximately consistent with propagation from the most distant point to South Africa along the great-circle Tonga-CONUS-Europe-South Africa-Tonga path ( Figure 1A). Returning TIDs occurred also across the CONUS over a prolonged period for~8 h initially at~0600 UT at 13,000 km (the US east coast) (see light-blue arrows in Figure 4B), following a longer path of Tonga-southern high latitudes-South Africa-Europe-CONUS.
In the near field, clear indications occurred of the wave returning to New Zealand ( Figure 3C) at 2,500-3,000 km distance from Tonga by~13:00 UT on 16 January (or to Tonga by~15:00 UT on 16 January), after traveling nearly 1.5 days along the complete great-circle. This timing is roughly consistent with a propagation at~350 m/s speed around the full great-circle.

Discussion
During the TID global propagation, the horizontal phase speed varied between 300-350 m/s depending on propagation direction. For example, Figure 2 marks a southward propagation (red arrow) at 300 m/s. However, these speed estimations are generally consistent with infrasonic detection of pressure wave arrival at individual stations around the world, e.g., in northern Europe (Norstar Website, 2022) using the network established to monitor compliance with the Comprehensive Nuclear Test Ban Treaty and over CONUS using the pressure altimeter observations (Iowa Environmental Mesonet, 2022). Lamb waves travel at the sound speed, typically 300-350 m/s in the troposphere, and can exist at any period. Although their energy is confined to the troposphere, their amplitudes increase exponentially with height due to decreasing density. Their wave energy can leak into the upper atmosphere when Lamb waves with~300 m/s horizontal phase speed are resonant with the atmosphere, as can be the case with acoustic and gravity waves (Bretherton, 1969;Lindzen and Blake, 1972;Nishida et al., 2014). Lamb waves with~319 m/s phase speed were previously identified as an atmospheric wave response to the Krakatoa eruption (Symons, 1888;Taylor, 1932). Similar Lamb waves were also detected by very sensitive microbarographs during the St. Helens eruption (Mikumo and Bolt, 1985) and the Sumatra-Andaman earthquake (Mikumo et al., 2008). A significant portion of our GNSS observations occurred inland, where direct tsunami wave contribution to these TID results may be ignored. Furthermore, an additional argument against the presumption of tsunami wave presence with 10-30 min periodicity and 300-350 m/s travel speed as noted earlier: Tsunami waves were reported to occur in the US west coast at 15:30-16:00 UT (NOAA DART and NOAA/NOS/ CO-OPS Data, 2022) consistent with an anticipated 210-220 m/s propagation speed across the Pacific ocean. This is clearly different from observed ionospheric wave propagation which arrived at least~4 h earlier than the tsunami waves in the US west coast. Nevertheless, continued community study is recommended to further clarify the roles of tsunami and gravity wave interactions and the factors that are potentially responsible for their different propagation speeds Kherani et al., 2016;Bagiya et al., 2017).
The earliest wavefronts could be seen traveling around the globe three times, passing six times over the CONUS over 100 h since the eruption ( Figures 4B,F); the pass 6 over CONUS occurred at~05:00 UT on 2022-01-19 at 13,000 km (but is not shown here due to the space limit). The travel time around the globe in the direction of Earth's rotation (eastward) was 34.8 ± 0.7 h and in the opposite direction 36.0 ± 0.7 h. The error standard deviations are roughly one order of magnitude estimates based on visual inspection. The measured propagation speed and the number of observed passes of the atmospheric wave are comparable to those reported with the Krakatoa eruption (Symons, 1888). Time periods before and after the eruption were processed with the same analysis, and no similar traces corresponding to ones shown in Figures 4B,F were observed before the eruption.

SUMMARY
The 2021 Tonga volcano eruption caused enormous and truly global perturbations in the ionosphere and thermosphere over an extended period. Ionospheric disturbances were observed traveling three times around the globe. They returned back to Tonga every 1.5 days. The ionospheric responses were characterized by an immediate supersonic (~1 km/s) plasma density impulse, and two shock waves with substantial amplitudes (3 TECu) and 600-700 m/s horizontal speeds in ≤5,000 km near-field regions including New Zealand to the south, and Hawaii to the north. Subsequently, persistent (lasting 8 h) slower-propagating TIDs developed with 10-30 min periods, most significantly in the near-field (≤ 5,000 km radius). Far-field wave effects include also 10-30 min periodical oscillations lasting for a few hours, but more significantly, TID shock fronts were clearly organized based on the distance from Tonga in the continental US (~12,000 km distances) and midlatitude west Europe (~18,000 km distances), with the two destination regions connected via the northern highlatitude region along the great circle with an origin in Tonga. A similar path of global wave propagation occurred in the southern hemisphere along a New Zealand-southern polar region-South Africa route. Other far-field regions including South America and Frontiers in Astronomy and Space Sciences | www.frontiersin.org March 2022 | Volume 9 | Article 871275 Japan saw the arrivals of these TID shock fronts. These far-field globally propagating waves had 500-1,000 km horizontal wavelengths. These disturbances resulted from eruption-induced perturbations at frequencies of acoustic waves (including Lamb waves) and gravity waves. Eruption-associated tsunami waves were slower than the main component of ionospheric waves that propagated globally and are therefore unlikely responsible for the TID global propagation. The presumption of Lamb wave global propagation at 300-350 m/s is consistent with our main observational results. These waves provide one of the main carriers for eruption energy leaking into the upper atmosphere because of atmospheric resonance to forcing provided by these waves at~300 m/s phase speed, equivalent to the speed of sound in the troposphere. These ionospheric propagation results are also consistent with data from infrasonic global detections and other pressure wave detections. Our multi-sensor investigation, based on 5,000+ world-wide GNSS receivers, reveals the unprecedented depth, severity, and extent of disturbances in the whole atmosphere in vertical and horizontal dimensions that occur during an extremely devastating geohazard impact. This is yet another demonstration of the ionosphere acting as a sensitive detector for atmospheric waves.

DATA AVAILABILITY STATEMENT
GNSS TEC data products and access through the Madrigal distributed data system [http://openmadrigal.org] are provided to the community by the Massachusetts Institute of Technology. Further inquires regarding the datasets generated for this study can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
S-RZ was responsible for differential TEC derivation, science analysis of the observational results, preparing the initial manuscript and organizing team efforts. JV contributed initial figures for comparisons, science analysis, discussion, and organizing team research. EA was responsible for Beidou GNSS data processing. LG conducted analysis of TIDs that could potentially be linked to meteorological disturbances. PE contributed substantially to the manuscript development. AC was responsible for GNSS data management. WR was responsible for software development of GNSS data processing and daily GNSS data processing. AS contributed examining satellite data. All members contributed to science discussion and manuscript development.

FUNDING
GNSS TEC data products and access through the Madrigal distributed data system are provided to the community by the Massachusetts Institute of Technology under support from US National Science Foundation grant AGS-1952737. MIT staff members were partially supported by NASA grants 80NSSC21K1775, 80NSSC21K1310, 80NSSC22K0171, and 80NSSC19K0834, AFOSR MURI grant FA9559-16-1-0364, NSF grant AGS-2033787 and ONR grant N00014-17-1-2186.