Multi‐Event Explosive Seismic Source for the 2022 $Mw$ 6.3 Hunga Tonga Submarine Volcanic Eruption

On 15 January 2022, the submarine volcano Hunga Tonga–Hunga Ha’apai (HTHH), located in the Tonga–Kermadec volcanic arc archipelago, entered a violent eruptive phase. This major eruptive phase followed a month of lower level volcanic activity that began on 20 December 2021 when the first Surtseyan explosions of the edifice were observed (Global Volcanism Program, 2022a). This month‐long activity occurred while the volcanic edifice was coming out of a dormant period of nearly seven years following the much smaller 2014 eruption.

At 04:15 UTC on 15 January 2022, a series of explosions was heard on surrounding islands (Global Volcanism Program, 2022b), accompanied by visible volcanic activity. An eruptive plume briefly reached an altitude of 55 km (Carr et al., 2022), crossing the stratosphere–mesosphere interface. The eruption produced an atmospheric wave that traveled multiple times around Earth and was audible as far away as Alaska, 9370 km to the northeast (Matoza et al., 2022). The eruption was accompanied by a tsunami that impacted the Kingdom of Tonga, causing damages in the vicinity of the volcano, three causalities in the islands of Tonga, and two in Peru (Global Volcanism Program, 2022b).

This powerful series of volcanic explosions generated globally recorded seismic signals, which form the basis for our study. We assume a point‐source approximation and estimate moment tensors (MTs) using both regionally recorded surface waves and teleseismic body waves. Our results from surface waves show an explosive source component with a total moment magnitude of $Mw$ 6.14 for the single‐event inversion and $Mw$ 6.34 when considering a set of four subevents. Our results provide fundamental information on the source physics of the HTHH event, complementing available geophysical observations worldwide.

The seismic moment tensor is a 3 × 3 symmetric matrix representation of seismic radiation from a point source. The tensor’s eigenframe allows expressing the MT magnitude and source type via its eigenvalues and its orientation as strike, dip, and rake angles (Tape and Tape, 2012). Although MTs are most commonly used as point‐source representations for earthquakes, they can also be used to describe subevents within a larger finite‐fault rupture (Tsai et al., 2005) or “exotic” events such as landslides, volcanic events, explosions, cavity collapses, and glacier events (see Alvizuri and Tape, 2016; Table 1).

MT inversion is typically formulated as a data‐fitting procedure, which seeks to infer the MT solution $M‾$ that minimizes the least‐squares distance between synthetic seismograms and observed data by exploring the 6D MT space. We then assume that the MT that best explains the data is representative of the actual source that generated the recorded seismic ground motion.

In the following, we adopt a simplified notation without time and space dependencies, in which $G$ is an abstraction of the entire set of Green’s functions for a given source and a network of stations, and $dobs$ are the observed seismograms in that network.

We use the moment tensor uncertainty quantification (MTUQ) software package (see Data and Resources) to carry out the seismic MT inversion. MTUQ is an open‐source, parallel MT inversion code that allows for scalable and efficient grid searches, enabling detailed uncertainty quantification. MTUQ includes an implementation of a “cut‐and‐paste” misfit function (Zhu and Helmberger, 1996; Zhu and Ben‐Zion, 2013) and takes advantage of Python seismology libraries such as ObsPy (Krischer et al., 2015) and instaseis (van Driel et al., 2015) for processing data and generating synthetic waveforms. We adopt the uniform MT parameterization of (Tape and Tape, 2015), which is well suited for grid searches and allows for source type visualization and graphical representation of the parameter space.

To characterize the HTHH source process, we search over depth, magnitude, source type, and orientation. Our main search involves $∼109$ evaluations of the misfit function, using the following grid discretization:

• 11 depths ranging from $zc−5km$ to $zc+5km$ with a $Δz$ interval of 1 km, in which $zc$ is an initial depth determined from a coarse grid search.

• 21 magnitudes ranging from $Mc−0.2$ to $Mc+0.2$ with a $ΔMw$ interval of 0.02, in which $Mc$ is an initial moment magnitude determined from a coarse grid search.

• 800 source types spanning a 20 × 40 grid on the lune (specifically, the grid is on a v‐w rectangle).

• 8000 orientations uniformly sampled from the strike–cosine(dip)–rake “brick” (Tape and Tape, 2012).

A typical grid search of this size required 45 min on 244 cores (183 CPU‐hr). As described in the following sections, these grid searches were used together with a stochastic optimization algorithm to estimate the source‐time function (STF).

The synthetic waveforms used in our inversion were generated from a Green’s tensor database for the 1D layered Earth model AK135F, computed at 5 s dominant period using the spectral element solver AxiSEM (Montagner and Kennett, 1996; Nissen‐Meyer et al., 2014). The database itself was obtained from the Incorporated Research Institutions for Seismology Data Management Center Syngine web service and parsed using Instaseis for efficient buffered I/O (van Driel et al., 2015; Krischer et al., 2017). Overall, this approach allows for rapid generation of synthetic seismograms, making it possible to carry out large‐scale searches of MT space.

The initial hypocenter and origin time provided by the U.S. Geological Survey are longitude −175.390°, latitude 20.546°, 0 km depth, and 2022‐01‐15 04:14:45 UTC; this hypocenter is near the summit of the HTHH submarine volcano. For our inversion, we use publicly available three‐component broadband seismograms from all stations within 3000 km from the hypocenter in a time interval from 400 s before to 3000 s after the origin time. We use ObsPy to remove the instrument response and to rotate horizontal components into the radial and transverse directions.

Figure 1a shows the regional seismic stations used in this study, centered on the HTHH Volcano. A record section of vertical‐component seismograms (Fig. 1c) immediately reveals an unusual feature of the event: it appears to last for hundreds of seconds—much longer than what would be expected for an earthquake of comparable magnitude. In Figure 1c, we qualitatively identify several surface‐wave arrivals that cannot be explained by a single explosive source. Our MT analysis follows two approaches for the STF: a single‐source pulse and multiple pulses.

There are no clear high‐frequency body waves in the regional data from the HTHH event, which prevented us from using regional P wave (polarity or waveforms). Therefore, we focus on long‐period surface waves to invert the MT associated with the main event. All data are downsampled to 5 Hz and bandpass‐filtered between 25 and 70 s. Later, we were able to exploit teleseismic body waves in our study to refine the event origin time and verify the plausibility of our solution source.

After analyzing the HTHH event, we performed MT inversions for eight $Mw>5.8$ shallow earthquakes close to HTHH (Fig. 1b; Fig. S5, available in the supplemental material to this article) to demonstrate that it is possible to discriminate MT source types, depth, and magnitude using only surface waves and the sparse network of available stations.

We adopted the same processing for the earthquake data as with the HTHH event, and we used the same set of stations, whenever possible, to allow for straightforward comparisons with the earthquake control group. We performed the same grid search over source parameters for each event, allowing for differences in the central magnitude (⁠$Mc$⁠) and central depth (⁠$zc$⁠).

The complete MT inversion results for the HTHH event and eight comparison earthquakes can be found in the supplemental material. The MT that best explains the surface waveforms generated by the HTHH event is a shallow source with all positive eigenvalues, corresponding to a mechanism with a significant explosive (i.e., positive isotropic) component.

We tested a range of magnitudes between $Mw$ 5.8 and 6.2 with a 0.02 increment and depths from 0 to 5 km; the minimum misfit is found for $Mw$ 6.14 at 2.0 km depth (Fig. S2). The depth uncertainty is approximately 2 km, indicating that a near‐surface hypocenter is consistent with the surface‐wave data used in our study.

The misfit map represented on the eigenvalue lune in Figure S4a displays a swath of low‐misfit area ranging from the upper left to the lower right. This is characteristic of shallow sources with azimuthally symmetric surface wave radiation patterns, such as a horizontal crack tensor or an isotropic‐like tensor having all positive eigenvalues (with energy being radiated outward in all directions). The pattern can be seen in Ford et al. (2010, 2012) and is discussed in Alvizuri et al. (2018). The probability density function for source type (Fig. S4) conveys the concentration of likely solutions within the upper portion of the lune.

Our results for the eight catalog earthquakes and the HTHH event are shown on the eigenvalue lune in Figure 2. We are able to discriminate the HTHH event, at the upper left of the eigenvalue, from the earthquakes, which are closer to the double couple at the center of the lune. Some of the inverted sources have source types that are far away, as measured by the $θ$ angle in Table S1 from the double couple; this amount of spread is typical for full MT estimates of earthquakes (see also fig. 2 of Boyd et al., 2015, and fig. 3d of Alvizuri et al., 2018). In addition, our depth search results in values that are generally within a few km of Global Centroid Moment Tensor catalog depths, which are based partly on body waves (Table S1).

Secondary wavepackets with a similar linear moveout as the first wavepacket (subevent S1) are clearly visible in the observed data, as highlighted by the second gray band in Figure 1c. Based on observations of long‐distance P‐wave arrivals (from 70° to 90° distance), Yuen et al. (2022) interpreted four similar‐type events within the first 400 s. Therefore, we chose to estimate a source‐time function having up to five subevents. Our formal exploration of model parameter space, described next and highlighted in Figures S8 and S9, resulted in only four subevents.

We invert the STF of the assumed HTHH explosive sequence by performing a search over onset time and relative amplitude of each of the subevents, defining a sequence of identical explosive sources. Each trial STF $Γ(μ,α)$ is defined as a sum of i = [1,2,3,4,5], fixed‐variance Gaussian functions with means $μi$⁠, which define the onset time for each subevent, and scaling factors $αi$⁠, which define the relative magnitude of each event. The variance for each of the five Gaussian function is fixed at $0.15s2$ to emulate short, impulsive sources. Finally, the sum of these five functions is normalized, such that the integral over the STF goes to 1 to preserve moment magnitude scaling.

To carry out the STF inversion, we employ the covariance matrix adaptation evolution strategy (CMA‐ES), which is a derivative‐free, iterative optimization method that has proven to be efficient on a large set of misfit functions (Hansen and Ostermeier, 2001). When inverting the STF, we limit the data to the vertical and radial components from a subset of the five closest stations, such that we can verify later the validity of our STF on the complete dataset. We fix the source mechanism $M‾$ to the one estimated for the S1 subevent; therefore, the model parameter space contains 11 parameters: the relative amplitude (ranging from 0 to 1, allowing for events to vanish or merge if not required) and onset time for each of the five subevents and the total magnitude of the sequence.

The STF that minimizes the waveform misfit is shown in Figure 3a, and a full breakdown of the four subevents (labeled S1–S4) can be found in Table 1. The resulting estimated STF reveals that four subevents are needed to explain the observed waveforms, and that subevents S2, S3, and S4 occurred 203.22, 238.83, and 266.03 s, respectively, after the first subevent (Table 1). Uncertainties in these onset times are approximately +/− 1.0 s, based on the long‐period waves used.

We then used the inverted STF in our MT inversion, with the same settings as in the previous section, except using a longer time window of 550 s to accommodate the two separate wavepackets (S1 and S2–S4). This means that we assume a source that is multi‐part in time (S1 + S2 + S3 + S4) but occurring with the same mechanism and hypocenter for each subevent, and now only the mechanism (including magnitude) is estimated. The result of this multi‐event MT inversion is shown in Figure 3b,c. We found that the complex waveforms can be explained by a repeated MT with a total moment magnitude of $Mw$ 6.34, and we achieve satisfactory waveform fits in all components across the network, even though the transverse data were not used to constrain the STF. The STF obtained with the CMA‐ES inversion exhibits similarities in relative timing and amplitudes with the body‐wave displacement waveform stacks seen in Yuen et al. (2022). The equivalent magnitude of S1–S4 is annotated next to each Gaussian peak on Figure 3a. Using the multi‐event STF, the magnitude estimate of S1 changes only slightly (⁠$Mw$ 6.14 to 6.04).

We considered two additional dimensions of the seismic analysis: (1) teleseismic P waves and (2) point force instead of (point) MT. Teleseismic P waves, such as those analyzed by Yuen et al. (2022), are shown in Figure S11b for stations with epicentral distances from 70° to 90°. As shown in Figure S12, we were unable to fit both the radial and vertical components of the P waves using any point force. We were also unable to fit the regional surface waves using any point force (Figs. S12c and S13). For both the regional surface waves and the teleseismic P waves, a point‐source MT provides better waveform fits than a force.

Figure S14 shows the best MT result for the teleseismic body‐wave inversion, with the same set of stations considered previously. Although the best‐fitting source type has a large isotropic component—similar to the surface‐wave inversion—there is a notable difference in the estimated magnitudes. Using teleseismic body waves, we obtain $Mw$ 6.48 for subevent S1 (Fig. S14); using regional surface waves, we obtain $Mw$ 6.14 (Fig. S3). Further work is needed to understand this discrepancy. One possibility is that the attenuation used in our 1D model (AK135F) needs to be much higher in the shallowest layer; this would have the effect of increasing the estimated amplitude. A second possibility is that the assumed STF for the body‐wave modeling is more complex than the simple Gaussian function that we used.

Our primary motivation is to estimate the depth, magnitude, moment tensor, and source‐time function of the HTHH event to provide fundamental characteristics of this exceptional event—exceptional in terms of tectonic, volcanological, acoustic, and societal impacts (Matoza et al., 2022). A secondary motivation arises from a second type of shallow, explosive, and societally relevant event—nuclear explosions. Recent studies have demonstrated the value of regional surface waves and teleseismic body waves for characterizing explosions (Ford et al., 2012; Chiang et al., 2014; Chiang et al., 2016; Alvizuri et al., 2018; Walter and Wen, 2018). Moment magnitude can be used to estimate the explosive yield of a nuclear explosion (Pasyanos and Chiang, 2022); similarly, it can be used to characterize the size of a volcanic eruption (Gudmundsson, 2014) and the conditions leading up to the eruption.

Our characterization of source parameter uncertainties is enabled by high‐performance computing (HPC) and effective visualization approaches. By running the MTUQ code on an HPC cluster, we can finely search the designated space of model parameters, including depth and magnitude.

Figure 4 provides intriguing examples of how the misfit function varies within a subset of MT space. Figure 4a shows how the best‐fitting MT changes for each depth. Although the binary coloring scheme of the focal mechanism suggests isotropic mechanisms for most depths, we know that there are amplitude variations for each of the all‐solid focal mechanisms, and we can see from the white mark for the 0 km depth focal mechanism that it has one negative eigenvalue. For the HTHH event, we see that the magnitude estimate changes slightly with depth, varying between 6.20 and 6.40 for depths below 10 km.

Figure 4b,c shows eigenvalue lune plots for the best‐fitting $Mw$ (Fig. 4b) and depth (Fig. 4c) and considering all possible orientations. These plots should be compared with the misfit function plot in Figure 3b. Taken together, the lune plots provide insights into the trade‐offs among the estimated source type (on the lune), magnitude, and depth. For example, we see a complicated pattern for $Mw$⁠, indicating that the highest magnitude estimates occur near the global minimum source type (−19°, 59°), which has $Mw$ 6.34. The pattern for depth is also complicated, with greater depths occurring near the global minimum source type, which has 3 km. The darker regions on the lune correspond to MTs associated with shallower depths and, more importantly, higher waveform misfit (Fig. 3b). The variations in Figure 4b provide a basis for our assignment of 0.10 (⁠$1σ$⁠), resulting in our final magnitude estimate of $Mw$ 6.34 ± 0.10.

The misfit function displayed in Figure 3b conveys a large range of MT source types that provide adequate fits to the observed waveforms. Nevertheless, it exhibits two features distinct from the earthquake results: (1) The misfit pattern on the lune is not centered on the double couple (Fig. S6), and (2) the best‐fitting MT is far from the double couple, which is at the center.

Having found that an explosion‐like mechanism is well constrained, we are interested in providing a more informative visualization than the all‐solid focal mechanism in Figure 3b. Two options are provided in Figure 4d. The first is a properly oriented ellipsoid for which semiaxes are proportional to the (positive) eigenvalues. We see a subvertically oriented T‐axis for the largest eigenvalue. The second display shows a sphere with shading proportional to the radial component and with integral curves (or flow lines) that represent the radiation pattern. In this depiction, the bending arrows convey that some shear motion will occur, as expected for any MT that is not exactly isotropic.

Visual representation of MTs with all‐positive eigenvalues is valuable for understanding differences among MTs in the uppermost and lowermost regions of the lune. Nevertheless, it is also important to realize that the differences between MTs in these regions—measured by the angle between matrix MTs—are relatively small, and therefore it is challenging to discriminate between MTs using regional seismic data. For the HTHH event, a variety of shapes of ellipsoids—from “pancakes” to “cigars”—are consistent with the recorded waveforms; furthermore, non‐all‐solid focal mechanisms are permissible, as indicated by the misfit plot in Figure 3b (see also Fig. S7).

Seismic surface waves generated by the HTHH event constrain fundamental aspects of the first minutes of the peak eruptive phase. Our MT analysis characterizes the seismic event as a sequence of at least four shallow explosions starting at 04:15:17 UTC, spanning ∼270 s, and having a cumulative magnitude of $Mw$ 6.34 ± 0.10 (Table 1). In comparison, the globally recorded acoustic waves from this event, subject to large uncertainties on acoustic‐elastic coupling and atmospheric properties, provide less direct constraints and do not meaningfully refine the seismically estimated source parameters. Constraints from the seismic analysis—notably the relative timing, size, and mechanisms of the subevents—should be beneficial in establishing a comprehensive characterization of the HTHH eruption.

Seismic data were obtained from the Incorporated Research Institutions for Seismology Data Management Center (IRIS DMC). We use stations from the Geoscope network (G, http://geoscope.ipgp.fr/networks/detail/G/), the Global Seismograph Network (II, DOI: 10.7914/SN/II, and IU DOI: 10.7914/SN/IU), the Australian National Seismograph Network (AU, http://auspass.edu.au/fdsnws/dataselect/1/), and the New Zealand National Seismograph Network (NZ, http://service.iris.edu/fdsnws/dataselect/1/). The moment tensor code used in this study, MTUQ, is openly available at https://zenodo.org/badge/latestdoi/483131450 and under development at https://github.com/uafgeotools/mtuq. The Global Centroid Moment Tensor (Global CMT) catalog (Dziewonski et al., 1981; Ekström et al., 2012) is accessible via https://www.globalcmt.org/. All websites were last accessed in September 2022. The U.S. Geological Survey (USGS) National Earthquake Information Center (NEIC) catalog origin time of 2022/01/15 04:14:45 UTC and magnitude of $Mw$ 5.8 were obtained from https://earthquake.usgs.gov/earthquakes/eventpage/us7000gc8r (last accessed August 2022). The supplemental material is a single pdf file containing three parts with supplemental text, 1 table, and 14 figures.

The authors thank two anonymous reviewers and editors Steven Gibbons and Keith Koper for feedback that improved this article. The authors thank the University of Alaska Fairbanks (UAF) Research Computing Systems for an allocation on the “chinook” high‐performance computing cluster. The authors benefited from discussions with Jeremy Pesicek (U.S. Geological Survey), David Fee (UAF), and Walt Tape (UAF). This project was supported by Air Force Research Laboratory Contract HQ0034‐20‐F‐0284 funded through the UAF Geophysical Detection of Nuclear Proliferation (GDNP) University Affiliated Research Center (UARC).