© T.G.W. Verhulst et al., Published by EDP Sciences 2022

1 Introduction

After a series of smaller disturbances starting at the end of 2021, on January 15, 2022, at 04:15 UT, the Hunga volcano in Tonga (20.54°S, 175.38°E) (Global Volcanism Program, 2013; Cronin et al., 2017) produced a powerful eruption (Global Volcanism Program, 2022). In terms of the total energy released, this was the largest volcanic eruption since the 1991 eruption of Mount Pinatubo. Large eruptions are known to produce wave-like disturbances up to ionospheric altitudes. Such travelling ionospheric disturbances (TIDs) have been observed, for instance, after the Pinatubo eruption of 1991 (Igarashi et al., 1994), as well as the 2004 eruption of the Asama volcano in Japan (Heki, 2006), the 2003 eruption of Soufrière Hills on Montserrat (Dautermann et al., 2009), and even earlier already in the context of the 1980 eruption of Mount St. Helens (Roberts et al., 1982). Nevertheless, due to the infrequent occurrence of large volcanic eruptions, fewer studies have been done on their ionospheric effects than, for instance, on the effects of earthquakes and tsunamis (Astafyeva, 2019) or strong tropospheric circulation (Šindelářová et al., 2009; Nishioka et al., 2013).

TID signatures were also observed after other explosive events (Huang et al., 2019). In particular nuclear tests have been known to cause ionospheric disturbances (e.g., Breitling & Kupferman, 1967; Hines, 1967; Kanellakos, 1967; Albee & Kanellakos, 1968; Park et al., 2013; Zhang & Tang, 2015). Other large explosive events, whether anthropogenic, such as conventional explosions (Barry et al., 1966; Fitzgerald & Carlos, 1997; Drobzheva & Krasnov, 2006) and rocket launches (Chou et al., 2018), or natural, such as superbolide meteor impacts (Pradipta et al., 2015; Luo et al., 2020), can also produce medium-scale TIDs (MSTIDs). Some historic events, such as the industrial accident in Flixborough in 1974 (Jones & Spracklen, 1974; Krasnov et al., 2003) and the conventional bombing campaign during World War II (Scott & Major, 2018), have recently been reanalysed from this perspective. However, it has been shown both from observations and from theoretical considerations that the MSTIDs generated through different mechanisms can be very distinct from each other (Kirchengast, 1997; Huang et al., 2019). Clear differences have been observed when comparing disturbances generated from highly localised explosive events such as nuclear explosions and volcanic eruptions (Roberts et al., 1982; Huang et al., 2019). Ionospheric signatures of less localised phenomena such as tsunamis and seismic events are even more different (Astafyeva, 2019; Huang et al., 2019). Because of the uniqueness and rarity of such violent eruptions, it is important to carefully analyse all observational data related to this event.

Two different physical mechanisms can cause TIDs to appear after a volcanic eruption: disturbances can be produced directly in the ionosphere at the location of the eruption and travel radially outwards at ionospheric altitudes, or ionospheric disturbances can result from the propagation of various types of waves through the lower atmosphere. The waves propagating through the lower atmosphere are further distinguished into acoustic waves, gravity waves and Lamb waves, depending on their frequency compared to the acoustic cutoff frequency (Yeh & Liu, 1974; Haaser et al., 2017). A volcanic eruption can produce waves of all these types, propagating outward from the eruption with different velocities and to different distances before dissipating. For this event, a Lamb wave was detected as the leading wavefront, followed by various disturbances of different natures (Burt, 2022; Kubota et al., 2022; Kulichkov et al., 2022; Saito, 2022). In this study, we look for disturbances in the ionosphere over Europe, close to the antipode of the eruption. MSTIDs propagating through the ionosphere from the site of the eruption are not expected to proceed to such distances (Haaser et al., 2017; Astafyeva, 2019). Themens et al. (2022) indeed find that large-scale TIDs produced directly in the thermosphere above the eruption dissipate after a few thousand kilometres. However, medium-scale TIDs travelling with the Lamb wavefront seem to propagate further (Themens et al., 2022; Zhang et al., 2022), and clear signatures of waves travelling around the globe in the lower atmosphere have been detected (Burt, 2022; Matoza et al., 2022; Wright et al., 2022). Waves might thus be expected to appear in the ionosphere at the same time as the ground-level pressure disturbances, for the Lamb waves, or with some delay, for disturbances propagating upwards from the ionosphere.

Various relatively dense networks of observatories are available in Europe to monitor the ionosphere using different instruments. This provides us with an opportunity to combine results from multiple data sources. We analyse data from both vertical and oblique ionogram traces, GNSS-derived total electron content (TEC), as well as in situ measurements from the Swarm C satellite when it passed over the region. In Section 2, we first describe the geomagnetic background conditions during the period of interest, followed by a description of the various observatories and data types used in this study. In Section 3, we systematically describe all the different observations from the various instruments, followed by Section 4 with a discussion of how the various observations relate to each other. Finally, in Section 5, we summarise the conclusions of our analyses. In this last section, we also give some ideas and recommendation for improving the observation of MSTIDs such as those seen here in the future.

2 Data and methods

2.1 Geomagnetic background conditions

The period before and after the eruption on January 15 saw some significant geomagnetic disturbances. It can be seen from the D_st shown in Figure 1 that this event took place during the recovery phase of a geomagnetic storm. The D_st reached a minimum of −91 nT on January 14 at 23 UT. The bottom panel of Figure 1 shows the local hourly K-index derived for the geomagnetic observatory in Dourbes, co-located with the ionosonde DB049. It is evident from this figure that some moderate geomagnetic disturbances were detected in Europe after 19 UT on January 15, until the early hours of January 16. This is precisely the period during which the various waves related to the eruption are expected to arrive in Europe. During periods of geomagnetic activity, large-scale TIDs are often produced in the auroral regions, travelling towards lower latitudes. A major challenge is, therefore, to distinguish TIDs produced by gravity waves in the lower atmosphere from those appearing at the same time from auroral sources.

Fig. 1 Top: real-time Dst provided by the Kyoto WDC for Geomagnetism; Bottom: local magnetic K-index from the Dourbes observatory co-located with the ionosonde DB049 at 50.1°N, 4.6°E. The vertical red line indicates the moment of the explosive eruption, while the purple band shows the interval during which volcano-induced MSTIDs are detected over Europe.

2.2 Vertical and oblique ionogram soundings

Table 1 lists the ionosondes from which data is used. Their locations are also shown in Figure 2. Together, these twelve ionosondes provide a relatively dense set of observations over Europe. It should be noted that the sounding cadences are not the same at all observatories. The highest time resolution is available at the five ionosondes operating at a five-minute sounding interval. Other observatories use intervals of up to 15 min. This is particularly important, considering the disturbances caused by a volcanic eruption are expected to fall in the MSTID range and might not be readily evident from observations with coarser time resolution.

Fig. 2 Location of the 12 ionosondes used in this investigation, as well as the 4 oblique sounding links.

The final two columns in Table 1 list the shortest distance over the Earth’s surface between the eruption and each observatory and the azimuth from the ionosonde to the eruption. It can be seen from the wide range of azimuths that disturbances arrive at various places in Europe from widely different directions. This is due to the studied region being close to the eruption antipode, located in southern Algeria. As a consequence, there can be some differences in the arrival times as different great-circle paths will pass through different atmospheric conditions, resulting in different average sound velocities.

For the included ionosondes, ionograms were obtained directly from the instrument operators or from the Global Ionosphere Radio Observatory repository (Giro) (Reinisch & Galkin, 2011). From 18:00 UT on January 15 to 06:00 UT on January 16, the period during which effects from the eruption are expected to be visible, all ionograms have been manually scaled. For each observatory at least the critical frequency f_oF₂, the height hmF₂ and the maximum usable frequency MUF(3000)F₂ (further simply called MUF(3000)) of the F₂ layer were scaled. The MUF(3000) is determined by fitting an empirical curve to the trace in the ionograms (Piggott & Rawer, 1978; Paul, 1984) and will be affected by variations both in f_oF₂ and variations in hmF₂. From previous work, for instance, Altadill et al. (2020), it is known that the MUF(3000) often exhibits a clear signature when a TID arrives. Furthermore, scaling the MUF(3000) for oblique traces is relatively straightforward. Because of this, in the current work, we primarily use the MUF(3000) from the ionosonde characteristics.

In addition to the parameters scaled from vertical ionograms, we include data from four oblique sounding paths. Oblique ionograms are produced by synchronising identical vertical ionogram soundings at different observatories. Oblique traces are then visible together with the vertical ones, allowing for the scaling of the oblique MUF at a distance between the ionosondes (see Verhulst et al., 2017, for more information). Details of the oblique-sounding paths used here are given in Table 2, and the paths are also shown in Figure 2. For each oblique sounding path, only the oblique traces at one of the ionosondes was selected for analysis, based on the relative signal quality and signal-to-noise ratio achieved in both directions.

The scaling of ionogram traces, both from vertical and oblique echoes, was, in some cases, hampered by the geophysical conditions. As discussed in Section 2.1, there were some geomagnetic disturbances throughout the period under consideration. This caused some instances of spread-F traces at the observatories at higher latitudes, particularly after January 16, 00:00 UT (corresponding to a slight increase in geomagnetic activity after midnight, see Fig. 1). In addition, the median night-time electron density at this time of year can become very low, making it difficult to confidently scale ionospheric parameters. This, again, mainly affects the most northern ionosondes. Parameters were only scaled when this could be done with high confidence, leading to some gaps in the time series.

For the vertical ionograms, we use the MUF at 3000 km for all observatories. In order to obtain the same parameter from the oblique traces, the virtual heights of the reflections should be determined. However, the distance at which the MUF(3000) is calculated has no bearing on the timing of possible peaks resulting from TIDs. Such conversions can even be a source of additional uncertainty, stemming from the uncertainty in determining the virtual heights. We consider here for each oblique sounding path the MUF at the length of the path, avoiding the need for conversion to a different distance.

Detrended iso-ionic data dh(f) provide information about local disturbances in the height of HF frequency reflections due to TIDs. The basic characteristics that can be calculated with this method are the oscillations’ periods and the amplitudes. Depending on the cadence of the vertical incidence ionograms, the method, in general, can be employed to detect both medium- and large-scale TIDs. The disturbances studied in this contribution are driven by lower atmosphere acoustic waves, and therefore the highest possible cadence is required in order to detect TIDs of medium scale. In the results presented in Section 3.2, the detrending is applied for some of the sounders providing data at a five-minute cadence, using a running window of 1 h. The iso-ionic data are extracted from manually scaled SAO files.

The MSTIDs expected to be observed in Europe due to a volcanic eruption on the other side of the world are caused by pressure waves travelling in the troposphere, which produce disturbances moving up to the ionosphere. Therefore, it is interesting to compare ionospheric observations with pressure measurements from barometers at the same locations as the ionosondes. In this work, we consider ground level pressure data from barometers co-located with the DB049 and EB040 ionosondes. In both cases, pressure measurements are made with a time resolution of 1 min, allowing the reliable identification of passing disturbances with periods of tens of minutes.

2.3 Ionosonde drift measurements

The DPS-4D ionosondes can be configured to produce Digisonde Drift Measurements (DDM) soundings, using a longer transmission on a small band of frequencies, in addition to the traditional ionogram soundings. From the angle of arrival and range of the registered echoes, the tilt of the reflecting layer can be deduced, while the Doppler shifts can be used to calculate the bulk plasma drift velocities.

It has been known for a long time that plasma drift observations can be used to detect TIDs (MacDougall, 1966). Paznukhov et al. (2020) showed that also variations in the ionospheric tilts can be produced by TIDs, and can be monitored for their detection and characterisation. However, obtaining drift and tilt observations from ionosonde soundings presents some difficulties. In particular, the calculation of tilts from DDM soundings is not straightforward. We looked at tilt data from a few observatories, but we were unable to draw reliable conclusions from it. Therefore, we do not include the tilt observations in this work.

A major challenge for these measurements lies in selecting a suitable transmission frequency. On the one hand, keeping the examined ionospheric height range small, but on the other hand, producing a sufficiently high number of echoes to obtain reliable drift data. The frequencies used during this event were the ones employed for routine night-time operation at each observatory. Because of the ongoing geomagnetic disturbances described in Section 2.1, these were not always optimal. We only included drift observations from a few selected observatories where a decent amount of good-quality data was available.

The Digisonde Drift Analysis software (Kozlov et al., 2008) extracts automatically the three-dimensional plasma drift velocity components from DDM records. However, as is the case with the automatically scaled ionograms, manual verification of all data is required (Kouba et al., 2008). For example, in Figure 3a, good quality skymap observation is shown on the left and a problematic one on the right. The former contains 399 echoes, all from ranges between 487 and 515 km. Although some echoes are spread over the sky, there is a clear clustering around one location, and the Doppler shifts are consistent. The latter skymap only comprises 70 echoes in total, with ranges spread between 267 and 415 km. Echoes with different Doppler shifts are present at different zenith and azimuth angles. This is indicative of a skymap containing echoes from several layers or multiple reflections. Data needs to be manually filtered and verified before calculating the plasma drift velocity to obtain a sensible result. Data is filtered to contain only echoes from between 190 and 400 km, which are the single reflections from the F-layer.

Fig. 3 Two skymap observations from January 15 produced by the EB040 Digisonde. The left image shows the result of the sounding at 23:04 UT and the right one at 23:14 UT. Each mark represents the direction of a recorded echo, with the colour indicating the Doppler shift. The arrows drawn on the image indicate the automatically calculated bulk plasma drift.

2.4 TEC and in situ data

GNSS data from different networks covering the European sector – shown in Figure 4 – have been used to investigate the impact of the Hunga eruption on the ionosphere in terms of TEC. To highlight the effect of the wave-like structure induced by the explosion, TEC data have been detrended using the Varion algorithm. Varion is an open-source, Python-based software (available at https://github.com/giorgiosavastano/VARION), described in Savastano et al. (2017) and Ravanelli et al. (2021). This algorithm is based on the computation of the integral over a certain interval of the time differences of geometry-free combinations of carrier-phase measurements from a stand-alone GPS station that reads:

(1)where dTEC is the detrended TEC, t_i represents the initial time of the considered period, is the geometry-free combination of the carrier phase measurements calculated considering the receiver R and the satellite S, and f₁ and f₂ are the GPS L1 and L2 signal frequencies, respectively (Ciraolo et al., 2006). Varion uses of the standard orbit and clock products, based on a thin layer approximation of the ionosphere located at 350 km altitude. As the time difference of the carrier phase measurements, the effect of the inter-frequency biases on TEC evaluation can be ignored since they can be considered constant along every single arc, if no cycle slips occur and no verticalisation is applied.

Fig. 4 Location of the GNSS receivers used for detrended TEC analysis.

The long-term trend is removed from the dTEC time series by computing the residuals with respect to a tenth-order polynomial fit. To remove the presence of wiggles at the arc boundaries, an elevation mask of 20° is then applied. Several algorithms are available in the literature for the detection of wave-like ionospheric structures with GNSSs measurements (see, e.g., Saito et al., 1998; Komjathy et al., 2005; Hernández-Pajares et al., 2006; Galvan et al., 2011; Belehaki et al., 2020; Maletckii & Astafyeva, 2021). Differences lie in the adoption of the slant or vertical TEC, the different use of calibrated TEC values, and the different employment of bandpass filtering methods. Varion has been validated against the TID detection method developed by Jet Propulsion Laboratory (Komjathy et al., 2005). We further compared Varion results against some of those techniques, finding no meaningful differences in the capability of detecting the presence of MSTIDs and of estimating its period. Therefore, we only discuss the Varion results in this work.

In situ measurements from the Swarm Charlie (C) satellite (Friis-Christensen et al., 2008) have also been considered to investigate the signature of the ionospheric variations related to the explosion. In 2022, Swarm A and C flew closely at around 435 km altitude (speed v_A = v_C = 7.2 km/s) and Swarm B at around 507 km (v_B = 7.0 km/s). For the purposes of our analysis, we consider only Swarm C because Swarm A data are not available for the period under consideration. However, significant differences in the detection of the perturbance between plasma density measured by the two spacecraft are not expected because they are closely separated (1.1° longitude at the equator). Swarm B, is located at higher altitudes and covering different local and universal time sectors with respect to Swarm C, recorded less evident signatures of the effect of the explosion. For our purposes, we consider the electron density (Ne) provided by the Langmuir Probe on board the spacecraft and available at a 1 Hz rate – downsampled from original 2 Hz observations – in the global Ionospheric Plasma IRregularities product based on Swarm measurements (Jin et al., 2022).

To highlight the possible signatures of the wave-like structure induced by the explosion in the Swarm C Ne measurements, we select tracks passing in the European sector and in time intervals that are compatible with the theoretical time of arrival. Figure 5 shows the selected tracks and the corresponding hours after the explosion. The longitudinal sectors of the two Swarm C passages are 28.9°E and 5.5°E, respectively.

Fig. 5 Location of the selected Swarm C tracks. Colour code reports the hours after the eruption.

Because the speed of the swarm satellite (about 7.2 km/s) is much larger than the speed of the ionospheric waves (about 300 m/s), the ionospheric waves are considered fixed with respect to the satellite (Kil & Paxton, 2017; Urbar et al., 2022). To identify wave-like signatures in Swarm Ne data, we make use of the novel Fast Iterative Filtering technique (Cicone, 2020; Cicone & Zhou, 2021). This procedure can be used to decompose a non-stationary, non-linear signal s(t) into simple oscillatory components, called intrinsic mode components or intrinsic mode functions (IMFs), according to the following formula:

(2)

In equation (2), N_IMF is the total number of IMFs, each having its own (quasi-stationary) frequency ν, from which we can derive the associated wavelength λ. res is the residual, which is assumed to be the background trend, as it contains no further oscillatory components. In the Swarm case, s(t) = Ne(t) = Ne(lat, lon). The advantage of using Fast Iterative Filtering comes from the fact that it is based on a complete theory, with its convergence and stability having been mathematically proven, and that it allows for a very fast calculation, being in the order of a hundred times faster than methods based on the use of Empirical Mode Decomposition (Huang et al., 1998).

Moreover, to further highlight the portion of the Ne signal driven by the passage through the ionospheric wave-like structure (hereafter dNe), we sum up the IMFs to provide a Ne representation as follows:

(3)in which Ne^Noise contains all the high frequency, small-scale variations and Ne^Trend is the non-oscillatory trend, as per equation (2).

Fast Iterative Filtering and methods based on the implementation of the iterative filtering concept have already been proven very effective for providing a fine time-frequency representation of signals in ionospheric physics applications (see, e.g., Materassi et al., 2019; Piersanti et al., 2018; Ghobadi et al., 2020; Spogli et al., 2019, 2021; Urbar et al., 2022).

3 Observations

3.1 Atmospheric pressure wave and TID arrival times

Figure 6 shows the MUF(3000) for ionosondes DB049 and EB040, together with ground level air pressure measurements from co-located barometers and median values of MUF(3000) taken over the whole of January 2022. At both locations, the time resolution for the ionosonde data is 5 min and for the barometer data 1 min. The clearest signature in the pressure data is the first wave arriving in Dourbes, which corresponds to the shortest of the four travel paths.

Fig. 6 Air pressure (black lines), MUF(3000) data (red lines), and monthly median MUF(3000) for January (orange line) for Dourbes (top) and Ebre (bottom) for the period from 12 UT on January 15 to 12 UT January 16. Note that there are some gaps in the MUF(3000) data, especially at Dourbes, due to unscalable ionograms.

The times of the highest peaks – positive or negative – in the pressure are at 19:29 UT and 01:40 UT for Dourbes and 20:40 UT, and 00:42 UT for Ebre. The average velocities for the pressure waves calculated from the arrival time and the distances given in Table 1 are 303 m/s and 304 m/s for the two waves arriving in Dourbes and 300 m/s and 304 m/s for the waves arriving in Ebre.

The background values for the MUF(3000) are clearly influenced by the geomagnetic disturbances described in Section 2.1. During the daytime, there is a clear enhancement, while the background value during the night is somewhat lower than the monthly median. Superimposed on this depleted background, the peaks associated with the MSTIDS can be seen around 21:30 UT in Dourbes and around 22:20 UT and 01:30 UT in Ebre. The increase in MUF(3000) in Dourbes before 01:00 UT is followed by a data gap due to spread-F conditions associated with the geomagnetic disturbances at this time and is likely also related to that. In addition, some data are missing for the period during which the second MSTID peak should be observed due to unscalable ionograms.

Note that the MUF(3000) only shows a single peak at the time of arrival of each TID. Especially at Dourbes, multiple periods for the pressure wave can be seen in Figure 6. The total time for which the pressure is showing a disturbance here is close to one and a half hours. For Ebre, because it is closer to the antipode of the eruption, the first pressure wave has not yet entirely subsided by the time the second wave arrives from the other direction. This results in a single, continuously disturbed period with two maxima in disturbances but no completely quiet period in between, such as can be seen in Dourbes after 22:00 UT.

The time between the highest maximum and the subsequent minimum in the pressure observed at Dourbes is 18 min. The acoustic cutoff depends on the atmospheric conditions but is normally around ω_a = 3.3 MHz in the lower atmosphere (Astafyeva, 2019), which corresponds to a period of about 5 min. Therefore, the waves detected here are well within the gravity wave regime. Gravity waves are expected to take 40–60 min to travel from ground level up to ionospheric heights (e.g., Astafyeva, 2019), which is in agreement with the data shown here.

3.2 Iso-ionic contours

Detrended iso-ionic true height data are shown in Figure 7. In the right side panels, the plots for the period during which MSTIDS are expected, 2022-01-15 18:00 UT to 2022-01-16 06:00 UT, are presented. On the left, the plots that correspond to the same time interval but occurred three days before, when the geomagnetic and auroral activity were at quiet level, are shown for comparison. The iso-ionic contours are created using data from the Athens (AT138) and Ebre (EB040) ionosondes. Data from these specific Digisonde stations are used here because of the completeness of the time series of observations and the five-minute cadence. Furthermore, the ambient ionospheric density, which is higher in these two stations compared to the higher latitude observatories, creates favourable conditions for the propagation of TIDs with higher amplitudes (Hunsucker, 1982; Reinisch et al., 2018).

Fig. 7 Detrended iso-ionic true height dh(f) plots for a geomagnetically quiet interval – 2022-01-12 18:00 UT to 2022-01-13 06:00 UT – are presented in the left column for the Athens (AT138, top) and Ebre (EB040, bottom) Digisonde stations. The corresponding plots for the period of interest here, from 2022-01-15 18:00 UT to 2022-01-16 06:00 UT, are shown in the right column.

While some weak fluctuations of random nature are observed during the quiet period, systematic oscillations in heights occur during the disturbed period for all frequencies ranging from 2.0 to 4.0 MHz. These oscillations have an average amplitude of approximately 50 km for AT138 and 75 km for EB040, and average periods of 45 and 40 min, respectively. The small differences in amplitudes and periods between the two locations can be understood considering that the electron density over EB040 is higher than at AT138. Regarding the periods, it should also be kept in mind that the sounding cadence is 5 min, so the difference is close to the limit of the time resolution. The most important observation from Figure 7 is that the characteristic parameters of the oscillations at both observatories are indicative of MSTIDs, which are most likely triggered by lower atmosphere forcing.

For both AT138 and EB040, it is clear that disturbances are appearing almost at the same time as when the pressure waves in the lower atmosphere arrive. This is different from the data displayed in Figure 6, where the MUF(3000) peaks are seen to appear with a clear delay compared to the pressure waves; this will be further discussed in Section 4.

3.3 MUF from vertical and oblique ionogram traces

Figure 8 shows the MUF(3000) time series for all twelve ionosondes between January 15, 18:00 UT and January 16, 06:00 UT. The MUF(3000) data are offset according to the distances – listed in Table 1 – from the site of the eruption to the various ionosondes. Some ionosondes show much larger disturbances than others. The largest effects are generally seen at the lower latitude observatories: EB040, GM037, VT139, AT138 and EA036. This might be associated with the different travel paths of the tropospheric disturbances, but it could also be influenced by the different background electron densities at different latitudes. This is especially true for the period after midnight when the higher latitude ionosondes – JR055, PQ052, DB049, etc. – show almost flat time series while the lower latitude sounders still show a clear peak – especially obvious in EA036 and EB040. Before midnight there are still clear effects observable at the higher latitudes as well (e.g., for JR055 and PQ052). Note again that the MUF(3000) only exhibits a single peak at the onset of the TID, as was seen already in Figure 6.

Fig. 8 Manually scaled MUF(3000) time series for all ionosondes. The ionosonde closest to the eruption site is JR055 at 15,931 km. The values for the other ionosondes are shifted vertically by 2 MHz per 100 km additional distance from the eruption (see Table 1 for the precise distances). Points are connected by a line if values could be scaled from subsequent ionograms. Whenever for some reason, an ionogram is missing, or the MUF(3000) cannot be reliably scaled, the line is interrupted. The green and orange lines show the arrival times at various distances for the tropospheric acoustic wave, calculated respectively from the barometer data at Dourbes and Ebre. The dashed lines correspond to the waves travelling along the shortest great-circle sector, while the dotted lines correspond to the wave travelling along the longer sector. The grey bands indicate the expected arrival time for the TIDs assuming a delay of 40 min to 1 h from the ground level pressure wave and using the velocities calculated from the Dourbes pressure data for the latter.

The precise timing of the ionospheric disturbances is complicated by the varying time resolutions at different observatories. Only JR055, DB049, AT138, EB040 and EA036 operate using a five-minute sounding cadence, which should provide sufficiently fine time resolution to assure a medium-scale disturbance will be clearly identifiable. At some other observatories, the presence or absence of a clear peak might be affected by coincidence between the TID arrival time and the sounding schedule. Nevertheless, all ionosondes show some disturbance about 1 h after the passage of the lower atmosphere pressure disturbances.

The green and orange dashed lines in Figure 8 indicate the times at which the acoustic gravity waves arrive at each location, assuming they travel with constant velocities of 303 m/s (green) or 300 m/s (orange). Similarly, the dotted lines show the waves travelling along the longer great-circle segments arriving at the same locations (with speeds respectively of 304 m/s and 303 m/s). It can be seen from this figure that the disturbances in the ionosphere appear only sometime after these tropospheric waves. The grey bands indicate the expected interval for the arrival of the ionospheric disturbance, assuming a constant delay of 40–60 min from the arrival times for the tropospheric waves, using the tropospheric wave speeds after Dourbes. These expected times for the appearance of the TIDs indeed correspond fairly well to the onset of the MUF(3000) peaks for those ionosondes where a peak is clearly discernible. Especially for the second wave, many observatories do not show any clear peak at all. Only for the lower latitudes – at EA036, GM037, AT138, and EB040 – is the second peak in the MUF(3000) clearly visible. This might be due to observational difficulties, as discussed in Section 2.2.

The MUF values scaled from oblique ionogram traces are shown in Figure 9. As can be seen from Table 2, the midpoints for the oblique paths PQ052→JR055 and DB049→JR055 are at similar distances from the eruption site. The midpoint of the DB049→EB040 path is somewhat more distant, and the midpoint of EA036→EB040 is the farthest away, at 17,824 km. The limited time resolution and occasional missing data limit the accuracy of TID arrival times that can be determined from these data. Nevertheless, a progression of the peak in MUF can be discerned, with the two closest locations seeing a peak around 20:45 UT, the DB049→EB040 path around 22:15 UT, and the farthest one at 22:50 UT. These times correspond well with those seen from the vertical ionogram data presented in Figure 8. For the EA036→EB040, a second clear peak, can be seen starting at 01:00 UT and reaching a maximum at 01:30 UT on January 16. This agrees well with the appearance of the second wave in the vertical ionogram data from the lower latitude observatories.

Fig. 9 Manually scaled MUF time series from oblique ionosonde traces. Note the different ranges for the vertical axes of the panels. The irregular spacing of the data points for the sounding paths DB049→JR055 and EA036→EB040 is due to the alternating intervals of 5 and 10 min between synchronised ionogram soundings.

The time from the onset of the MUF(3000) increase to its return to the background value is between 50 min and 1 h. Although the limited time resolution makes a precise determination of the start and end of the MUF(3000) increase difficult in some cases, this duration for the event is consistent among all observatories. This includes both the vertical and oblique ionogram data, as well as the second wave, wherever it is visible.

For those time series with a higher time resolution, it seems from Figure 8 that the increasing phase is systematically a little longer than the return to the background value. This can be seen in the JR055 data – 35 min increase and 20 min decrease – in both peaks at EB040 35 min and 20 min for the first wave and 40 min versus 25 min for the second wave, and possibly also the second wave at EA036, although this is less clear because some data are missing. Also, in the oblique data from the EA036→EB040, path this asymmetry can be seen (bottom panel of Fig. 9). However, the time resolution of the data is not sufficient to clearly draw any conclusions on this possible asymmetry.

3.4 TID signatures in individual ionograms

Travelling disturbances can usually be detected from time series of ionospheric characteristics derived from ionograms. However, in some cases, already in a single ionosonde sounding, the effect of a TID can be discerned. For example, in Pradipta et al. (2015), forking, splitting and folding of ionogram traces resulted from the passing of TIDs caused by the 2013 Chelyabinsk superbolide, including at some of the same ionosondes used here. Similarly, we show here some examples of ionograms exhibiting TID signatures in Figure 10. These comprise only a limited number of representative examples from the EB040 and AT138 observatories. Additional ionograms, including from other observatories, can be found in the Giro repository.

Fig. 10 On the left, a sequence of observations from the EB040 ionosonde, including ionograms for 23:00:01 UT, 23:05:01 UT, and 23:10:01 UT. On the right is a sequence of ionograms from AT138. These ionograms are for times 21:00:00 UT, 21:05:00 UT, and 21:10:00 UT. In the ionograms, red points represent O-polarised echoes and green points represent X-polarised ones. Other colours indicate oblique reflections. All panels show ranges from 80 km to 800 km on the vertical axis and frequencies from 0.0 MHz to 6.0 MHz on the horizontal axis. The echoes in the EB040 ionograms starting around 3 MHz and above 440 km are the oblique trace received from the EA036 transmissions, while those at ranges above 600 km are the oblique signals from DB049.

The sequence of three consecutive vertical incidence ionograms recorded by the Athens Digisonde at five-minute intervals, from 21:00 UT to 21:10 UT, is indicative of the tilted ionospheric layers that resulted in O- and X-traces appearing closer together than expected. This is associated with multiple reflections of the sounding signals from irregular layers, several of them having an oblique arrival angle. For instance, in the ionogram for 21:10:00 UT, on the bottom right of Figure 10, both the O- and X-traces are shown to be received obliquely from the north. This indicates a tilted ionospheric layer. In the ionogram of 21:05 UT (middle row on the right), it can be seen that both the O- and X-traces exhibit an unusual kink (between 3.0 and 3.5 MHz in the O-trace). The ionograms recorded in Ebre at five-minute intervals between 23:00 UT and 23:10 UT show split and folded traces which are also indicative of passing MSTIDs.

3.5 Plasma drift velocity observations

In Figure 11, the vertical plasma drift component v_z on January 15, obtained from Digisonde Drift Measurements soundings at DB049 and AT138, are shown, while Figure 12 shows the same for the EB040 observatory for both January 15 and 16. Note that we only discuss v_z here because the behaviour under quiet conditions is known for this component (Kouba & Knížová, 2016), allowing the identification of irregularities; effects of zonal winds, measured closer to the eruption site, were reported in Harding et al. (2022). It can be seen in Figure 11 that the v_z component at both AT138 and DB049 starts showing quick variations from around 17:00 UT. Comparing this to the vertical and oblique MUF and iso-ionic contour plots shown in Sections 3.1–3.3, the v_z variations can be seen to start around the time of arrival of the tropospheric pressure disturbance and the Lamb wave signature, but before the big peak seen in the MUF.

Fig. 11 Vertical drift observed by the ionosondes AT138 (top panel) and DB049 (bottom panel) on January 15. DDM soundings were performed every 5 min, but occasionally no usable data was obtained. This is indicated by gaps in the line.

Fig. 12 Vertical drift from the EB040 ionosonde for January 15 and 16. Similar to Figure 11, gaps in the line indicate skymaps from which no vz could be reliably obtained. Note the disturbances with long period in the early hours on January 15, likely due to large-scale TIDs of auroral origin.

The v_z at EB040 shows a continuous disturbance throughout the period from around 20 UT on January 15 until after 04 UT the next day. The continuous disturbance is again consistent with the tropospheric pressure data shown in the bottom panel of Figure 6, but the disturbances in v_z continue after the disturbances in the pressure have already subsided.

Both in Figures 11 and 12, large vertical drifts can also be observed in the early hours of January 15. These are likely the signatures of the geomagnetic disturbances during the night preceding the eruption. In particular, in the data from EB040 (Fig. 12) it can be seen that these variations have a period of about an hour. The variations seen in the evening hours have shorter periods, indicative of medium-scale TIDs.

3.6 Detrended TEC

Figure 13 presents the latitudinal distribution of the dTEC defined in equation (1) in a range from 30°N to 70°N in the time interval from 18:00 UT to 24:00 UT on January 15, 2022 by considering all the Ionospheric Piercing Points (IPPs) in a longitudinal range from 12.5°E to 17.5°E. Only values of dTEC > 0.1 TECu are shown. Solid black lines represent wavefronts propagating at 310 m/s, highlighting the presence of waves propagating southward from about 19:00 UT to 21.30 UT and northward after 22:00 UT. The former is compatible with an ionospheric perturbation originating from the pressure wave generated by the explosion travelling northward from the eruption, over the North Pole, and finally reaching the European sector from the north about 15 h after the eruption. The second disturbance is related to the pressure wave travelling first southward, reaching the European sector and going northward 18 h after the event. The width of the blue and red bands showing the wave-like perturbations suggests that the southward wave has a smaller wavelength compared to the northward one, as further discussed below. The presence of waves returning from the antipode further confirms that the identified TIDs are not of auroral origin, i.e., not related to the geomagnetic storm occurring around the considered period.

Fig. 13 Keogram reporting the latitudinal distribution of the dTEC from 30°N to 70°N for the time interval 18:00 UT to 24:00 UT on January 15, 2022. To guide the reader’s eye, the solid black lines indicate the wavefronts and the related velocity as inferred from the dTEC in the figure.

Figure 14 shows the hodocrone (i.e., the travel time diagram) of dTEC over Italy, from 16,800 km to 18,500 km from the Hunga volcano as a function of time after the explosion. The solid and dotted black lines represent the theoretical time of arrival calculated by considering waves propagating at the velocity reported in the legend southward and northward, respectively. Dashed and dot-dashed black lines represent the corresponding error band of ±5%. The southward wavefronts travel away from the site of the eruption when passing over Italy, while northward wavefronts travel towards the eruption when passing over Italy. This behaviour is reflected in the opposite slope of the theoretical curves represented by the various black lines in Figure 14. In the case of the first (southward travelling) wavefronts, the measured time of arrival represented by the red and blue bands fits with the theoretical time of arrival within the 5% error band. For the second (northward travelling) wavefront, the measured time of arrival is shorter than the calculated one.

Fig. 14 Hodocrone reporting dTEC over Italy from 16,800 km to 18,500 km from the Hunga volcano as function of time since the eruption (13–19 h). The solid and dotted black lines represent the theoretical time of arrival calculated by considering waves propagating at the velocity reported in the legend northward and southward, respectively. Dashed and dot-dashed black lines represent the corresponding error band of ±5%.

3.7 Swarm C measurements

Figure 15 shows the time profiles of Swarm C Ne (black solid line), Ne-trend (black dashed line), Ne-noise (red dashed line) and dNe (red solid lines) for the periods around 16.7 h (Fig. 15a) and 18.2 h (Fig. 15b) after the explosion. The values of dNe, dNe-trend and dNe-noise (i.e. higher frequency/smaller spatial scales oscillations) are those obtained through the application of Fast Iterative Filtering techniques, according to equation (3). The timing of the main dNe peaks is also indicated in the plots. The dNe behaviour illustrates an intensification of wave-like perturbation after 20:54 UT (Fig. 15a) and after 22:28 UT (Fig. 15b). The timing of the dNe peaks allows estimating the wavelength of the wave-like structure projected along the Swarm C track. Given that, v_c = 7.2 km/s >> v_sound, we can calculate the wavelength by taking into account: (i) for the first Swarm passage, by averaging the wavelength estimated by the time difference between the two maxima and the two minima; (ii) for the second passage, by estimating the distance between the maximum and the minimum. According to this, the first wavelength equals 360 ± 44 km, and the second is 418 km. This agrees with what is seen in the keogram (Fig. 13), in which the blue-red bands representing the southward perturbations appear narrower than those for the northward perturbations.

Fig. 15 Time profiles of Swarm C Ne (black solid line), Ne trend (black dashed line), Ne noise (red dashed line) and dNe (red solid lines) for the periods around 16.7 h (a) and 18.2 h (b) after the eruption. The longitudinal sectors of the two Swarm C passages are 28.9°E and 5.5°E, respectively. The timing of the main dNe peaks is also indicated in the plots.

In addition, we can estimate the mean speed of the propagation of the wave-like structure to the Swarm track by considering the distance between the largest dNe peak (in absolute value) and the Hunga volcano. This results in an estimated velocity of 300 m/s for the first and 340 m/s for the second track. These values are again in agreement with the theoretical expectations.

4 Discussion

It is known that during earthquakes, vertical displacements of the ground or of the ocean floor induce perturbations in the atmosphere and ionosphere. The Rayleigh surface waves generated by earthquakes propagate along the Earth’s surface and induce acoustic waves that 8–9 min later can be observed in the ionosphere (Astafyeva, 2019). Zel’dovich & Raizer (2002, pp. 464) showed that earthquake-induced disturbances, being of acoustic origin, are N-shaped, with an initial overpressure half-cycle with a relatively fast rise time and a slower pressure decay followed by a half-cycle of rarefaction. The observed perturbations were suggested to be the manifestations of long-period ducted acoustic-gravity waves emitted into the ionosphere close to the epicentre. Chum et al. (2012) on Continuous Doppler Sounding System records observed ionospheric disturbances over the Czech Republic excited by the 2011 Tohoku earthquake. The authors pointed out that individual wave packets recorded on the ground had different observed horizontal velocities and corresponded to different types of seismic waves.

On the other hand, tsunamis propagating along the ocean surface generate internal gravity waves that propagate obliquely upwards. Because of the low vertical velocity of about fifty meters per second, these gravity waves reach the ionospheric heights about 45–60 min after their generation on the surface. The tsunami-related ionospheric disturbances are usually characterised as quasi-periodic structures with typical periods between 10 and 30 min, falling in the MSTID range. Such ionospheric disturbances match the period, velocity and propagation direction of the tsunamis from which they originate.

Shults et al. (2016) used data from ground-based GNSS receivers and observed quasi-periodic ionospheric TEC oscillations following the two Calbuco volcano eruptions in April 2015 in southern Chile. The TEC response was registered about 15 min after the beginning of the first eruption and about 40 min after the second eruption. The authors explained such a time delay in ionospheric responses by different source waves emitted by the eruptions. Most likely, the first eruption was accompanied by a shock acoustic wave, followed by the gravity waves generated by the ash emission. During the second eruption, only an ash plume was emitted, producing only a more delayed response in the ionosphere. The apparent velocities of the observed ionospheric disturbance were in the range of 900 to 1200 m/s. It was also noted that the amplitude of ionospheric disturbances seems to scale with the intensity of volcanic eruptions, as they do for earthquakes. TEC changes caused by the five different volcanic eruptions that occurred between 2004 and 2015 – the Asama, Shin-Moe (two eruptions), Sakurajima, and Kuchinoerabu-jima volcanoes – analysed by Nur Cahyadi et al. (2021) showed similar N-shaped disturbances with periods of about eighty seconds propagating outward with the acoustic wave speed in the F region of the ionosphere. The authors believe that such a uniformity suggests its origin in the atmospheric structure rather than in the characteristics of the volcanic eruptions. Infrasound records observed by ground sensors associated with explosive volcanic eruptions have more power in periods much shorter than 1.3 min. However, only those with periods of 1.3–4.0 min can reach the ionospheric F region without strong attenuation (Blanc, 1985). The amplitudes of the detected impulsive TEC disturbances were found to be a few percent of the background absolute vertical TEC, in each case appearing 8–10 min after the eruption.

Saito (2022) analysed TIDs observed over Japan, about 7800 km away from the eruption, using GNSS receiver network data after the eruption of the Hunga Tonga-Hunga Ha’apai volcano. Two types of TIDs with different characteristics were reported as affecting the TEC. The first TID arrived in Japan about 3 h after the eruption. The amplitude was about ±0.5 TECU, the wavelength was estimated as 400 km and the velocity at 250 m/s. The second TID arrived about 7 h after the eruption. The recorded amplitude was about ±1.0 TECU and velocity at 270 m/s. However, the wavelength was longer than in the case of the first TID and was estimated as 800 km. The second-time derivatives of ten-minute interval images provided by the Himawari-8 satellite analysed by Otsuka (2022) clearly showed tropospheric waves propagating at about 310 m/s. The signals in the satellite images well matched the surface pressure observations in Japan. Kulichkov et al. (2022) analysed various characteristics of acoustic-gravity waves induced by the eruption detected at different infrasound stations of the Infrasound Monitoring System and by a network of low-frequency microbarographs in the Moscow region. Using the correlation analysis of the signals at different locations, six arrivals of disturbances emanating from the volcano, which made up to two revolutions around the Earth, were detected.

Fedorenko et al. (2013) showed that TID parameters such as amplitudes, horizontal spatial periods and the TID front inclination angle in the vertical plane increase as the distance between the gravity wave and the excitation source increase. They validated their model using literature data on disturbances generated by about 20 surface and high-altitude nuclear explosions, 2 volcanic eruptions and 1 earthquake, as well as by energetic proton precipitation events in the magnetospheric cusp of the northern hemisphere.

In summary, the literature shows that powerful impulsive events such as this volcanic eruption produce simultaneously, through different mechanisms, ionospheric disturbances of various kinds and with various propagation characteristics. It is, therefore, of great interest to compare the TID signatures observed through the various independent methods presented in Section 3, as some observation techniques can reveal traces of different kinds of disturbances than others.

The arrival times of the disturbances determined from vertical and oblique ionosonde soundings and TEC measurements all agree with the predictions based on a single acoustic wave in the troposphere travelling around the globe. The velocities derived for the disturbances in the lower atmosphere from the pressure data in Dourbes and Ebre, between 300 m/s and 305 m/s, are also in good agreement with each other, as well as with values found in the literature (Haaser et al., 2017; Astafyeva, 2019). Although the period under consideration here was affected by some geomagnetic disturbances, we can therefore be reasonably sure that these medium-scale disturbances were indeed the result of the eruption of the Hunga volcano. Closer to the site of the eruption, disturbances can immediately propagate vertically to the ionosphere and spread radially at ionospheric altitude. However, these disturbances dissipate within a few thousand kilometres from the eruption and do not propagate to the long distances involved in this study, while Lamb waves and infrasound propagating in the lower atmosphere propagate to larger distances (Themens et al., 2022).

The TIDs observed in the MUF, both from vertical and oblique traces, are delayed from the arrival times of the tropospheric gravity wave, as seen in Figure 8. The delays are somewhat longer than the typical 40–60 min mentioned, for instance, in Astafyeva (2019) for TIDs produced through violent disturbances in the lower atmosphere. One possible reason is that the disturbances arrived in Europe at night when the F-layer electron density peak was at the highest altitude. The disturbance, therefore, takes longer to travel from ground level up to hmF₂. Around the same time as the MUF peak, indicators for medium-scale TIDs have also been observed directly in ionograms and skymaps. This confirms that the MUF signature is indeed associated with an MSTID. The length of this peak varies between observatories but is, in all cases, longer than the period of the pressure waves in the lower atmosphere. This is likely due to changes in the wave structure while travelling upwards through an increasingly thin medium.

The arrival times of TID signatures seen in the vertical plasma drift obtained from DDM soundings coincide with those seen in the iso-ionic contours and TEC. However, v_z disturbances were seen to persist after the sea-level pressure wave have passed entirely but before the main peak in MUF. Thus, we conclude that there is a wave propagating upwards from the troposphere, which causes the largest MUF variation with some delay compared to the pressure wave (while, on the other hand, the ionospheric disturbances caused by the Lamb wave coincide with it).

As mentioned in Section 3.1, the waves in the lower atmosphere have frequencies below the acoustic cutoff. Nevertheless, acoustic waves with longer periods can propagate through the upper atmosphere and ionosphere. The earlier onsets of the disturbances seen in the iso-ionic contour plots presented in Figure 7 are most likely due to acoustic waves generated by the lower atmosphere disturbance since gravity waves take longer to travel up to the ionosphere. Acoustic waves can reach ionospheric altitudes in less than 5min. Since the shortest ionosonde sounding cadence is 5 min, the delays for acoustic waves will not be discernible in ionosonde data.

It should be noted that the speed of sound in the lower atmosphere depends strongly on the weather conditions such as temperature, air pressure and background winds, which can be very different along different paths of propagation around the globe. As can be seen from Table 1, the azimuths of the great circles connecting respectively Dourbes and Ebre to the location of the eruption are 17.2° and 13.7°. This indicates that the propagation paths of the acoustic waves to these observatories are relatively close to each other. The azimuths to other ionosondes are more different, with AT138 at 62.3° and EA036 at −6.8°. The paths connecting these observatories to the eruption location are thus very different, and can therefore be expected to have different average sound velocities as well. Similarly, the atmospheric conditions through which the waves propagate upwards are not identical over the entire region. This explains some of the differences in the TID delay times in Figure 6. Harrison (2022) found a series of six pulses in the ground level air pressure observed in the United Kingdom, with the odd and even pulses travelling respectively at speeds of 309 m/s and 314 m/s. This small difference in speed of the disturbances propagating along the two sectors of the great circle can also be seen in the GNSS TEC data we show in Section 3.6.

A clear agreement between the different data sets is observed in the propagation directions of the TIDs. In particular, Figures 8 and 9 for the ionosonde data and Figures 13 and 14 for the TEC show a first circular wavefront contracting towards the antipode of the eruption, followed by a second wavefront seemingly expanding from the antipode.

To verify the correspondence between the signatures in the GNSS and Swarm C measurements, Figure 16 provides a comparison between dTEC and dNe. Specifically, Figure 15a shows the time profile of dNe (blue) from Swarm C, dTEC_att (black dots) and the corresponding spline fitting curve (red line) between 22:24:53 and 22:29:35 UT on January 15, 2022, i.e. the second considered track. Inspired by what was reported in Spogli et al. (2021), the values of dTEC_att are the along-the-track (att) values of dTEC obtained by considering the dTEC values having the shorter spatio-temporal distance with each Swarm C measurement.

Fig. 16 (a) Time profile of dNe (blue) from Swarm C, dTECatt (black dots) and corresponding spline fitting curve (red line) between 22:24:53 and 22:29:35 UT on January 15, 2022. (b) Map of dNe (blue line) from Swarm C on top of dTEC (with dTEC > 0:1 TECu) red and blue points in the same time interval of panel a. Deviations from dot-dashed line represent the reference (dNe = 0 cm−3) is proportional to the positive (right) and negative (left) dNe values.

Figure 16b shows the map of dNe (blue line) from Swarm C on top of dTEC (with dTEC > 0.1 TECu) red and blue points in the same time interval of Figure 16a. Deviations from the dot-dashed line representing the reference (dNe = 0 cm⁻³) are proportional to the positive (right) and negative (left) dNe values. We did not apply the same analysis to the first track, as it passes in a sector scarcely covered by GNSS observations (see Fig. 4 and 5). Notwithstanding the different geometry and nature of the observations, dTEC and dNe are also in agreement after 22:28 UT, i.e., the time after which the wave-like perturbations are enhanced.

5 Conclusions

Atmospheric pressure waves from the eruption of the Hunga volcano travelled around the globe. These waves were still sufficiently powerful when arriving in Europe to produce medium-scale travelling ionospheric disturbances. Signatures of these MSTIDs were detected in the data from different instruments, all giving mutually consistent values for the arrival times, wave periods, and propagation velocities. The ionospheric disturbances were detected with a delay of about one hour after the pressure signatures observed at ground level. This is consistent with gravity waves travelling upwards from the troposphere (Astafyeva, 2019). Therefore, we can be confident that the detected TIDs are indeed associated with pressure waves produced by the eruption.

A significant difference was observed between the onset times of the disturbances seen in the MUF series, described in Sections 3.1 and 3.3, and the iso-ionic contours, shown in Section 3.2 and TEC in Section 3.6. The MUF shows sensitivity to gravity waves, while in the iso-density contours, the signatures of acoustic waves can be detected. This illustrates a major advantage of combining multiple types of data in order to study different facets of an event like this volcanic eruption.

This is the first volcanic eruption for which medium-scale travelling ionospheric disturbances have been observed in such detail and through a variety of complementary instruments on the other side of the world. This illustrates the continuing importance of studying TIDs and developing techniques for their detection and characterisation, as major natural events anywhere on the planet can globally affect the ionosphere.

In recent years, significant progress has been made in developing methods for the automated detection and characterisation of travelling ionospheric disturbances. Several such systems are now working in real-time (Belehaki et al., 2020). Different approaches are used for this purpose, based for instance on observations from ionosondes (Kutiev et al., 2016; Reinisch et al., 2018; Altadill et al., 2020) or using GNSS-based TEC data (Hernández-Pajares et al., 2006; Borries et al., 2017). However, most of these techniques are optimised for the detection of large-scale TIDs and sometimes also assume a planar wavefront. Although we could find clear signatures of MSTIDs in various data sets, this required careful human analysis. The results from several automated TID detection systems can be found online at the URL https://techtide-srv-pub.space.noa.gr/techtide, including archived results for the 15th and 16th of January 2022. Indeed, these archives show that no MSTIDs were identified by these automated systems.

It is therefore important in the future to consider observational methods that can potentially be used to automatically detect also medium-scale TIDs to supplement the existing systems for large-scale disturbances. One such technology is based on continuous Doppler sounding (Laštovička & Chum, 2017). An example of the results of such a system can be seen in Figure 17. Here, a clear signature of an MSTID can be seen precisely at the time it is expected. However, at this time, there is only one such observatory in Europe, in Czechia. Thus, this technology can not be used to monitor the ionosphere over the whole region unless a more comprehensive network of observatories is implemented. For this reason, we did not systematically include these data in the current paper.

Fig. 17 Doppler shift spectrogram recorded in Czechia from 20:30 to 21:30 UT on January 15 at about 4.65 MHz (the transmission frequencies differ by a few Hertz). Each trace is the signal from a different transmitter obtained by the same receiver, with slightly different reflection points.

We can also draw from this event some lessons about the use of ionosonde data for studying MSTIDs. It can be seen from the data shown in Figures 8 and 9 that clearly detecting medium-scale TIDs requires a high cadence of soundings, preferably producing an ionogram every 5 min. From the latter figure, it is also evident that oblique ionogram traces are a valuable source of data and can serve as a virtual observatory midway between two ionosondes. Systematically synchronising the operations of different ionosondes to produce oblique traces can, therefore, significantly increase the coverage of a region. Since medium-scale TIDs are often the result of disturbances in the lower atmosphere, it can also be useful to systematically install pressure sensors collocated with ionosondes. Finally, it is worth noting that the effects of the MSTIDs can be seen in ionosonde data outside the standard Ursi parameters. In Section 3.4, we showed some examples of MSTID signatures in the shapes of ionogram traces and in Section 3.5, we showed how the TIDs manifest in drift observations produced by ionosondes. These sources of data are currently not systematically exploited but might, in the future, prove valuable tools for detecting and studying MSTIDs. The interpretation of DDM-sounding data is currently not as well understood as that for traditional ionogram soundings, making it difficult to draw clear conclusions from these data alone. However, in this case, we have seen several feature present simultaneously in skymap and ionogram data, giving us confidence in our conclusion that these are indeed the signatures of MSTIDs originating in the Hunga eruption.

Overall, this study demonstrates that different types of data from independent network instruments can be combined and compared to obtain a coherent picture of disturbances passing over Europe. Despite the moderately disturbed geomagnetic background conditions, combining data from different observatories allows to clearly identify those TIDs associated with the volcanic eruption.

Acknowledgments

Raw ionograms and manually corrected ionospheric characteristics are available from the Giro repository: http://spase.info/SMWG/Observatory/GIRO. This paper includes some ionospheric data from the USAF NEXION Digisonde network, the NEXION Program Manager is Annette Parsons.

The study is partially funded by the PITHIA-NRF Horizon 2020 Grant Agreement 101007599 of the European Commission and uses resources deployed in the frames of the TechTIDE Horizon 2020 Grant Agreement 776011.

A. Belehaki acknowledges the financial support provided by the AFRL grant award FA9550-19-17019. The Swarm IPIR data can be obtained through the official Swarm website ftp://swarm-diss.eo.esa.int and through the Swarm IPIR web page http://tid.uio.no/plasma/swarm/IPIR_cdf/.

The Fast Iterative Filtering code is available at http:// www.cicone.com. The authors are grateful to Antonio Cicone (University of l’Aquila, Italy) for his valuable help with the FIF algorithm and useful discussions.

GNSS data are provided by the National Land Survey Finland NLS (https://www.maanmittauslaitos.fi/en/maps-and-spatial-data/positioning-services/rinex-palvelu), SWEPOS Sweden (https://swepos.lantmateriet.se/), Norwegian Mapping Authority (Kartverket, https://ftp.statkart.no/), and by the “Archivio Dati GNSS Centralizzato” of INGV, which collects GNSS data from eighty networks in the Euro-Mediterranean region and Africa. The editor thanks Chris Scott and an anonymous reviewer for their assistance in evaluating this paper.

References

All Tables

All Figures

	Fig. 1 Top: real-time Dst provided by the Kyoto WDC for Geomagnetism; Bottom: local magnetic K-index from the Dourbes observatory co-located with the ionosonde DB049 at 50.1°N, 4.6°E. The vertical red line indicates the moment of the explosive eruption, while the purple band shows the interval during which volcano-induced MSTIDs are detected over Europe.
In the text

	Fig. 2 Location of the 12 ionosondes used in this investigation, as well as the 4 oblique sounding links.
In the text

	Fig. 3 Two skymap observations from January 15 produced by the EB040 Digisonde. The left image shows the result of the sounding at 23:04 UT and the right one at 23:14 UT. Each mark represents the direction of a recorded echo, with the colour indicating the Doppler shift. The arrows drawn on the image indicate the automatically calculated bulk plasma drift.
In the text

	Fig. 4 Location of the GNSS receivers used for detrended TEC analysis.
In the text

	Fig. 5 Location of the selected Swarm C tracks. Colour code reports the hours after the eruption.
In the text

	Fig. 6 Air pressure (black lines), MUF(3000) data (red lines), and monthly median MUF(3000) for January (orange line) for Dourbes (top) and Ebre (bottom) for the period from 12 UT on January 15 to 12 UT January 16. Note that there are some gaps in the MUF(3000) data, especially at Dourbes, due to unscalable ionograms.
In the text

	Fig. 7 Detrended iso-ionic true height dh(f) plots for a geomagnetically quiet interval – 2022-01-12 18:00 UT to 2022-01-13 06:00 UT – are presented in the left column for the Athens (AT138, top) and Ebre (EB040, bottom) Digisonde stations. The corresponding plots for the period of interest here, from 2022-01-15 18:00 UT to 2022-01-16 06:00 UT, are shown in the right column.
In the text

Fig. 8 Manually scaled MUF(3000) time series for all ionosondes. The ionosonde closest to the eruption site is JR055 at 15,931 km. The values for the other ionosondes are shifted vertically by 2 MHz per 100 km additional distance from the eruption (see Table 1 for the precise distances). Points are connected by a line if values could be scaled from subsequent ionograms. Whenever for some reason, an ionogram is missing, or the MUF(3000) cannot be reliably scaled, the line is interrupted. The green and orange lines show the arrival times at various distances for the tropospheric acoustic wave, calculated respectively from the barometer data at Dourbes and Ebre. The dashed lines correspond to the waves travelling along the shortest great-circle sector, while the dotted lines correspond to the wave travelling along the longer sector. The grey bands indicate the expected arrival time for the TIDs assuming a delay of 40 min to 1 h from the ground level pressure wave and using the velocities calculated from the Dourbes pressure data for the latter.

In the text

	Fig. 9 Manually scaled MUF time series from oblique ionosonde traces. Note the different ranges for the vertical axes of the panels. The irregular spacing of the data points for the sounding paths DB049→JR055 and EA036→EB040 is due to the alternating intervals of 5 and 10 min between synchronised ionogram soundings.
In the text

Fig. 10 On the left, a sequence of observations from the EB040 ionosonde, including ionograms for 23:00:01 UT, 23:05:01 UT, and 23:10:01 UT. On the right is a sequence of ionograms from AT138. These ionograms are for times 21:00:00 UT, 21:05:00 UT, and 21:10:00 UT. In the ionograms, red points represent O-polarised echoes and green points represent X-polarised ones. Other colours indicate oblique reflections. All panels show ranges from 80 km to 800 km on the vertical axis and frequencies from 0.0 MHz to 6.0 MHz on the horizontal axis. The echoes in the EB040 ionograms starting around 3 MHz and above 440 km are the oblique trace received from the EA036 transmissions, while those at ranges above 600 km are the oblique signals from DB049.

In the text

	Fig. 11 Vertical drift observed by the ionosondes AT138 (top panel) and DB049 (bottom panel) on January 15. DDM soundings were performed every 5 min, but occasionally no usable data was obtained. This is indicated by gaps in the line.
In the text

	Fig. 12 Vertical drift from the EB040 ionosonde for January 15 and 16. Similar to Figure 11, gaps in the line indicate skymaps from which no vz could be reliably obtained. Note the disturbances with long period in the early hours on January 15, likely due to large-scale TIDs of auroral origin.
In the text

	Fig. 13 Keogram reporting the latitudinal distribution of the dTEC from 30°N to 70°N for the time interval 18:00 UT to 24:00 UT on January 15, 2022. To guide the reader’s eye, the solid black lines indicate the wavefronts and the related velocity as inferred from the dTEC in the figure.
In the text

Fig. 14 Hodocrone reporting dTEC over Italy from 16,800 km to 18,500 km from the Hunga volcano as function of time since the eruption (13–19 h). The solid and dotted black lines represent the theoretical time of arrival calculated by considering waves propagating at the velocity reported in the legend northward and southward, respectively. Dashed and dot-dashed black lines represent the corresponding error band of ±5%.

In the text

	Fig. 15 Time profiles of Swarm C Ne (black solid line), Ne trend (black dashed line), Ne noise (red dashed line) and dNe (red solid lines) for the periods around 16.7 h (a) and 18.2 h (b) after the eruption. The longitudinal sectors of the two Swarm C passages are 28.9°E and 5.5°E, respectively. The timing of the main dNe peaks is also indicated in the plots.
In the text

Fig. 16 (a) Time profile of dNe (blue) from Swarm C, dTECatt (black dots) and corresponding spline fitting curve (red line) between 22:24:53 and 22:29:35 UT on January 15, 2022. (b) Map of dNe (blue line) from Swarm C on top of dTEC (with dTEC > 0:1 TECu) red and blue points in the same time interval of panel a. Deviations from dot-dashed line represent the reference (dNe = 0 cm−3) is proportional to the positive (right) and negative (left) dNe values.

In the text

	Fig. 17 Doppler shift spectrogram recorded in Czechia from 20:30 to 21:30 UT on January 15 at about 4.65 MHz (the transmission frequencies differ by a few Hertz). Each trace is the signal from a different transmitter obtained by the same receiver, with slightly different reflection points.
In the text