The January 2022 Hunga Volcano explosive eruption from the multitechnological perspective of CTBT monitoring

SUMMARY

The massive eruption of the Hunga Volcano on 15 January 2022 provides an ideal test case for reviewing established methods to discriminate and analyse source processes. Discriminating source mechanisms and identifying their origins is a key task when analysing suspicious events in the frame of the Comprehensive Nuclear-Test-Ban Treaty (CTBT). Earthquakes and explosions can be distinguished in some cases using well established methods such as inversion for the seismic moment tensor. In more complex cases the combination of analyses of the seismic, infrasonic and hydroacoustic waveform content can be of help. More challenging is the discrimination of the specific kind of explosive source such as a nuclear test and a volcano eruption based on the data from the three waveform technologies alone. Here, we apply standard techniques destined to analyse relevant events in the frame of the CTBT, that is all three waveform technologies (seismology, infrasound and hydroacoustic) and atmospheric transport modelling of radionuclides. We investigate the potential of standard analysis methods to discriminate a source and identify their possible weaknesses. We show that the methods applied here work very well to identify, investigate and discriminate an explosive event. During discrimination we could not only exclude a shear-source (i.e. earthquake) but also distinguish the volcanic explosion in contrast to a man-made explosion. However, some tasks remain difficult with the available methods. These tasks include the reliable estimation of the strength of a non-shear event and thereupon a yield estimation of a possibly CTBT relevant event. In addition to evaluating our methods, we could relate our results with specific phases of the eruption process providing a more detailed insight of what happened. Our investigations of the eruption details only provide a starting point for further in-depth analysis. However, they underline the importance of the Hunga eruption event for science. The huge amount of observations provide a unique opportunity for knowledge gain in several subdisciplines of the geosciences. In addition, although not being a nuclear test, it also provides a useful and important data set for further developing multitechnology analyses in the frame of the CTBT.

1 INTRODUCTION

The Comprehensive Nuclear-Test-Ban Treaty (CTBT) bans all nuclear explosions, regardless of purpose and location (CTBT 1996). For fulfilling the compliance with the CTBT, an International Monitoring System (IMS) has been established, equipped with sensors for three geophysical waveform technologies: seismology, hydroacoustics and infrasound (SHI). A fourth technology monitors radionuclides and noble gases in the atmosphere, supported by atmospheric transport modelling (ATM) to connect possible detections with potential source events (Pilger et al. 2017). This approach offers a reliable monitoring in combination with a comprehensive analysis of events, both vital to identify and to reveal potential nuclear tests. On 15 January 2022, the violent eruption of the Hunga Volcano triggered large parts of the global IMS network for waveform technologies. The tremendous strength of this event makes it an ideal natural test case to evaluate and discuss so far established routines and methods in the context of the CTBT’s verification regime. In this paper, we show and discuss the results of analysing data from all three waveform technologies (SHI) and atmospheric transport modelling of radionuclides. By doing so, we try to characterize the eruption as comprehensively as possible.

The Hunga Volcano is situated between the islands of Hunga Tonga and Hunga Ha’apai and is part of the Tonga-Kermadec volcanic arc, located about 65 km north of the main island and capital of the Kingdom of Tonga. Previous to the January 2022 eruption, the submarine volcano showed only weaker eruptive activity in 2009 and during the 20th century (Bohnenstiehl et al. 2013). During a 5-week-long eruptive phase starting in December 2014 (Garvin et al. 2018), the caldera reached the sea surface and a barren volcanic landmass was formed connecting the two neighbouring, uninhabited islands of Hunga Tonga and Hunga Ha’apai. Volcanic activity increased again in mid-December 2021 with continuous, remotely detected infrasound emission from the volcano until early January 2022 (Vergoz et al. 2022). Satellite imagery revealed a significant growth of the tephra cone (Global Volcanism Program 2022). Seismic and hydroacoustic activity was not observed at this stage (although possibly only because the nearest seismological station was more than 700 km away). After a few days of quiescence, infrasonic activity recurred during 13 and 14 January 2022, associated with a major eruption that removed most of the new land from 2014 to 2015. The largest and most violent phreatoplinian eruption of the whole sequence took place around 04:15 UTC on 15 January 2022. Strong seismic, hydroacoustic and atmospheric waves were observed all around the globe and the ash plume of the volcano raised up into the mesosphere to an altitude of almost 60 km (Matoza et al. 2022; Proud et al. 2022).

Probably the most striking feature of the volcanic eruption was a Lamb wave (Lamb 1924), an acoustic-gravity wave travelling at the speed of sound in the lowermost few kilometres of the atmosphere. It was clearly recorded as a worldwide pressure variation in the range of hundreds of Pascals (hPa) with infrasound and barometric sensors, at high accuracy weather stations and even on smartphone devices (Matoza et al. 2022; Díaz & Rigby 2022; Grabowski 2022) leading to resonant acoustic normal modes (Widmer-Schniedrig 2022) and oceanic meteotsunamis (Kubota et al. 2022). The Lamb wave was also detected in satellite observations (Matoza et al. 2022) and identified in infrasonic measurements even after circumnavigating the globe up to eight times (Vergoz et al. 2022).

Among the earliest historic recordings of natural infrasound is data of the volcanic eruption of Krakatoa volcano in 1883 (Symons 1888) and the Tunguska meteorite impact in 1908 (Whipple 1930), where absolute air pressure changes of these events were identified in precise barometers, providing resolutions of some Pascal and some tenths of Pascal, respectively. Both signals circumnavigated the globe multiple times and also generated pressure changes in the order of hPa (Matoza et al. 2022). Further strong volcanic eruptions, accompanied by significant acoustic waves recorded as air pressure variations and far-reaching infrasound signals are reported from Mount St Helens 1980 (Donn & Balachandran 1981), El Chichôn 1982 and Pinatubo 1991 (Widmer & Zürn 1992). The strongest man-made pressure changes, documented in barometric records, are the nuclear tests starting in 1945 (Donn & Shaw 1967), which culminated in the Tsar bomb detonation in 1961, a hydrogen aerial bomb with an explosive yield of approximately 58 Mt of TNT equivalent (Stevens et al. 2002). The Tsar bomb explosion resulted in peak-to-peak pressure changes of up to 2 hPa at distances of 2000–10 000 km, whereas the Hunga Volcano reached changes of up to 5 hPa within the same range (Matoza et al. 2022). First estimates on the yield of its eruption resulted in 100–200 Mt of TNT equivalent (Vergoz et al.2022).

The IMS has several dedicated hydrophone and seismic T-phase stations to detect nuclear underwater explosions. Naturally, also signals from underwater earthquakes and underwater volcanoes are recorded. According to event analyses of the CTBT organization (CTBTO), underwater volcanism is especially prevalent in the Pacific, where up to several thousand seismic events can be recorded during a single active volcanic phase. Observation points are also located close to previous nuclear weapon test sites, such as the Marshall Islands and French Polynesia (Talandier & Okal 1987). Broad-band seismic records are used to study the time dependent source processes of volcanism such as the 2022 Hunga eruption phase due to their high resolution in time and global station coverage (Zobin et al. 2006; Poli & Shapiro 2022; Yuen et al. 2022). Still, underwater volcanism and its seismic activity can remain undetected, as initially happened during the Mayotte earthquake swarms (Fallou et al. 2020), where very low frequency and monochromatic earthquakes have occurred (Cesca et al. 2020). These kinds of underwater sources, which are not directly visible, make CTBT-relevant source discrimination even more challenging.

The dynamic source processes during volcanic activity in general and underwater in particular can repeatedly excite low-frequency harmonic surface waves (Talandier & Okal 1996). Modelling shows that volcanoes can radiate low (∼1–3 Hz) and high (∼5–10 Hz) frequency signals related to gas bubbles moving through a conduit filled with magma and breaking at the liquid surface, respectively (Ripepe et al. 2001). For underwater eruptions, bubble pulse and water column reverberation can cause spectral scalloping, that is distinct spectral peaks, which can be observed in both hydroacoustic and seismic data (Koper et al. 2001; Ping et al. 2022). Observations of seismic records of underwater explosions are dominated by low-frequency content (<1 Hz; Ping et al. 2022). Hence, underwater explosions as well as low-frequency underwater volcanic eruptions can have similar spectral signatures owing to the fact that the origin of both potentially involves bubble type sources, making discrimination an important and challenging task of manual and automatic data analysis in the CTBT framework (Del Pezzo et al. 2003).

Even though the seismic and acoustic signals of the 2022 Hunga event have been recorded globally (Yuen et al. 2022; Poli & Shapiro 2022; Matoza et al. 2022), this eruption is not related to exceptionally strong seismic magnitudes (see Table 1). However, the continuous volcanic activity has been seismically observed to last about 12 hr (Yuen et al. 2022). In general, standard seismic methods and concepts such as magnitude-yield relations do not apply to volcanic eruptions. Therefore other metrics are used to describe the force of a volcanic eruption: based on the seismic data records, a volcanic explosivity index (VEI) of about six has been inferred for the main eruptive phase (Poli & Shapiro 2022). This value is comparable to the 1883 Krakatau and the 1991 Pinatubo eruptions, being two of the largest eruptions ever measured (NOAA National Centers for Environmental Information 2022).

In this study, we focus on three distinct and impulsive events with clearly identifiable P-phase onsets in the seismological records allowing the recognition of the three events as individual particular sources amidst the complex source processes of the Hunga eruption. We identify the events, as described in the data and method section, as a smaller event at 04:01 UTC, the strong event at 04:15 UTC and another strong event at 04:18 UTC, which is the first of a sequence of events between 04:18 UTC and 04:21 UTC. We consider all four CTBT-relevant technologies for interpreting the complex source processes.

2 DATA AND METHOD DESCRIPTION

Within this paper we concentrate on those volcanic processes of the Hunga eruption, which are related to the main eruptive phase on 15 January between 00:00 and 12:00 UTC. Both, the catalogues of the European-Mediterranean Seismological Centre (EMSC) and the United States Geological Survey (USGS) only list an eruption at 04:15 UTC. Whereas the Incorporated Research Institutions for Seismology (IRIS) catalogue itemizes six events with m_b ≥ 4.0 until 12:00 UTC, two before (04:07 and 04:13 UTC) and three after the main eruption phase (04:40, 05:30 and 09:51 UTC). We skip the 09:51 UTC event because of its depth of about 80 km. The catalogues of the USGS and IRIS have been merged at a later stage. The Reviewed Event Bulletin (REB) provided by the International Data Centre (IDC) of the CTBTO is the most comprehensive event bulletin for Hunga’s eruptive period as it accounts for all three waveform technologies to form an event: seismology, hydroacoustics and infrasound. A summary can be found in Table 1, where we distinguish the REB entries between the considered technology as some are infrasound-only or seismology-only events.

2.1 Seismology

2.1.1 Seismic data and array processing

We use three-component waveforms of stations up to an epicentral distance of 10 000 km (90°, see Figs 1 and 2). We use two different seismic data sets, one obtained via FDSN service from IRIS and GEOFON (Scripps Institution of Oceanography 1986; Hanka & Kind 1994) and one considering only the IMS seismological station data.

The data are converted to displacement by removing the instrument response. The horizontal components are rotated from the geographical coordinate system (ZNE) into the event-centric coordinate system (ZRT). For inverting moment tensors and source time functions, the data are resampled to 0.5 Hz to match the calculated Green's function stores sampling rate. It is filtered with a 4th order Butterworth bandpass between 0.01 and 0.09 Hz. For other applications used in this paper, the frequency range might change slightly (see corresponding method sections).

Furthermore, we use the stations of the German Regional Seismological Network (GRSN) and the Gräfenberg Array (GRF) in southern Germany for array analyses (see Supplement 1 for a list of stations used). Currently, these networks consist of 48 and 13 permanent broad-band stations, respectively. The epicentral distance of both networks to the Hunga volcano varies slightly around 16 700 km (150 ± 3°). This distance is in the centre of the range of 140° to 160° where the different branches of the PKP caustic arrive (PKPdf, PKPbc and PKPab). We perform a frequency–wavenumber (fk) analysis using two different simulation filters, a short-period WWSSN-SP and a long-period SRO-LP filter. Usually, the WWSSN-SP filter (response frequency ∼1 Hz) is used for such an analysis because it emphasizes the PKP energy radiated by teleseismic earthquakes neatly. In the Tonga case, the SRO-LP filter (response frequency ∼0.05 Hz), designed to analyse surface waves, performed superior to the WWSSN-SP filter due to the signal including lower frequencies. Applying a beam-forming method, we estimate slowness and azimuth of the incoming phase and use them to locate the event. Considering GRSN stations and applying the SRO-LP simulation filter we determine a surface wave magnitude M_S for the event at 04:15 UTC of the eruption.

Moreover, we calculate spectrograms of the vertical ground velocity for selected stations (see Supplement 2 for station list) up to 5000 km (45°) distance (Fig. 3). Using this selection we want to ensure both, high signal-to-noise ratios and avoiding potentially dominating local earthquakes. For the spectrograms, the power spectral density (PSD) is calculated in sliding time windows of 51.2 s (2048 samples) overlapping by 70 per cent for a time range between 03:00 UTC and 11:30 UTC. In order to estimate the duration of the main eruption phase, we compute the spectral power (SP, summation of PSD) between 0.02 and 0.09 Hz. Assuming the signal is caused by surface waves, we align these SP curves using a velocity of 3680 m s⁻¹ (Zhang et al. 2022), compensate for geometrical spreading of surface waves (multiplying with distance as power is considered) and linearly stack them resulting in a mean spectral power over all stations.

2.1.2 Waveform inversion for moment tensor

For computing waveform inversions to estimate seismological moment tensors, we create two data sets, one with FDSN stations, including IMS stations which are available via FDSN, and one with IMS stations only. We select stations at epicentral distances between 24° and 93°. Thus, we exclude near-field wave effects and possible phase triplications in the far-field. In total, we consider data from 119 seismological broad-band stations for the FDSN station only data set and 29 stations for the IMS only data set after manual exclusion of noisy records. See Supplement 3 for a list of the IMS stations and Supplement 4 for the FDSN stations used.

We use the open-source software Grond based on the Pyrocko package (Heimann et al. 2017, 2018; Kühn et al. 2020). Grond applies a Bayesian bootstrap-based probabilistic joint inversion scheme. This provides the advantage that uncertainties of the solution are retrieved along with the best-fitting moment tensor. In addition, different input data types can be combined by a flexible design of the objective function. Green’s functions are provided by a pre-computed database to facilitate the fast calculation of synthetic waveforms (Heimann et al. 2019). The Green’s function database is calculated based on the AK-135 1-D velocity model (Kennett et al. 1995; Montagner & Kennett 1996) using the QSSP code (Wang 1999). We invert simultaneously for the six moment tensor components, origin time, source time duration and depth. The epicentral location of the events is fixed to 20.546°S and 175.39°W in order to reduce the number of inversion parameters. We determine a first rough guess of the origin time based on manually picked P-arrival times and station distances. We use the same stations and the misfit configuration for all analysed events. Geometric spreading is compensated by applying a distance related weighting scheme to all stations. We assume a triangular source time function and allow for a depth variation between the surface and 3000 m below the seafloor (Shane Cronin personal communication; Colombier et al. 2018).

We invert the P- and SH-phases simultaneously in the time domain between 0.025 and 0.09 Hz, a frequency band which is suited for point-source teleseismic inversions of body waves, as it is between the microseism dominated low frequencies and the corner frequency of the explosion. For the P-phase, we fit only the vertical components in a tapered time window of 94 s around the theoretical phase onset. For the SH-phases, the vertical and transversal components are used in a tapered time window of 262 s around the theoretical phase onset. We allow a shift of the onset time of the source of ±25 s around the estimated origin time.

Bootstrapping for seismological waveforms in the Grond framework is realized through Bayesian random choosing of the station weights (Rubin 1981). The inversion in the time domain and with an L2 norm is performed in 100 parallel so-called ‘bootstrap chains’ for different combinations of misfit functions, where each bootstrap chain has a different set of pseudo-random station weights. Therefore each bootstrap chain is an independent optimization process. We use 5000 uniform distribution sampler iterations, followed by 180 000 iterations using the directed sampler of Grond, where the sampled distribution shrinks in relation to the variance of the best performing models in the bootstrap chains (Kühn et al. 2020).

2.1.3 Inversion for source time function

To separate the source time function from other parts constituting the seismograms, we perform a linear least-squares inversion assuming different point sources (Gaebler et al. 2019). Seismological waveforms can be assumed to be a convolution between source time function s, Green’s function G, noise and site effects. While considering a global data set of waveforms, we assume that local noise and site effects are averaged out and can therefore be neglected. If u(t) is the seismogram at a station with coordinates r_r of a source with coordinates r_s, the P-wave train (i.e. P, pP and sP phase) at time t is represented by

where * is a time convolution and the summation convention applies. M_jk is the moment tensor component who’s force starts at axis k and points in the direction j. G_j, k are the spatial derivatives of the Green’s function indicated by the comma before index k. The displacement in the far-field is controlled by the time derivative of m(t), which is declared as moment rate function |$\dot{m}(t)$|⁠. This explains that the far-field body wave pulses from earthquakes are single-sided pulses (Dahm & Krüger 2014). Because only the far-field is considered, the Green’s functions in eq. (1) are replaced by far-field Green’s functions G^(ff), and m(t) is replaced by |$\dot{m}(t)$|⁠. The inversion is set up for |$\dot{m}(t)$| using a set of K time-shifted triangular basis functions h_l(t):

where β_l is a weighting factor for each h_l(t).

The convolution in eq. (3) can be written in discrete form, leading to an overdetermined matrix system for unknown weighting factors to be solved in a least-squares sense using the L2 norm (Dahm & Krüger 2014). We set up the inversion using a set of point sources with an amplitude of 1 N located at the surface (seafloor) and a pulse width of 4 s. The sources are shifted in increments of 10 s. In a next step, the respective synthetic waveforms are calculated for all vertical components and for each of these sources. The synthetic traces are truncated to 5 s before and 50 s after the theoretical P-wave arrival and zero-padded to the same length.

As source M, we use the inverted moment tensor, an explosion, and a vertical single force. Obviously, an explosion represents the typical scenario targeted in the field of CTBT verification (Gaebler et al. 2019). The single force source is a valid assumption for volcanic eruptions (Kanamori et al. 1984), assuming the bulk of the material involved in the volcanic process can be approximated by a perfect gas (Nishimura 1998). The force acts vertically downward and represents the counter force to the thrust of the volcanic jet. We use stations between 4400 km (40°) and 10 000 km (93°) distance (Fig. 1 and supplement 5 for a list of the stations) and the same Green’s function store as in the inversion for the moment tensor. The inversion is performed in a time window from 03:00 to 07:00 UTC. We align all waveforms with respect to the theoretical P-wave arrivals. To estimate the uncertainty stemming from potential data error, we pseudo-randomly bootstrap the weights between 0 and 1 for all station records 300 times. Therefore, in each inversion run, the contribution of each station varies in a similar fashion to the Bayesian bootstrap used in the inversion for the moment tensor (Rubin 1981).

2.2 Hydroacoustics

2.2.1 Hydroacoustic data and array processing

The hydroacoustic component of the IMS consists of eleven stations (de Groot-Hedlin & Orcutt 2001). Five of them are so called T-phase stations, which are seismometer sites near the coasts of deeper ocean areas. The other six stations, which are the only ones considered in this study, are hydrophone arrays. Each of them consists of a triplet of hydrophones mounted in a water depth of about 1 km (i.e. within the acoustic wave guide described in the next section). All but one of these six stations (HA01) are island stations comprising a northern and a southern triplet to avoid shadowing of a signal due to the island itself.

Waveform data of the hydroacoustic stations is available in a time resolution of 250 Hz. The separate instruments of each hydrophone triplet can be analysed, for example using array processing to derive backazimuth and apparent velocity information of coherent sources within the observations as well as their time-dependent amplitude and frequency content. We apply the Progressive Multi-Channel Correlation Method (PMCC, Cansi 1995) for this task and identify signatures related to the 2022 Hunga eruption, but also to other hydroacoustic sources like earthquakes, in the corresponding array data of two of the six hydrophone stations (HA03 near Juan Fernandez Island and HA11 near Wake Island, see Fig. 1).

2.2.2 Underwater propagation and wave coupling

Hydroacoustic waves most efficiently propagate over large distances within an underwater acoustic duct at approximately 1 km water depth, which is called SOFAR (sound fixing and ranging) channel. It is an acoustic waveguide which is formed due to a minimum of the hydroacoustic sound speed (around 1500 m s⁻¹) in the water column. The actual sound speed value as well as the depth of the SOFAR channel depend on the ambient water pressure (increasing with water depth) and water temperature (decreasing from the surface but stabilizing at a certain depth). This produces a waveguide and hydroacoustic propagation within that waveguide can be modelled using, for example ray tracing or parabolic equation methods (Oliveira et al. 2021). For this study, a direct back-projection from the hydroacoustic stations along the observed backazimuth directions is performed. This is adequate and sufficient to identify a coupling region, most likely along a steep ridge or trench that reaches the necessary water depth of the SOFAR channel, where seismic waves from the volcano are converted into hydroacoustic waves propagating in a SOFAR channel towards the station.

2.3 Infrasound

2.3.1 Infrasonic data and array processing

The infrasound component of the IMS is planned to consist of 60 stations (Christie & Campus 2009), of which 53 were certified and operational in January 2022. Each station is constructed as an array of at least four (and up to 15) microbarometers distributed over one or more square kilometres with individual wind noise reduction systems at each array element. Waveform data of the infrasound arrays is available at a 20 Hz resolution. The frequency range of 0.02–4 Hz is the IMS passband minimum requirement, and in this band the sensor response is considered flat. It covers the core range of interest for explosion monitoring and other anthropogenic sources as well as many natural sources. Similar to hydroacoustics, the waveform data of the IMS infrasound arrays can be analysed using array processing methods like PMCC. Here, we process the data with one-third octave frequency bands between the centre frequencies 0.02 and 5 Hz. The window lengths decrease with frequency from 636 to 25 s, with 95 per cent overlap. For capturing atmospheric wave modes beyond that frequency range, the waveform data had to be restituted by removing the instrumental response (Vergoz et al. 2022).

2.3.2 Atmospheric propagation and origin time inversion

Infrasound can propagate through the atmosphere over long distances if a stratospheric waveguide is present. This waveguide evolves when the effective sound speed (speed of sound plus along-path wind speed) in the upper stratosphere exceeds the speed of sound near the surface. Strong winds at these atmospheric layers (between 30 and 60 km) mainly control this effective sound speed, as their direction changes twice per year. Downwind conditions therefore allow a stratospheric waveguide, in which the attenuation is low compared to upper layers. For the back-projection of signals detected at the IMS infrasound stations, we use modelled celerities (propagation velocities of atmospheric waves over ground) defined on a global grid, the calculation of which is described in detail in the supplement S1 of Vergoz et al. (2022). The infrasound celerities result from 2-D acoustic ray tracing simulations for a full orbit. The Lamb wave celerities were calculated from the average effective sound speed for the first 4.8 km in the troposphere (Vergoz et al. 2022). Since the source location was well known shortly after the eruption, we focus the location procedure on the estimation of the origin time(s), based on the distances and celerities. From the celerity models for infrasound and Lamb wave propagation, the mean celerity along the propagation path (azimuth and distance) is considered while the source coordinates are fixed.

For the Lamb wave, we account for the short-orthodrome arrivals at stations within 90° range. We apply a second-order Butterworth bandpass filter between 1800 and 3600 s to the restituted waveforms and pick the time of the maximum amplitude (Vergoz et al. 2022). Then we determine the period of the Lamb wave and subtract half of the period from the picked time of the Lamb wave’s maximum amplitude. The resulting time is assumed to be the onset time of the Lamb wave. The arrival times of the later infrasound signals are based on PMCC detections (0.02–5 Hz) with backazimuths consistent with the Hunga volcano. Here, we constrain the inversion for the origin time to the infrasound arrivals at stations IS22 (New Caledonia), IS05 (Tasmania), IS07 (Central Australia) and IS04 (Western Australia). These stations are located within the strongest downwind conditions of the westward stratospheric waveguide [see supplement S1 of Vergoz et al. (2022)] and thus we assume comparable propagation paths with low uncertainty.

2.4 Source time functions of all waveform technologies

For investigating the consistency of the source processes recorded by the different technologies, we apply a second method of determining the source time function by linearly stacking the records of the three waveform technologies. Thus, we determine a (quasi-) source time function of hydroacoustics, infrasound and seismology. Here, stacking is a simple and straightforward approach assuming that unwanted path effects are averaged out while SNR is enhanced. However, sufficient data from a large variety of directions is mandatory for obtaining reliable stacks.

For the seismological data stacks we use the approach similar to Houston & Vidale (1994). We select the vertical component traces of stations between 4400 km (40°) and 10 000 km (93°) to ensure coherence between data of different stations and avoid bias due to near-field effects. We create a data set with 109 FDSN stations (see Supplement 4 for the list of stations) and another separate data set with 27 IMS stations (see Supplement 5 for the list of stations). These two separate data sets for stacking the seismic data allow us to assess the necessity for additional data besides the IMS network for investigation of the Hunga eruption. We filter the data with a bandpass filter of 4th order between 0.025 and 0.1 Hz. The traces are aligned based on the slowness of the P wave based on the AK-135 model (Kennett et al. 1995; Montagner & Kennett 1996) and normalized to the maximum before linear stacking and dividing by the number of stations.

For stacking the hydroacoustic records we use the southern hydroacoustic array at Wake Island, H11S and the northern and southern arrays at Juan Fernandez Island, H03N and H03S. Only two stations at H03N can be used, as H03N3 was malfunctioning. The northern array at Wake Island cannot be used for the stacking, because the direct propagation path to the source is blocked by the island itself and therefore the recorded signal potentially includes complicated path effects. The hydroacoustic records are first converted to Pa and filtered with a 4th order Butterworth bandpass between 2 and 13 Hz and aligned assuming a velocity of 1500 m s⁻¹ in a homogeneous medium. We then stack the absolute waveform records at each array individually and normalize them before stacking linearly all resulting three stacks together. We finally divide by the three contributing stacks to normalize the result between 0 and 1. In total eight sensor records are used.

In a similar manner, we linearly stack the envelopes of infrasound records using expected arrival times for the Ig1 phase that are based on the distance and celerities (Vergoz et al. 2022). We perform the stack for three different subsets of stations to investigate the uncertainty associated with the complex atmospheric paths and different phase arrivals. We choose stations which are at maximum 10 000 km (63°) away. One stack is done based on stations in the easterly propagation direction (see Figs 4a and e) where low celerities prevail. For the same reason as for the origin time inversion (see Section 2.3.2), we also perform two stacks of infrasound waveform records of stations in the westerly propagation direction, one with a larger subset of stations (see Figs 4b and f) and one with only four stations (see Figs 4c and g). We filter the data after restitution with a bandpass filter between 0.05 and 0.5 Hz, above the Lamb wave and mountain wave frequencies and below the high frequency band of 0.5–5 Hz dominated by anthropogenic sources. We first stack the envelopes of the individual waveform records at each of the arrays. We then normalize them before stacking linearly all resulting stacks together, dividing by the number of contributing stacks. To visually smooth the resulting function we low-pass filter the envelope of the resulting stacks of all three technologies at 30 s.

2.5 Radionuclide measurements supported by atmospheric transport modelling

For proof of the nuclear origin of an explosion the radionuclide (RN) network of the IMS monitors tiny traces of radioactive fission and activation products in the atmosphere (Pilger et al. 2017). After network completion, 80 IMS RN stations with high volume air samplers collect aerosol particles which are then analysed by high resolution gamma spectroscopy with a germanium detector. As of 2022, 72 RN particulate stations are installed and certified. They are extremely sensitive with minimum detectable activity concentrations in the range of a few μBq per m³ due to the large sampled volume and the long acquisition time.

A subset of the stations is equipped with noble gas systems measuring radioactive xenon isotopes. These isotopes are more likely to escape than particulates from underground nuclear explosions and have convenient tracer properties in the atmosphere as they are not subject to wet and dry deposition or chemical reactions. Minimum detectable concentrations are in the range of a few tenths mBq m⁻³. The only relevant sink is radioactive decay, which makes these isotopes also suitable for event timing and source discrimination using isotopic ratios (Kalinowski et al. 2010). There are also civilian nuclear facilities releasing radioxenon; the most significant background contribution is emitted alongside the production of radioisotopes for medical purposes.

To assess the possible source regions of radionuclide samples atmospheric transport modelling (ATM) is applied. The standard approach is to use Lagrangian Particle Dispersion Models in backward mode to calculate source-receptor-sensitivity (SRS) fields which denote the linear contribution of each grid cell at a certain source time step to the activity concentration in the sample (Becker et al. 2007). If there is some knowledge about a potential release location such as a nuclear facility or a suspicious waveform event, the atmospheric dispersion can also be simulated in forward mode to predict expected concentrations at radionuclide stations. As the IMS RN network aims for the detection and characterization of rather short-lived anthropogenic radionuclides in the atmosphere, it is not generally expected to detect radionuclides from volcanic eruptions. However, in the context of CTBT monitoring the radionuclide component is crucial and so the potential contribution has to be included in a multitechnology analysis of any potential treaty violation.

For simulations of atmospheric transport we use the Lagrangian Particle Dispersion Model HYSPLIT which was developed and is maintained by the Air Resources Laboratory of the US National Oceanic and Atmospheric Administration (NOAA) (Stein et al. 2015). We operate the model driven by NOAA-NCEP (National Centres for Environmental Prediction) meteorological data. For the model runs discussed in this study we chose GDAS (Global Data Assimilation System) meteorological analysis data with 1° horizontal and 3 hr temporal resolution. This resolution has been proven to be sufficient to model long range transport particularly over oceans in an ATM inter-comparison experiment (Maurer et al. 2018). In each model run for our assessment here, a release of one million model particles is simulated and resulting concentrations (or sensitivities, in backward mode, respectively) were sampled to a 0.5° × 0.5° grid.

3 ANALYSIS AND DISCUSSION OF SHI AND ATM RESULTS

Through a first rough inspection of the seismological waveforms we identify three distinct events comprising transient signals with visible and clearly distinguishable P-wave onsets, and one later event with an extremely emergent and hardly identifiable P-wave onset but clearly separated in time by more than 4 hr. They are listed in Table 2 as ‘BGR-SE’. The same events were also identified in IMS hydroacoustic data, complemented by a fifth event at 04:12 UTC; all hydroacoustic events are listed as ‘BGR-HA’. From IMS infrasound data, we obtain four infrasonic events between 04:00 and 08:30 UTC plus the origin time of the Lamb wave at around 04:29 UTC, all listed as ‘BGR-IS’. The estimated uncertainty for the infrasound events is up to 6 min as detailed later. Considering also Table 1, the fact that different catalogues and observers identified such a variety of events clearly illustrates the complexity of the processes related to the eruption. Within this paper, we concentrate on the first three transient events that we identified from seismological data. These are the ones which are more closely related to the main eruption phase. In general, we use data of all three waveform technologies (SHI) as well as ATM for forward and backward tracking of plumes containing radionuclides of interest.

3.1 Localization of the event in space and time

In case of detecting a possibly explosive event, the first processing step is the localization of the event in time and space. This step is usually carried out on determined arrival times of seismological phases and appropriate location algorithms. In the case of the Hunga eruption the coordinates of the origin were very clear from the beginning; for example, due to timely available GOES-17 satellite images showing the impressive ash cloud (NASA Jet Propulsion Laboratory 2022a). Also, the hypocentral depth was clear within the usual range of uncertainty (a few kilometre, depending on station geometry). Regardless, the localization of events with such large magnitudes no longer poses a challenge nowadays. Therefore, we did not determine the coordinates and depth from seismological data by applying classic localization routines again. Just as a by-product of PKP-phase analysis on data of the GRSN array in Germany (16 700 km, 150°), we determined the azimuth and slowness value of this phase to 14.9° and 2.41 s deg⁻¹, respectively, for the 04:15 UTC event and used these to backtrack the ray path. As expected, the signal is assigned to originate from the Hunga volcano within 150 km uncertainty which is a reasonable value.

However, we refined the origin time using various data and methods. For the event at 04:15 UTC, the origin times provided by different catalogues vary in a range of up to 14 s (see Table 1). According to our seismological waveform inversion for the moment tensor, the main shock occurred at 04:14:54 UTC (see Table 2). For the two events for which we also inverted the waveforms for moment tensors, we estimated the origin times to 04:01:40 UTC and 04:18:33 UTC, respectively. There is a discrepancy between IDC’s REB, publicly available event catalogues, and our own event identification. For example, the events at 04:01 UTC and 04:18 UTC are not included because the USGS, IRIS and EMSC seismic catalogues do not list individual volcanic events but only main eruptions. Additionally, the complexity of the physical processes involved during the eruption makes it difficult to distinguish single events due to low body-wave energy and mingling of phases of the single events.

Regarding the hydroacoustic event processing, our findings of Hunga eruption signatures in IMS hydrophone data are in agreement with the REB as well as other published information (Matoza et al. 2022; Le Bras et al. 2022). Four hydroacoustic station triplets (H03N, H03S, H11N and H11S) with a total of 11 sensors (sensor H03N3 was malfunctioning) registered the eruption sequence. Detections begin at HA11 (at 4800 km distance to the Hunga Volcano) at 04:56:46 UTC, whereas the larger distance to HA03 (at 9350 km) resulted in a starting time of eruption signatures at around 05:46:46 UTC. Backazimuth directions from H03S were found to be 248°, while they were 149° from H11S and 150° from H11N. This direction information, the traveltime, and the observed hydroacoustic apparent velocities (1480–1500 m s⁻¹) clearly indicate a connection to the Hunga eruption. Moreover, our hydroacoustic findings are consistent with the general pattern of observed seismic waves and also correspond to the events at 04:01, 04:15, 04:18 and 08:25 UTC (see Table 2). We also identify a signal peak at 04:12 UTC, which is consistent to the REB entry at 04:12:54 UTC and the so-called T1 phase in Le Bras et al. (2022).

A back-projection from the six IMS hydroacoustic array locations towards the region nearby the eruption site indicates whether the direct propagation paths are free or blocked. The propagation is generally blocked toward HA01 (Cape Leeuwin, Western Australia), HA04 (Crozet Islands), HA08 (Chagos Archipelago) and HA10 (Ascension Island) due to the nearby bathymetry and more remote continental shelves of Australia, New Zealand and South America. This inhibits the detection of the volcanic eruption at these four stations. Also, despite our observations, the direct source-to-receiver propagation paths towards HA03 (Juan Fernandez Island) and HA11 (Wake Island) are confined by the island chain of Tonga and the Marshall islands, respectively. Therefore, we suspect that seismic waves of the volcanic eruption are converted into hydroacoustic waves at a near-source location, where a sufficient water depth for direct coupling into the SOFAR channel and a free direct propagation path towards the respective hydroacoustic station exists. Bathymetry as well as the precise backazimuth and traveltimes indicate that the region 200 km southeast of the Hunga volcano, where the Tonga-Kermadec-arc transitions to water depths larger than 1 km, is the most likely conversion region for the eruption signals observed at HA03. A similar seismic-to-hydroacoustic conversion region for propagation towards HA11 is situated around 600 km northeast of the Hunga volcano, which points at the island shelf of the country of Niue.

Broad-band infrasound signals of the eruption were observed over all recorded frequency octaves, including even gravity and Lamb waves below 0.003 Hz (Matoza et al. 2022; Vergoz et al. 2022). All 53 stations of the global IMS network detected the pressure waves of the Hunga eruption (Vergoz et al. 2022), which is an unprecedented occurrence as it exceeds by far the second-largest number of 20 stations (out of 42) that detected the Chelyabinsk fireball in 2013 (Le Pichon et al. 2013). Some detections of the Hunga main eruption phase were made days after the signal circumnavigated the globe multiple times.

The origin times determined from infrasound data in Table 2 denote the mean of the back-projected source times of all considered stations. The origin time uncertainty can be reflected by the time residuals of the single-station inversions. We first determined the origin time of the later eruption (08:27 UTC, Table 2) because this transient event was well detected by at least 20 IMS infrasound stations and is clearly recognizable as a broad-band signal in PMCC results. The time residuals are below 1 min for the closest three of the considered stations and around 2 min for IS04 (at 6782 km). In contrast, PMCC detects a series of signatures from the main eruption phase that are not that clearly distinguishable from each other due to seemingly complex source processes and propagation effects. Based on the arrivals with the maximum amplitudes, the 04:29 UTC event is formed, which differs from the independently estimated Lamb wave origin time by only a couple of seconds (our infrasound-based times in Table 2 are rounded to the minute). The residuals for the 04:29 UTC infrasound event range from 1 min (IS22 and IS04) to 4 min (IS07). For the Lamb wave, 28 IMS infrasound stations (those within 90° distance, that is within the range of the seismic stations depicted in Fig. 1 including IS23) are considered due to clear arrival onsets. The time residuals vary with a median of 6 min. When taking additional stations at westerly directions from Hunga into account for estimating the origin times of the 04:29 and 08:27 UTC infrasound events, the results differ by ±2 min for these events, where the median of the absolute residuals is also 5–6 min.

Two earlier events are found by inspection of both the PMCC detections and the waveforms (10–100 s bandpass filter). The 04:15 UTC event corresponds to the arrival times of the first broad-band signals (i.e. the PMCC configuration range). The time residuals at the four considered stations are within ±2 min. The first infrasound event at 04:01 UTC represents the first arrivals at the stations. These are higher frequency arrivals (>0.1 Hz), which could either indicate small pre-eruptions or fast arrivals of the 04:15–04:30 UTC sequence. According to Table 2, this infrasound event matches the first seismic event in time. We even find earlier high-frequency arrivals that indicate an event around 03:45 UTC, but these seem not to be correlated well at the stations (residuals between 3 and 8 min). Therefore, we do not list it in Table 2.

Regarding the accuracy of origin time of explosive events there is a difference based on which waveform technology is used. From seismological data the time can be estimated with an accuracy of a few seconds. The accuracy depends on how many detecting stations are available and how well they are distributed with distance and azimuth around the source. The IMS was expected at the time of its design to detect reliably nuclear test explosions in any medium with an explosive yield above 1 kt TNT equivalent although this threshold was not officially stated in the treaty for political reasons (Dahlmann et al. 2011). This yield corresponds roughly to a seismic body-wave magnitude of about 4.0 according to standard relations, the actual sensitivity of the IMS is much better in most regions (Gaebler & Ceranna 2021). However, it cannot be excluded that test explosions of about 1 kt TNT or slightly more could result in events with smaller magnitudes depending on the surrounding and in particular if the conducting party attempts to evade detection by decoupling measures (Council 2012). Nevertheless, we are confident that CTBT relevant seismic events are well captured by the IMS in its current shape. In addition, the propagation of the signal through the Earth can be modelled very precisely because of its structural radial symmetry. Both temporal and spatial localization based on hydroacoustic data are more challenging, because the source to receiver propagation is a combination of seismic and hydroacoustic waves with a coupling region in between (the latter is not always on the direct path between source and receiver, as stated above for HA11). Therefore the traveltime and observed backazimuth direction are not necessarily that of a direct (back-)propagation path towards the source with hydroacoustic celerities.

At almost all 53 IMS infrasound stations, signatures from the backazimuth corresponding to Hunga were detected that can be associated with the 04:15–04:30 UTC sequence. However, infrasound detections strongly depend on atmospheric propagation conditions between the source and the receivers. In contrast to the solid Earth, the atmosphere is a medium that strongly varies in four dimensions (latitude, longitude, altitude, time); the traveltime, signal amplitude and general detectability at a station depend on daytime, season, location and station azimuth. Hence, some stations only detected later arrivals of the main event (e.g. IS01 after one circumnavigation). Additionally, even within a stratospheric waveguide, which is generally known for long-range ducting and remote detection of infrasonic waves, such waves can travel different paths, potentially resulting in multiple detections of the same event. Since PMCC detects a series of coherent broad-band signals with similar wave parameters around the arrival times of signatures from the main eruption phase, distinguishing between different phases or processes such as the seismic 04:15 UTC and 04:18 UTC events is challenging. Due to different travel paths and dispersion effects over long ranges through the atmosphere, arrivals of subsequent events might overlap. Therefore, we aimed to focus on onset times of clearly separated arrivals and the detections with maximum amplitude of the dominant wave train within the frequency range of the PMCC configuration (0.02–5 Hz) to determine the origin times. Uncertainty is also introduced by the atmospheric wind and temperature profiles taken at a specific time step, which were incorporated during the ray tracing for estimating the celerities. Although the waveguide conditions along a propagation path are not expected to vary significantly within short time, the computed celerities represent only a model snapshot of the atmospheric state. For some stations located at easterly directions from Hunga, PMCC detections indicate earlier arrival times than expected for upwind propagation conditions. We therefore limited the location procedure on stations located to the west from the source, corresponding to reasonable celerities between 290 and 310 m s⁻¹ for a stratospheric waveguide.

In contrast, the Lamb wave can be easily identified within the restituted waveforms as a long-period wave with large amplitude. Therefore, we used all stations within 10 000 km (90°) for estimating the origin time. Since the exact onset times are not as clearly recognizable as for seismic P waves, for instance, we chose to consistently pick the maximum pressure peak and subtract a quarter of the dominant period as a robust estimation of the onset times, which corresponds to the starting pressure increase. Remarkably, the time residual of the closest station IS22 is the largest (17 min) following the computed celerities, indicating that the chosen method is more robust at remote stations where the absolute residuals are mostly below 10 min. When we pick the onset of the Lamb wave at IS22 directly (i.e. irrespective of the maximum amplitude and period as described above), the onset time is 7–10 min later and thus better fitting the origin time estimate (04:29 UTC). Overall, our estimation of the Lamb wave origin time agrees well with the one published by Wright et al. (2022) (04:28 UTC ±2 min).

3.2 Finding indications for explosive mechanisms in the waveforms

We examine the waveform signals to obtain indications for whether the source is of tectonic or explosive nature and start with the seismic data. Fig. 2 shows the seismic waveforms for stations up to 10 000 km (90°) with theoretical arrival times of body and Rayleigh waves based on the AK-135 structural model. The figure shows the superposition and mingling of phases which complicates the identification of single events. Almost no body-wave energy is visible. The first clearly visible onset corresponds to the theoretical arrival time of Rayleigh waves, at about 04:05 UTC for the closest IMS station at 756 km distance (AS031 or II.MSVF, Monasavu, Fiji). Possible body-wave arrivals of events at 04:15 and 04:18 UTC are obscured by the dominant Rayleigh-wave trains of the earlier events. Only at distances larger than 5000 km (45°) the theoretical onsets of body-waves from the 04:15 UTC event are earlier than the Rayleigh waves of the 04:01 UTC event. Still, the body-wave phases are hard to see as they contain very low energy, especially the S waves.

High P- to S-wave energy ratios are an important indicator for characterizing explosion-like events, which are potentially CTBT relevant. Generally, low P- to S-wave energy ratios are generated by shear rupture components and thus identify a tectonic origin for an earthquake. Theoretically, explosive sources do not radiate shear energy although small amounts exist in real data. These can be caused by reflections on structural boundaries (scattering) along with phase conversions or are due to small amounts of shear components accompanying the explosion. However, we cannot distinguish between an explosive source or a volcano eruption from the amount of S-wave energy content within the waveforms alone. Inversion of the waveforms for the source mechanism is needed (see Section 3.5).

The Rayleigh waves show in general two distinct wave trains, especially in the closer distance range up to 5000 km (∼45°). They can be connected with the events at 04:15 and 04:18 UTC, respectively. The generally dominant Rayleigh waves, next to the low S-wave energy in comparison to the P-wave energy, are a first indicator for an explosive mechanism.

The frequency range of the eruption signal is very low, oscillating between 0.02 Hz and roughly 0.15 Hz, depending on the distance of the station to the volcano (Fig. 3). However, the signal could be best presented in the waveforms if the bandpass filter is further constrained to a maximum of 0.06 Hz, thus clearly focusing on surface wave energy. Above roughly 0.1 Hz the signal of the secondary oceanic microseism (0.06–0.08 Hz and 0.11–0.17 Hz for 1st and 2nd microseism, respectively, Friedrich et al. 1998), obscures the signal of the volcanic eruption.

In contrast to the waveforms, the spectrograms in Fig. 3 show the onset of the eruption signal with the theoretical P-wave arrival of the event at 04:01 UTC only at the closest station Monasavu on Fiji (IMS station AS031 or II.MSVF, 754 km). At the two later stations, only the Rayleigh wave of the 04:01 UTC event produces significant amount of energy. This energy is in the same frequency range as that of the most prominent and obvious events around 04:15 UTC, thus suggesting a common nature and indicating that the eruption sequence started already with the 04:01 UTC event. The energy of the Rayleigh waves concentrates around 0.05 Hz. Around this frequency, we see varying energy strength over time which is the on- and off-rising of the Rayleigh wave energy of the different events overlapping each other.

While for explosions the duration plays a minor role, a detailed analysis of the source time function may provide information for discriminating between anthropogenic and natural events, such as a volcanic eruption. By considering the waveform shape and its spectral content in time and frequency domain, respectively, it has been demonstrated for exemplary comparisons between explosive and natural sources, that for all three CTBT waveform technologies such differences are indicated (Schwardt et al. 2022). From analysis of mainly infrasound along with hydroacoustic and seismological data Vergoz et al. (2022) estimated a total duration of the main eruptive phase of about 1.5 hr from 04:15 to 05:45 UTC with 15 min of increasing activity starting at 04:00 UTC. Following this time of strong activity, Vergoz et al. (2022) determined an additional time span of lower activity until 08:15 UTC and state that only after the event at about 08:30 UTC the activity decreased below the level of detection. We estimate the duration of the main eruptive phase from the seismic spectrograms, averaged over the frequency range of the eruption signal (0.02–0.09 Hz). Fig. 5 shows the spectral power (SP) for stations up to ∼4100 km distance aligned by the Rayleigh wave velocity of 3680 m s⁻¹ and compensated for geometrical spreading, as well as their stack. As expected the amplitudes distinctly increase with the first event at 04:01 UTC and reach their maximum around 04:15 UTC. Subsequently, they drop again but maintain a constant, relatively high level until 05:30 UTC before they start to decrease again. At about 06:15 UTC the decrease even accelerates. Therefore, we also estimate a length of the main eruptive phase of 1–2 hr. In addition, we see that the stacked SP-curve decreases to an amplitude minimum around 07:18 UTC before it reaches a second maximum at 08:25 UTC and then quickly decreases again to the level before the eruption. Overall, the long duration of the event as observed in all three waveform technologies excludes a classification as an anthropogenic explosion-like event.

An explosive mechanism in infrasound observations is indicated when the waveform of the direct arrival shows an impulsive N-shape (e.g. Fig. 2, middle panel, Vergoz et al. 2022). This distinct, high amplitude pressure pulse of short duration with broad-band frequency content and a follow-up pressure settling is indicative of a blast wave, well described by a Friedlander waveform (Reed 1977), also applicable to volcanic eruptions (Marchetti et al. 2013). A non-explosive waveform shape like, for example an earthquake normally lacks the impulsive, high frequency components of this effect and its waveform is smoother, often lens-shaped over a longer duration. The infrasound observations of the Hunga Volcano eruption, especially the largest-amplitude Lamb-wave signal components, clearly show the N-shape characteristic, representative of a strong blast wave.

Furthermore, the global-scale detection of a solitary Lamb wave in terms of hPa pressure variations is a further indication that an extremely strong point-like—likely explosive—source and not an earthquake was recorded at infrasound arrays. Only the largest anthropogenic sources like nuclear explosions in the Mt TNT equivalent range or similar events like large volcanic eruptions have the necessary eruptive potential to generate such waves (Matoza et al. 2022). This typically excludes earthquakes, which do not directly produce an atmospheric blast wave of such intensity as an above ground infrasound source.

Within hydroacoustic monitoring, another effect can be investigated that is indicative of an underwater explosion. This is the spectral scalloping caused by bubble pulse and water column reverberation (see also Section 1). It is typically identified within the frequency spectrum as distinct peaks (Koper et al. 2001; Vergoz et al. 2021). For both the hydroacoustic waveform records at the HA03 arrays and the southern array at HA11 we cannot identify clear indications of spectral scalloping, for neither the 04:15 nor the 04:18 UTC events. In addition, no distinct proof of this effect is provided in corresponding hydroacoustic data and spectral analyses, as shown in Le Bras et al. (2022). Hence, the lack of evidence of spectral scalloping might indicate that the explosion associated with the eruption did not occur freely in the water column, but either (partly) above the water or (partly) inside the volcano structure, as supported by bathymetry measurements indicating a 150 m deep water column above the caldera (NASA Jet Propulsion Laboratory 2022c). It must be considered however that the Hunga eruption is an extended source in space and time compared to ideally point-source explosions which are better known for scalloping effects. The stronger amplitudes of the 04:15 UTC event in comparison to the 04:18 UTC event in the hydroacoustic data suggests a better coupling of the first event to the SOFAR channel.

3.3 Determining magnitude

Determining the magnitude and thus the strength of an event under investigation is key in characterizing its CTBT relevance, because magnitude usually is the base to estimate how much kilotons of TNT equivalent were detonated (see subsection on yield estimation). The main issue lies in the definition of the various magnitude scales. In general, they are derived from hard rock on-land environments and are developed and scaled for tectonic earthquakes. Explosive sources differ significantly in their characteristics from earthquakes which influences the magnitude estimation. Most importantly, exothermic contributions during an explosion (among other unconsidered effects) need to be taken into account. Nevertheless, the common magnitude scales for earthquakes, especially the body-wave magnitude m_b, are used in the frame of CTBT due to a lack of alternatives. For the six North Korean nuclear tests since 2006, magnitudes between m_b 3.8 and m_b 6.4 were estimated with reasonable accuracy (Gaebler et al. 2019).

The Hunga eruption can be appraised as a complex process, as it is (i) neither a tectonic earthquake nor a classic explosive source and (ii) neither on- nor offshore but occurred in shallow waters more in between. This intermediate location leads to eruption energy being partly radiated laterally into the atmosphere and also into the water column, coupling into the SOFAR channel, and reducing the amount of energy coupled into the ground. Thus, the above mentioned difficulty is even increased and all magnitude estimates are more or less underestimated. For the 04:15 UTC event, the IDC determined a magnitude estimate of m_b 4.2 based on data of 11 stations at distances of 2900 km (26°) to 11 000 km (96°). In contrast, the IRIS and USGS catalogues provide estimates of M_S 5.8 and m_b 5.3, respectively, whereas the EMSC catalogue provides a moment magnitude estimate of M_w 5.8.

We did not determine a body-wave magnitude because the P-wave signals on the global seismological data are hard to analyse due to a low signal-to-noise ratio. As an alternative, we used the data of GRSN stations to estimate a magnitude based on surface waves M_S and estimated a value which varies between 5.8 and 6.3 with a median of 6.0, wrapping up all surface waves generated by the event at 04:15 UTC. In addition, based on waveform inversion for moment tensors we determined moment magnitude estimates of 5.0, 6.2 and 6.0 for the 04:01, 04:15 and 04:18 UTC events, respectively. For the 04:15 UTC event, our estimate is larger than the body wave magnitude estimate from EMSC. Still, we assume that the seismic magnitude estimate potentially underrates the full energy released into the atmosphere during these complex events, as indicated by the impressive effects this eruption produced on a global scale (e.g. Lamb waves circumnavigating the globe four times). This might be because the forces of the eruptive processes partly do not couple well to the solid Earth.

In summary, the eruption demonstrates that we need magnitude relations (i) for non-earthquake events based on seismic data and (ii) calibrated specifically for infrasound and hydroacoustic data.

3.4 Estimating the explosive yield

The complex nature of volcanic processes and the long duration of the energy release during an eruption are not comparable to the energy released during a nuclear explosion, which is for all intents and purpose a point source. The long duration and complex source processes of the Hunga event (see Figs 6 and 7) are clear indicators of the volcanic eruption nature of the event. A volcanic explosive eruption is also not an explosion, but a mixture of source processes, possibly dominated by a single force (Kanamori et al. 1984). For this reason the assumptions and relations typically used to estimate the explosive yield from the data of all three waveform technologies do not apply. We do however offer in the following some yield estimates based on some assumptions and relations typically applied in the CTBT context to highlight the discrepancies in yield estimates between the technologies and the specific difficulties arising for each waveform technology.

Typically, yield estimates are obtained based on relations to body-wave magnitude m_b (Murphy 1981; Bowers et al. 2001; Bormann et al. 2009). In general, it is not possible to state one single relation between magnitude and yield because it is influenced by various effects, such as geological setting and wave propagation. In the case of the Hunga volcanic eruption, the difficulties in determining the magnitude in the first place (see section above) make it challenging to provide a yield estimate based on seismological data. The broad variety of m_b estimates from different sources (see Table 1) would estimate the yield to below 1 kt and up to 50 kt TNT. In our case, the M_w 6.2 for the event at 04:15 UTC relates to an m_b 6.0, which would result in an estimate of roughly 300 kt TNT (Das et al. 2013; Bowers et al. 2001). All these values underestimate the fierce energy of the volcanic explosion drastically.

In contrast, Vergoz et al. (2022) estimated the yield to 100–200 Mt of TNT from infrasound data using Lamb wave amplitudes. While empirical relations for nuclear explosions like the AFTAC period-yield relationship (ReVelle 1997) are not applicable to extended eruption sources in general and the frequency regime of the Hunga eruption in particular (see supplement section 2.8 of Matoza et al. (2022) for a discussion on this), an approach modified from Pierce & Posey (1971) focusing on the eruption’s Lamb wave is applied by Vergoz et al. (2022) to provide an approximate yield estimate. Díaz & Rigby (2022) used an empirical pressure–distance relation calibrated by data from the 1980 Mount St Helens eruption (∼35 Mt TNT) and the 1961 Tsar bomb (∼57 Mt TNT) to obtain an estimate for the energetic output of the Hunga eruption. This results in 61 Mt TNT as a lower boundary value. However, from their Fig. 5, also values around 100 Mt TNT would be explained by the measured pressure amplitude. In comparison to the Tsar bomb, the Mount St Helens eruption, or the 1883 Krakatau eruption (∼100–150 Mt TNT), these values seem to be reasonable estimates for the Hunga eruption. As a reverse conclusion, assuming a yield of 100 Mt TNT and applying the relation of Bowers et al. (2001), a corresponding magnitude equivalent would be m_b 8.0. A similar discrepancy due to non-ideal coupling into the ground was already found in the investigation of the 1961 Tsar (aerial) bomb whose magnitude was estimated to m_b 5.0 (corresponding to ∼10 kt TNT) whereas the yield estimate was 50–58 Mt (Khalturin et al. 2005).

Underwater explosions generate hot gases that quickly rise. The gas bubble oscillates with a dominant period, which should be clearly visible in hydroacoustic and, to some extent, also in seismic data. Based on relations between bubble period, depth of the explosion, and yield strength (e.g. Chapman 1985), it would be possible to estimate a yield value. However, we do not see indications for such a signal, neither in the hydroacoustic nor the seismic data.

3.5 Discrimination based on source mechanism

For further verification of the source process, we perform a seismological waveform inversion. This is an important step in classifying potential events of CTBT relevance, as it provides information on a full moment tensor point source. Considering a volcanic source process, targeting, for example a volume point source or a combination of moment tensor and single force is probably a more appropriate assumption on the underlying source process (Kanamori et al. 1984; Cesca et al. 2020). However, in the literature a moment tensor source is often used to analyse volcanic processes as well (Shuler et al. 2013; Gudmundsson et al. 2016). Here, in addition to the moment tensor, we also invert for a vertical single force source. We carry out five different moment tensor inversions. One with FDSN data only for the 04:01 UTC event and two for each of the strong events of interest at 04:15 UTC and 04:18 UTC, with one inversion for each event with only FDSN data used and one where only IMS data are used. This approach allows for the assessment of the necessity of additional non-IMS data.

The inversion results, which are targeting a moment tensor for the 04:01 UTC event, indicate an origin time of 04:01:30 UTC and a moment magnitude of M_w 5.0 for the best fitting model. For the entire parameters of the solution ensemble and further result details see the uploaded inversion reports (Donner et al. 2022). The visualization of the solution within the Hudson space in Fig. 8 a illustrates the explosive character of the 04:01 UTC event. The Hudson plot also shows the non-negligible component of crack opening of the source. The spread of solutions within the Hudson plot shows the ensemble of the determined full mechanisms from bootstrapping, providing an estimate of the robustness of the solution. The full moment tensor decomposition in Fig. 8(b) clearly shows again the explosive character of the source (positive isotropic part, or abbreviated as ISO) with only minor double-couple (DC, shear component) but a significant positive subvertical compensated linear vector dipole (CLVD) part.

Fig. 9 visualizes the distribution of the parameter estimates from the bootstrap ensemble for each of the six moment tensor components smoothed by a Gaussian kernel. All components are well resolved for the event at 04:01 UTC in Fig. 9(a) indicating a stable moment tensor solution. Especially the components M_nn, M_ee and M_dd, which are the trace (isotropic) components of the moment tensor providing information on the volume change of the source, are very well confined. Interestingly, the component M_dd is about four times the value of M_nn and M_ee. According to Kumagai (2009) this is a typical pattern of a vertical tensile crack assuming a relation between the elastic moduli of λ = 2μ. Shuler et al. (2013) discussed that such mechanisms can be best explained by ring faulting (a ring around the crater walls); in case of strike angles larger than 180° and positive rake angle it is oriented as outward dipping reverse.

The components M_nd and M_ed are less well resolved due to the fact that the source is close to the (seafloor) surface. It is commonly known that for shallow sources the inversion for the moment tensor is ill-conditioned due to the free surface condition (Dufumier & Rivera 1997; Julian et al. 1998). However, the components M_ne and M_nd do not contribute much to the solution, as visible from the high PDF values for the parameter value at zero. In contrast, the component M_ed does contribute with a positive value but with larger uncertainty. It may be interpreted as deviation of the volcano geometry from the vertical axes. Also, it can likely be explained by interaction with the local stress field or by unmodelled complexities in the source process.

Fig. 10(a) shows examples of P-phase waveform fits between observed and synthetic data for the event at 04:01 UTC. The complete P- and S-waveform fits can be seen in the Grond report uploaded to a repository (Donner et al. 2022). The fit is in general quite good and shows a positive first motion on all stations. This is another clear indicator for a positive isotropic source mechanism, that is an explosion. Interestingly, the second positive pulse of the P wave is underestimated in amplitude at some stations. This might be a hint that the underlying source process favours a single force or a set of single forces instead of a full seismic moment tensor as a source. Not shown here but available in the repository Donner et al. (2022) are the fits between observed and synthetic S-wave waveforms. They are not well fitted for all three events. This is in agreement with the events having a mostly isotropic character and therefore not exciting strong S-wave energy.

The inversion results for the event at 04:15 UTC indicate a moment magnitude estimate of M_w 6.2 for the best fitting model. Overall, the general results are very similar to the ones from the event at 04:01 UTC. We see a clear explosive source with significant positive vertical CLVD component, representing a tensile crack opening. The solution is much more stable (Fig. 9b) than for the event at 04:01 UTC, visible from the strong focus in the Hudson plot in Fig. 8(b) for FDSN stations and in Fig. 8(d) for IMS only stations. The effect of underestimating the second positive pulse of the P wave in Fig. 10 is even stronger than for the event at 04:01 UTC. Both the FDSN and IMS data set based moment tensor inversions resolve the same source processes, with the FDSN data set performing slightly better in terms of resolution (Figs 8b and c).

For the event at 04:18 UTC we determine a moment magnitude estimate of M_w 6.0 for the best fitting model (Fig. 9c). Both moment tensor inversions results of the FDSN data set (Fig. 8c) and IMS data set (Fig. 8e) are comparable, however again the FDSN data set performs better in terms of resolution. In contrast to the two events before, we infer an event with implosive character (negative ISO part) and a negative vertical CLVD, indicating a combination of downward contraction and horizontal expansion. This is a common observation of volcanic calderas during collapse, implying failure of support structures, either directly above or within the magmatic reservoir (Gudmundsson et al. 2016). According to Shuler et al. (2013) such mechanism with strike angle 180° and larger and negative rake angle can be interpreted as outward-dipping normal ring fault.

This switch in polarity sign of the causative process can be observed even visually in the waveforms in Fig. 10. For the event at 04:15 UTC (Fig. 10b) we see a positive first motion, while for the event at 04:18 UTC (Fig. 10c) it is negative. Regardless of the sign switch, again the amplitude of the P-pulse is underestimated by the synthetics, but in contrast to the other two events it is the first instead of the second P-pulse which is underestimated.

Due to the indications that the source processes of all three analysed events may favour a single force over a moment tensor, we repeated the inversion. This time we target a vertical single force mechanism which might be a better model of a volcanic eruption (Kanamori et al. 1984). For the event at 04:01 UTC the results are highly unstable. The best fit for the vertical single force amplitude for events at 04:15 UTC and 04:18 UTC yields 4.68 × 10¹⁵ N and −6.89 × 10¹⁵ N, respectively. Thus, the force of the event at 04:15 UTC is directed upward, while the force of the 04:18 UTC event is directed downward. This pattern is congruent with the polarity switch between events during waveform inversion for moment tensors. However, the upward force contradicts the classic eruption model of Kanamori et al. (1984). In addition, the underestimation of P-wave amplitudes between observed and synthetic waveforms is evident once more. Again, for the event at 04:15 UTC the second (positive) P-pulse, while for the event at 04:18 UTC the first (negative) P-pulse is underestimated in amplitude by the synthetic waveforms. We assume that in an inversion setting where the waveform inversion targets both single forces and subsequent events jointly, these second order underestimations would be solved. However, within this paper we concentrate on the application of a routine processing as would be applied within a CTBT context and only in the second instance on details of the eruption process. Within this framework the current single source inversion method proves sufficiently capable at first order for the use in the CTBT context to discriminate source types and invert for the acting forces.

The seismic waveform inversion for the source mechanisms clearly shows that a tectonic earthquake can be excluded as possible source. With respect to CTBT that means the suspicion against the investigated event has been increased. Though, at this stage it is unclear whether an explosive source is of nuclear character or not. Nevertheless, the Hunga eruption has also underlined that more sophisticated tools are required for characterizing the real nature of an explosive source. Based on the waveform inversion alone, we would have difficulties to recognize the volcanic eruption, although the high CLVD part is a hint in that direction. The location of the source in (shallow) water is complicating things even further.

3.6 Source time function

In order to conduct a clandestine nuclear test, it is theoretically possible to add slightly time-shifted explosions in the immediate neighbourhood to the actual test. As a result, all three waveform data sets would show a complex pattern of signals overlaying and obscuring each other. A possible approach to identify such a case would be the analysis of the source time function.

Theoretically, a pure explosion can be sufficiently described with a single triangular function. From the waveform inversion, we inverted source time functions based on the moment tensor of the 04:15 UTC event, an explosion, and a single force (Fig. 6). The narrow shaded area of uncertainty, obtained by repeating the inversions with different station weighting, demonstrates the stability of the resulting source time functions. The origin times of events at 04:01, 04:15 and 04:18 UTC (vertical green lines) correlate with rather the source time function peaks instead of their starting point.

Poli & Shapiro (2022) determined a source time function by a simple and straightforward approach in the frequency domain. Here, the amplitude source spectra can be determined by dividing the Green’s function amplitude spectra from the seismogram amplitude spectra. The final force spectrum is then the average of the thus calculated single station spectra and is constrained to positive values. In contrast to Poli & Shapiro (2022), we invert directly in the time domain using a least-squares approach. The resulting source time function can therefore be positive or negative. Hence, for example, in the case of the vertical single force, it can act vertically from above or below the (seafloor) surface on a point. Furthermore, instead of surface waves, we focus on body waves and thus avoid contamination due to scattering or conversions. Nevertheless, the peaks in our source time functions correlate well with peaks in the source time function from Poli & Shapiro (2022) shown in Vergoz et al. (2022).

While the event at 04:01 UTC is weakly pronounced, the events at 04:15 and 04:18 UTC show very distinct peaks. Again, the peaks switch from positive to negative sign between the events as already seen in the waveform inversion for moment tensors and single forces.

In general, the three inverted source time functions based on different source types are remarkably similar in their basic features. This similarity might be an effect of the chosen time window shifted over the data and concentrating on the P-wave. Thus, the inversion might be stabilized by suppressing other effects, such as surface wave reflections from atmospheric interaction. Nevertheless, the similarity demonstrates that the basic physical processes of the eruption are dominant enough to be independent of the assumed source mechanism (i.e. single force, moment tensor, or explosion). There are, however, a number of small and interesting differences. Directly after the event at 04:15 UTC and before the event at 04:18 UTC, the single force based source time function shows a strong positive peak of the source time function, which is less pronounced in the other two source time functions. Moreover, in the first 3 min after the event at 04:18 UTC there are two weaker negative peaks visible in the source time function for the moment tensor and explosive sources, but only a single one for the single force source.

At 04:27 UTC, another peak is visible in all three source time functions, aligning maybe with events detected from infrasound data (04:29 UTC in Table 2). An even more distinct signal is visible at 05:30 UTC corresponding to a seismic event in the IRIS and the REB catalogues. At about 06:15 UTC, the signal amplitude falls to the background level, which is about 20 min later than from the source time function of Vergoz et al. (2022). However, caution is needed in general when interpreting the tail of a source time function, because surface wave energy is dominant in the low frequency band used by both, us and Vergoz et al. (2022). Therefore, an overestimation of the duration of the actual source processes might be possible.

The array beam on the PKP-phase determined from GRSN data at 16 700 km (∼150°) distance (Fig. 11) also confirms the patterns found in the source time functions. The PKP-phase needs roughly 20 min to travel this distance through the mantle and the outer core. The first PKP-onset arrives at 04:35 UTC corresponding to the event at 04:15 UTC. After a gap of about 220 s, a second wave train arrives at 04:39 UTC corresponding to the event at 04:18 UTC. The second wave train shows at least two distinct peaks with first positive, then negative polarity. The second peak is not as distinct as the one before. It could be assumed that more than one PKP-phase arrives here, mixing with each other. The first peak in this later wave train is, interestingly, not constant in amplitude over all GRSN stations, in contrast to the other peaks. The first intuition might be to assume a radiation pattern. However, the backazimuths of all GRSN stations with respect to the Hunga Volcano vary only between 5° and 17°.

The more straightforward approach of stacking allows us to determine source time functions not only from seismic but also from hydroacoustic and infrasound data (Fig. 7). Hydroacoustic data need a positivity constraint because the polarity information of the phase is lost due to the unknown number of reflections during travel in the SOFAR channel. This is similar for the infrasound data, where the information on phase is difficult to obtain due to the complex structure of the atmosphere and the propagation of the signal therein. For consistency reasons, we applied the positivity constraint to all three waveform data during stacking.

The seismic waveform data stacks of the FDSN and IMS stations are remarkably similar (Fig. 7b), highlighting the good quality of the IMS stations. For the investigation of the Hunga eruption processes the IMS station coverage is therefore enough to produce analysis results of similar quality. The events at 04:15 and 04:18 UTC are much more pronounced compared to the event at 04:01 UTC. Surprisingly, the event at 04:01 UTC is very clearly visible in the stack of the hydroacoustic data, suggesting a better coupling to the water column than to the solid ground. The peak for the event at 04:15 UTC has an impulsive increase and decrease in the seismic stack. In contrast, the peak for the event at 04:18 UTC starts impulsive but is followed by a lot of energy reverberations.

The peaks for the respective events can be found in the stacks of the seismic data and hydroacoustic data and as well in the infrasound stacks. However, the prevalence of the peaks in the infrasound stacks are dependent on the subset of stations used for stacking (Fig. 4). In all three infrasound stacks done for different subsets of arrays the onset at 04:01 UTC is visible and the maximum of the infrasound energy release is at around 04:30 UTC, corresponding with the localized events (Table 2). In the easterly propagation direction stack (Fig. 4a) the 04:18 UTC event is roughly indicated. The stack for the larger subset of arrays in the westerly propagation direction (Fig. 4b) clearly indicates the 04:18 UTC as an impulsive source. Only the stack for the smaller subset of arrays in the westerly propagation direction (Fig. 4c) also clearly indicates the 04:15 UTC event. Despite the comparable propagation conditions towards these two westerly propagation direction stacks, propagation and dispersion effects increase with range and likely decrease the quality of the stack. We therefore choose the stack of the smaller subset of arrays in the westerly propagation direction for the comparison with the stacks of the other three waveform technologies in Fig. 7, but advise caution that the 04:15 UTC event from infrasound data alone cannot be established without doubt.

Yuen et al. (2022) also estimated a source time function from stacking seismic data between 7800 km (70°) and 10 000 km (90°). Similarly, they see a clear pulse for the event at 04:15 UTC (E1 in their fig. 6) but a more complex pattern for the event at 04:18 UTC (E2 in their fig. 6). The latter shows three pulses within 100 s. They correspond to several higher peaks in our seismic stack and energy reverberations in the infrasound and hydroacoustic stack. Also, Thurin et al. (2022) inverted for a source time function with positivity constraint and found the same four subevents as Yuen et al. (2022) in similar time ranges.

From the hydroacoustic stack, the eruptive activity seemingly comes to rest at 04:40 UTC. However, sporadic peaks of activity (e.g. at 04:45 and 05:30 UTC) are evident until 06:00 UTC. Likewise, the elevated activity on the infrasound and seismic stack lasts until 05:50 UTC and 06:05 UTC, respectively. That is a duration of roughly 1.5–2.0 hr, similar to the estimate of Vergoz et al. (2022).

In summary, the interpretation of the source time functions is complex for both, from inverted and from stacked waveforms, including a lot of detail on the eruption process. For investigating the potential presence of a nuclear explosion non-compliant with the CTBT the more relevant information is limited to not being an explosive source because it is too complex. Otherwise, following radionuclide detections combined with atmospheric transport modelling for radionuclides would be the only sure discriminator for a forbidden nuclear test.

3.7 Detection of radionuclides and modelling atmospheric transport

The final step in the CTBT verification chain is to search for traces of radioactive isotopes in the atmosphere. In the case of the Hunga Volcano obviously no CTBT relevant short lived fission products were expected to be detected. Indeed, IMS radionuclide stations detected no significant anomaly in activity concentrations of radioisotopes which we could attribute to the volcanic plume, neither natural nor anthropogenic. Fast uplift of the plume generally reduces the chance of ground based radionuclide measurements. However, in the fictive case of a nuclear detonation near the Earth’s surface during the Hunga eruption, a large part of the radioactive debris would have been dispersed. Even if only a small portion of effluents would have remained in or returning to atmospheric layers close to the Earth’s surface, they would have been detectable due to the high sensitivity of radionuclide aerosol samplers. In addition, radioactive xenon isotopes could be measured also at larger distances as they are removed from the atmosphere by decay only. Elevated sulphur-dioxide concentrations above Fiji led to warnings of acid rain in the week after the eruptions (Sun 2022). It has to be kept in mind that these traces could also originate from the pre-eruption starting on 13 January 2022. Nevertheless, taking those as plume proxy, it is justified to assume that for a nuclear explosion during the Hunga eruption, with direct venting of particles into the atmosphere, traces of radioactive fission products would likely have been detectable at RN26, Fiji. The actual dispersion and evolution of the Hunga Volcano plumes, for the 15 January eruption mainly in higher altitudes, are well tracked by satellite observations of sulphur-dioxide, water vapour and aerosols (NASA Jet Propulsion Laboratory 2022b)—with remanence in the stratosphere over months and not yet fully assessed climate effects (Taha et al. 2022).

Applying ATM is challenging for the Hunga Volcano as the explosive eruption itself heavily dominated atmospheric dynamics in the source region. Therefore the typical models are not able to adequately represent the development of the atmospheric state there. Usually only the troposphere up to about 11 km altitude is considered for the transport of radioactive traces from nuclear explosions. The Hunga volcanic plume reached about 55 km altitude and the blast caused plume wind speeds directed away from the source in the range of 50 m s⁻¹ (Carr et al. 2022). Another indicator for the extraordinary vertical dynamics in the atmosphere are the extreme lightning rates observed at the volcanic cloud column, actually the highest ever observed. Even a Terrestrial Gamma-ray Flash (TGF), an invisible highly energetic gamma-ray pulse generated by the acceleration of charged particles within lightning, was captured from satellite measurements and is considered to be the first observed TGF produced by volcanic lightning (Briggs et al. 2022). Monitoring electromagnetic pulses from space is also used as national monitoring technology for the detection of signals from nuclear explosions (Sandia National Labs 2016). Also, the future use of ground based lightning detection networks as additional technical mean for CTBT monitoring could be investigated as they are working in real time with high location accuracy.

Independent of a fictive nuclear test, it is interesting to see whether natural radionuclides from the volcanic ash could be measured. Volcanic ash is known to have an abundance of natural radionuclides, mainly products of the uranium and thorium chains including radon decay products. Atmospheric evidence of such isotopes in volcanic ash was reported, for example for Mount St Helens (Kuroda et al. 1984) as well as Pinatubo (Sato et al. 1994). For the investigation of magmatic gases the ²¹⁰Po/²¹⁰Pb ratio is routinely used and an excess has been found in local measurements at Eyjafjallajökull in 2010 (Sigmarsson et al. 2015). IMS stations are not sensitive to ²¹⁰Po. Attributing other measurements of volcanic ²¹⁰Po at distance is difficult due to largely variable background emissions of biomass burning and coal-fired power plants. Within the CTBT-IMS at the stations RN71 (Sand Point) and RN76 (Salchaket) in Alaska, untypical detections of the potassium isotope ⁴²K and the uranium child protactinium ^234mPa occurred in coincidence with the passage of ash plumes from the Aleutian volcano Bogoslof in 2016/2017 (Kuśmierczyk-Michulec et al. 2021). In contrast to these earlier studies, for the period after the Hunga Volcano eruptions, we could not detect a significant elevation of natural radioactive isotopes above the typical concentration range at adjacent IMS stations.

With respect to CTBT issues another interesting aspect is the question of the ‘readiness’ of the IMS. We assume a hypothetical scenario of a nuclear explosion at the location and time of the Hunga eruption (ignoring that it was an eruption) with leakage of CTBT relevant trace substances. For such a scenario, backward ATM revealed considerable sensitivity to the source region for the IMS station RN26 on Fiji. In this case the source areas of air masses relevant to the sample were traced without considering specific radionuclide properties as decay or deposition. For the sample taken on 18 and 19 January 2022, the sensitivity result to the source time step after the main eruption phase (15 January, 06:00–09:00 UTC, Fig. 12a) shows a sensitivity >10⁻¹⁵ m⁻³ for the ground layer (0–200 m altitude). The forward modelling for a hypothetical release at the time of the main eruption phase gives comparable relative concentrations as shown in Fig. 12(b) for 18 January 2022. As the minimum detectable concentrations of key fission products for CTBT monitoring are in the range of a few μBq m⁻³ in this scenario an activity release of several 10⁹ Bq would have been detectable. Depending on the decay and deposition of the respective radioactive aerosols the releases have to be higher. Typically assumed generated quantities for various isotopes produced by nuclear explosions of about 1 kt TNT are in the order of 10¹⁴ Bq and therefore well above the detection threshold if largely released.

Furthermore, it can be expected that traces of radioactive noble gases would have shown up at stations at greater distance. The station RN26 is not equipped with a noble gas system. However, full containment of radionuclides is possible for emplacements of the nuclear explosive device underground or under seabed. In cases of full containment it is intrinsically justified to expect a better seismic coupling which makes the characterization of the explosion by seismological methods more convenient than in the quite complex and ‘open’ Hunga event series. For nuclear test explosions underwater, the expected radionuclide releases into the atmosphere could also be minimized. However, we would expect a stronger focus on hydroacoustic recordings.

4 ERUPTION INTERPRETATION

Zimanowski & Büttner (2003) determined the detailed physical processes during phreatomagmatic explosions in subaqueous systems from laboratory experiments. They summarized the processes in four phases which we are able to clearly relate to specific signals in our data.

In phase 1 water mingles with magma under distinct hydrostatic pressure conditions keeping this pre-explosive mixture (premix) stable. At the contact stable film boiling conditions (Leidenfrost phenomenon) set in. The reason for the mingling may be ascending magma or intruding water. This process can last for seconds to minutes on length scales of millimetres to decimeters, depending on the viscosity of the magma. An important controlling factor for this phase is the hydrostatic pressure of the water column, which should be large enough to provide sufficient pressure for stabilization but shallow enough to not suppress a blast. Bathymetric investigations showed that the caldera of the Hunga Volcano was covered by a water column of about 150 m providing quite perfect hydrostatic conditions (Amos 2022). We identify the start of this phase at 04:01 UTC where smaller explosions lead to fracture opening, hence, allowing first amounts of sea water to enter the magmatic system. This activity occurred submarine with little atmospheric coupling. Hence, we see the distinct onsets in the stacked source time functions of the hydroacoustic and seismic data but only a minor increase in the stack of the infrasound data (Fig. 7). Also the source time functions of the inverted seismic waveforms show clear signal onsets for all three source mechanisms (single force, moment tensor, explosion) at that time (Fig. 6). Moreover, the moment tensor solution for the 04:01 UTC event supports our interpretation, showing an explosive mechanism including vertical crack opening, which might be the opening of a vertical vent, as discussed above (Fig. 8a). During the following ∼15 min the mingling of water and magma under stable conditions seems to continue, as visible by increased energy and partly with distinct peaks in the source time functions from inversion (Fig. 6) and stacking (Fig. 7). From the infrasound and hydroacoustic stacks (Figs 7a and c), we see that this process phase is accompanied by possibly minor explosions and significant seismic unrest, respectively.

Phase 2 is characterized by a destabilization of the water-magma mingling due to high mingling energy. A sharp pressure pulse triggers the break-down of the vapour film in the premix. The pressure pulse could be caused, for example, by brittle failure of rock walls, rapid condensation of vapour, or a shock wave. This trigger phase 2 directly leads to the fragmentation phase 3, a process by which a continuous liquid or solid phase is transformed to a continuous gas phase with suspended particles. It is caused by the direct contact between cool sea water and very hot magma leading to (repeated) detonations and converts thermal to mechanic energy. Both phases last only micro- to milliseconds before they lead to phase 4, the abrupt expansion of the magma-water-vapour-mixture, that is the phreatomagmatic explosive eruption of the volcano.

In our data, we see changes in the physical processes between 04:10 and 04:14 UTC. The inverted source time function for the single force shows stronger variability, while the ones for the moment tensor and explosion show peaks (Fig. 6). Also, the source time function stacks show signals in the same time frame, especially in the hydroacoustic stack (Figs 7c and f). We interpret that we see here the transition of phase 1 to the phases 2 and 3 with distinct pressure build-up. The REB of the IDC and the IRIS catalogue list a stronger seismic and hydroacoustic event at about 04:13 UTC (Table 1). Maybe the ground shaking due to this event was the trigger to the destabilization of the water-magma mingling phase. It could be interpreted as the cracking of the caldera rock caused by thermal stress. The large explosive event at 04:15 UTC is then connected to the explosive eruption of phase 4, visible by a clear, impulsive, positive peak at all source time functions and supported by the moment tensor solution (Figs 8b and d). Based on the source time function from inversion, we may interpret a secondary weaker explosive event at about 04:17 UTC.

Around 04:18 UTC to 04:21 UTC another sequence of events with impulsive phase onsets is visible in all three waveform technologies as at least one distinct peak in the infrasound and hydroacoustic stacks and as several peaks in the seismological stack (Fig. 7). The first event of this seismic event sequence at 04:18 UTC is characterized mostly by forces in the opposite direction than the 04:15 UTC event, as visible in the moment tensor solution (Fig. 8) and all source time functions (Figs 6 and 7). Interestingly, the two opposed pulses have approximately the same strength, indicating that about the same volume amount and/or explosive strength is involved in both processes. The following events of the 04:18 UTC to 04:21 UTC subsequence are however indicated by the source time function inversion to be again of a more explosive characteristic (Fig. 6). The seismic moment tensor analysis is not reliably possible for the individual later events of this 04:18 UTC to 04:21 UTC subsequence, due to the coda. However our tests with longer duration of fitting windows of the moment tensor inversion for the first event at 04:18 UTC indicate that both implosive and explosive characteristics can be inferred. This implies that the extended 04:18 UTC to 04:21 UTC subsequence is made of complex source processes. Bathymetric measurements showed that the caldera dropped in depth to more than 850 m below sea level (Amos 2022; O’Callaghan 2022; Stern et al. 2022), indicating that implosive mechanisms are plausible.

In the following 20–30 min we see continued activity in all source time functions with partly more or less clear peaks. We interpret this time frame from the seismological and hydroacoustic observations and analysis as further (minor) explosion activity and/or partial collapses as sea water continues to enter the cone and gets in contact with hot magma. The interpretation is strengthened by testimonies of eye witnesses who reported ‘crackling and noise like artillery fire’ (O’Callaghan 2022). This is in contrast to the (stacked) infrasound activity peaking around 04:30 UTC. It is in good agreement with our origin time localization of the maximum-amplitude infrasound event and the estimated Lamb wave origin time, both around 04:29 UTC (see Table 2). From the stacked source time functions, the main activity phase comes to a temporary end at about 05:50 UTC. At the inverted source time functions the signal decreases to a background level at about 06:10 UTC. That makes a total eruption duration of about 2 hr in agreement with the 1.5–2 hr of Vergoz et al. (2022).

5 CONCLUSION

With this paper we made use of the natural test case of the large Hunga volcanic eruption in January 2022 to evaluate established routines and methods in the context of the CTBT’s verification regime. We used all three waveform technologies (SHI) and atmospheric transport modelling (ATM) of radionuclides to characterize the eruption as comprehensively as possible and discussed the results in the context of CTBT issues.

The Hunga volcanic eruption triggered large parts of the global IMS network. Thus, it was straightforward and uncomplicated to perform the most basic analyses, such as localization and investigating the waveforms in the time and frequency domain to search for explosion characteristics. The IMS network, meanwhile, is well established for all four technologies to sufficiently analyse events of interest even in such remote areas as the Tonga archipelago. It has been proven that this is also the case for less forceful events. We show that the IMS seismic station coverage for Hunga eruption is enough to produce consistent analysis results of similar quality in comparison to the more expansive publicly available FDSN network. Nevertheless, the additional use of data from stations not part of the IMS, especially from additional seismic stations, is of high value. They can complement the IMS data conveniently. Furthermore, data of other technologies (so-called national means in CTBT technical terms) such as GOES-17 satellite images and topography/bathymetric data (used here) or GNSS and InSAR data (not used here, due to most of the processes occurring underwater) cannot be valued highly enough to support the analysis of CTBT relevant events.

For National Data Centres, the CTBT does not stipulate specific analysis methods but describes key parameters (e.g. signal frequency, amplitude, depth, etc.) and analysis goals to stay flexible in case new and more efficient methods develop over time. However, with various test cases which occurred in the past, a set of forensic methods and tools evolved. These analyses include localization in space and time with additional waveform data, strength estimation via magnitude determination and thereupon a yield estimation, as well as a determination of the source mechanism for discrimination. If applicable array methods are popular for these processing steps. This is reasoned in their capability of improving the signal-to-noise ratio and acting as a directional antenna, hence, being a powerful tool even in case of inauspicious station coverage.

Here, we could show that the used methods are already well suitable to fulfil the tasks in general. There is no doubt in the readiness of the CTBTO verification regime, though not entered into force yet. However, specific issues in method details remain. They are well known from earlier studies but did not prevent reliable results due to application on comparably more simple source mechanisms. However, the specific mechanism of the Hunga eruption, having a well constrained point-like origin but being highly complicated in time and temporal-spatial source processes, highlighted that we could improve our methods. Especially, interactions and coupling between the different media solid Earth, water, and atmosphere are not very well understood in detail yet and we are missing appropriate models.

The issue of using magnitude relations calibrated for shear instead of explosive events used to be a minor one in previous studies but showed its drawbacks strongly for the Hunga eruption. It would be favourable to have magnitude/strength relations specifically adjusted for explosive events from all the three waveform techniques. Preferably, they would consider possible coupling of energy between the different propagation media. First approaches exist and could be further developed based on data from the Hunga eruption.

The discrimination of the source mechanism, crucial with respect to the CTBT tasks, can be at first order approached by seismic waveform inversion for the moment tensor. It is simple to distinguish a tectonic earthquake from an explosive source. However, the discrimination of the real nature of the explosion (e.g. nuclear test versus volcanic eruption) is not that simple, for example due to a trade-off between the non-shear parts of the moment tensor (ISO and CLVD). Moreover, the seismoacoustic coupling at the Earth’s surface is not yet considered sufficiently during waveform modelling, a joint problem of seismic and infrasound data analysis. An interesting task would be the joint inversion of SHI data to exploit their different sensitivity for source processes and balance their drawbacks. So far, this task is hampered by open questions concerning Green’s functions for infrasound and hydroacoustic data as well as details of the inversion setting, such as weighting schemes between data, applicable frequency range, etc.

The tremendous strength of the Hunga eruption and, hence, the vast amount of various data types with very good signal-to-noise ratio is an ideal opportunity to further improve existing methods and routines. Especially, when it comes to multitechnology analyses as it is established when investigating relevant events with respect to the CTBT, such a high-quality global data set provides ideal conditions for more in-depth examination.

Not only for method development, the Hunga eruption is a unique opportunity. A common feature of all presented technologies is that they observed the event from hundreds of kilometres distance. Yet, many details of the eruption processes, occurring on a relatively short timescale of a few hours, could be resolved. The results revealed a very complex pattern of processes that we could disassemble in remarkable detail and from different viewpoints based on different data types. We are confident that this unique data set will lead to an immense increase in knowledge-gain concerning volcanic eruptions.

ACKNOWLEDGEMENTS

All authors thank the CTBTO and the IMS station operators for guaranteeing the high quality of the data. We thank Shane Cronin for sharing details on his impressive work on the Tonga bathymetry before the data were published. We thank the editor Robin Matoza for efficient and supportive handling of the paper. We highly appreciate comments by two anonymous reviewers.

DATA AVAILABILITY

Data from regional seismometers are available via FDSN services from GEOFON and IRIS. Data from global IMS infrasound, seismic and hydroacoustic stations and arrays as well as radionuclide stations are available to National Data Centers of the CTBT and to others upon request through the virtual Data Exploitation Centre (vDEC) of the IDC at https://www.ctbto.org/specials/vdec. ECMWF products, including the atmospheric model analysis, being no longer valid for forecasting are made available via https://www.ecmwf.int/en/forecasts/dataset (last accessed 21 August 2022) under CC-BY 4.0 License. The NCEP-GDAS meteorological data as input for ATM simulations can be retrieved via https://www.ready.noaa.gov/archives.php. The dynamic Green’s function store is uploaded on the Pyrocko project Green’s mill repository https://greens-mill.pyrocko.org/ as ‘global₂s’ for reproduction. We make the Grond reports of our MT inversion results available as a zenodo repository (Donner et al.2022).

References of operators of seismological networks and stations used in this study

Alfred Wegener Institute For Polar And Marine Research (AWI). (1993). AW – AWI Network Antarctica. Deutsches GeoForschungsZentrum GFZ. https://doi.org/10.14470/NJ617293

Instituto Nicaraguense de Estudios Territoriales (INETER). (1975). Nicaraguan Seismic Network [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/NU

Hugh Glanville. (2019). Beetaloo Seismic Monitoring Project [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/2O2019

Arizona Geological Survey. (2007). Arizona Broadband Seismic Network [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/AE

Federal Institute for Geosciences and Natural Resources. (1976). German Regional Seismic Network (GRSN). Bundesanstalt für Geowissenschaften und Rohstoffe. F

Istituto Nazionale di Oceanografia e di Geofisica Sperimentale. (1992). Antarctic Seismographic Argentinean Italian Network - OGS [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/AI

Alaska Earthquake Center, Univ. of Alaska Fairbanks. (1987). Alaska Regional Network [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/AK

NOAA National Oceanic and Atmospheric Administration (USA). (1967). National Tsunami Warning Center Alaska Seismic Network [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/AT

H., G., & Geoscience Australia. (2021). Australian National Seismograph Network Data Collection (Version 2.0, September 2018) [Data set]. Commonwealth of Australia (Geoscience Australia). https://doi.org/10.26186/144675

Universidad de Chile. (2012). Red Sismologica Nacional [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/C1

Institut de physique du globe de Paris (IPGP), & École et Observatoire des Sciences de la Terre de Strasbourg (EOST). (1982). GEOSCOPE, French Global Network of broad band seismic stations. Institut de physique du globe de Paris (IPGP), Université de Paris. https://doi.org/10.18715/GEOSCOPE.G

Natural Resources Canada (NRCAN Canada). (1975). Canadian National Seismograph Network [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/CN

Albuquerque Seismological Laboratory (ASL)/USGS. (1980). US Geological Survey Networks [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/GS

Albuquerque Seismological Laboratory (ASL)/USGS. (1993). Global Telemetered Seismograph Network (USAF/USGS) [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/GT

Scripps Institution of Oceanography. (1986). Global Seismograph Network - IRIS/IDA [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/II

Albuquerque Seismological Laboratory (ASL)/USGS. (1992). New China Digital Seismograph Network [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/IC

Albuquerque Seismological Laboratory/USGS. (2014). Global Seismograph Network (GSN - IRIS/USGS) [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/IU

Rutgers University. (2013). Ocean Observatories Initiative [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/OO

Regional Integrated Multi-Hazard Early Warning System (RIMES Thailand). (2008). Regional Integrated Multi-Hazard Early Warning System [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/RM

Alberta Geological Survey / Alberta Energy Regulator. (2013). Regional Alberta Observatory for Earthquake Studies Network [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/RV

Australian National University (ANU, Australia). (2011). Australian Seismometers in Schools [Data set]. AusPass. https://doi.org/10.7914/SN/S1

UC Santa Barbara. (1989). UC Santa Barbara Engineering Seismology Network [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/SB

Información de la Red Sismológica Nacional de Costa Rica. (2017). https://doi.org/10.15517/TC

Bureau of Economic Geology, The University of Texas at Austin. (2016). Texas Seismological Network [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/TX

Albuquerque Seismological Laboratory (ASL)/USGS. (1990). United States National Seismic Network [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/US

University of Utah. (1962). University of Utah Regional Seismic Network [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/UU

University of Washington. (1963). Pacific Northwest Seismic Network - University of Washington [Data set]. International Federation of Digital Seismograph Networks. https://doi.org/10.7914/SN/UW

SUPPORTING INFORMATION

suppl_data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

References