Episodic Magma Hammers for the 15 January 2022 Cataclysmic Eruption of Hunga Tonga‐Hunga Ha'apai

Understanding the forces and magma system dynamics on timescales of seconds to minutes remains challenging. In the January 2022 phreatoplinian Hunga Tonga‐Hunga Ha'apai eruption, four remarkably similar seismic subevents within a 5‐min interval occurred during the intensifying early eruptive phase. The subevents are similar in waveforms and durations (∼25 s each). Each subevent begins with an unusual negative P‐wave polarity which is inferred, using full‐wave seismic modeling, to be caused by an upward single‐force mechanism at the volcano created by a magma hammer likely in response to magma flow blockage/constriction during the early part of the eruption as discharge rapidly increased over orders of magnitude with concomitant conduit geometry evolution and instability. Our proposed episodic magma hammer model is consistent with thermodynamic and phase properties of the magmatic mixture, and yields an estimate of conduit mass flow in agreement with vent discharge rates derived from satellite imagery of plume heights.


Timings of Various Volcanic Phenomena
Historically, HTHH is a large submarine edifice that has experienced small-scale submarine and Surtseyan activity (Cronin et al., 2017;GVP, 1988GVP, , 2009aGVP, , 2009bGVP, , 2015GVP, , 2022a) ) between large, caldera-forming eruptions (Brenna et al., 2022) (see Figure 1).The latest period of activity began on 20 December 2021, ultimately culminating in the cataclysmic 15 January 2022 eruption (Figure 1b).An eruption on 13 January at 15:20 UTC (04:20 14 January local time) produced a 20 cm tsunami and a 20 km tall plume, accompanied by ∼190,000 lightning flashes (GVP, 2022b)-and even partially destroyed a land mass formed during the 2014-2015 eruption that connected the islands of Hunga Tonga and Hunga Ha'apai (Cronin et al., 2017).An overnight eruption totally removed this land mass in the early morning hours of 15 January, before the final, most intense eruptive phase began.
The 15 January volcanic eruption did not begin abruptly.Instead, an ensemble of time-transgressive phenomena overlapped (Figure 1c).The USGS reported a M5.8 seismic event with an origin time of 04:14:45 15 January 2022 (UTC), located at 20.546°S, 175.390°E with 0 km depth.This time cannot be the start time of the eruption, as satellite imagery captured at 04:10 UTC shows a rapidly growing volcanic plume already 18 km high-eventually peaking at mesospheric heights of ∼58 km around 04:30 UTC (Yuen et al., 2022).Extrapolating the infrasound signals in time and distance back to the volcanic location gives an eruption inception time around 04:02 ± 1 UTC (Yuen et al., 2022), which is more consistent with observed plume heights and the recent IMS infrasound analysis (Vergoz et al., 2022).Residents of nearby islands reported ashfall persisting throughout the night, suggesting an eruption duration of ∼12 hr (Yuen et al., 2022).Although the precise timing of the symphony of phenomena remains uncertain, an uptick in lightning intensity around 04:12 UTC is consistent with a strong eruptive pulse starting at 04:07-04:09 UTC (Figure 1c), allowing for an estimated 5-min delay between plume ascent and lightning onset, consistent with previous studies of volcanic lightning (Behnke et al., 2013;Hoblitt, 1994).The rapid increase in lightning after 04:12 UTC is consistent with high mass flow at the vent around 04:08 UTC.These phenomena indicate that the eruption was well underway before the 04:15 UTC seismic event analyzed in detail below.

Data and Results
Our seismic analysis involves minimal processing of the raw records.We first downloaded seismic waveforms recorded worldwide from Incorporated Research Institutions for Seismology (IRIS), removed their instrumental responses, and high-pass filtered traces with a corner frequency 0.01 Hz.Theoretical travel times of the direct P waves for all the stations using the USGS reported origin time and location and a crosscorrelation were used to align the waveforms (Figure 2b) (Yu et al., 2017).We selected 417 stations (Figure 2a) based on the signal-noiseratio (SNR > 1.6) defined as the ratio of the root-mean-squares amplitude in the P-wave signal window (within 60 s following the P time) and the noise window (within 60 s before the P time).Shifted waveforms were stacked to obtain the seismic wavelet used to study the source process (Figure 3).

Unusual Repeating Episodic Seismic Signals Recorded Worldwide
Within a time of ∼300 s immediately after the reported seismic event, four subevents-E-1, E-2, E-3, and E-4-are visible in both the stacked ground displacement and velocity seismograms (Figures 3a and 3b).Each subevent has a similar duration of ∼25 s.If these were regular tectonic earthquakes, a 25 s source duration would result in an event of around Mw7.5, far greater than the reported magnitude M5.8 by USGS.A 204.6 s time interval occurs between subevents E-1 and E-2, ∼39.8 s between E-2 and E-3, and ∼25.4 s between E-3 and E-4; note that the tail end of E-3 overlaps the beginning part of subevent E-4 (Figure 3).Each subevent can be further divided into four consecutive eruption stages A/B/C/D (Figures 3c-3e) with respective durations 3 s/15 s/3 s/3 s.We call the subevents E-1 to E-4 "episodic" because they have similar waveforms Figures S1 and S2 in Supporting Information S1), indicating that similar eruption dynamics likely governed these subevents.Information S1).Secondly, the azimuthally stacked P waveforms are highly similar which shows the radiated P-wave is mostly azimuthally symmetric (Figure S2 in Supporting Information S1).Thirdly, this symmetry property is further supported by a lack of visible energy in the transverse component in the stacked S waves (see Figure S2 in Supporting Information S1).These characteristics suggest that the subevents were not caused by earthquake dislocation faulting or landslides on the volcano flank, but rather by forces with mostly azimuthal symmetry to the first order.
An azimuthally symmetric source consistent with negative polarity could be: a single force in the vertical direction, or an implosion, or a CLVD source (Knopoff & Randall, 1970) at the volcano.Using the stacked wavelet as the source time function (Figure 3a), our detailed full-wave seismic modeling (see Text S1 in Supporting Information S1) showed single forces in the vertical direction are preferred over either implosion or CLVD because (a) the synthetic waveforms can fit the observed data for both body waves (Figure S3 in Supporting Information S1) and surface waves (Figure S4 in Supporting Information S1) in both vertical and radial components; (b) waveform synthetics due to an implosion/explosion or a CLVD source show strong P-to-S converted waves at ∼600 s, which are not observed in the seismic records (see Figures S5 and S7 in Supporting Information S1), despite the fact that they can fit the surface waves well (Figures S6 and S8 in Supporting Information S1).
We also investigated if the single forces were truly vertical by examining the subtle azimuthal variations in the locally averaged P-wave amplitudes of E1, E2, and E3 (Figure S9 in Supporting Information S1).Similarities in the azimuthal trends of E1, E2, and E3 amplitudes strongly suggest that the four subevents have an identical or closely related physical origin.We find that a tilted force fits the observed azimuthal P-amplitudes better (Figure S9 in Supporting Information S1).The tilted single force is along the line containing the vector [1, 0.08, 0.46] in [Upward, North, West] directions, respectively.This force has a ∼25.0-degree tilt with respect to the true vertical direction.Furthermore, regional surface wave modeling using this tilted single-force and stacked source wavelets (Figure 3) provided a satisfactory match between modeled and observed waveforms in all three components (Figure S10 in Supporting Information S1).Note that we used the 1-dimensional PREM (Dziewonski & Anderson, 1981) model in the full waveform modeling (not taking into account the heterogenies in the Earth) so perfect matching between the data and the modeled waveforms is not expected.The explosion mechanism explored by Thurin et al. (2022) is not consistent with data in two aspects: the negative P-wave first motion polarity and the above-mentioned converted S waves in the modeled waveforms ∼600 s.Thurin et al. (2022) were not able to model multicomponent regional surface waves and body waves using single forces whereas we are able to fit both types of multicomponent data.Their poor waveform fitting could be caused by their use of Gaussian wavelets which have only one polarity.However, it is evident that the stacked source wavelets are not Gaussian (Figure 3).Poli and Shapiro (2022) used single force to estimate the Volcanic Explosivity Index.Garza-Girón et al. ( 2022) set the single-force history to be downward a priori and they considered our E2∼E4 not as independent episodic events.We are not to say that during the ∼12 hr HTHH eruption explosions did not take place, but only that for the records studied here, explosions can be excluded.(GVP, 1988(GVP, , 2009a(GVP, , 2009b(GVP, , 2015(GVP, , 2022a(GVP, , 2022b) ) and by our analysis of available satellite imagery (e.g., Yuen et al., 2022).For days where no eruptive activity was noted, volcano symbols show no eruption.Days with eruptive activity but no substantial eruption column are marked by lava on the volcano symbols with Surtseyan eruptions, and eruption columns are to-scale for those days where substantial plume development was reported.Volcanic lightning data are also shown in black dashed line with lightning icons (Yuen et al., 2022).10.1029/2023GL102763 4 of 13

Eruption Stages and the Magma Hammer Mechanism
As the lack of subsurface information precludes obtaining well-constrained parameters for any HTHH dynamic model, the purpose of our simplified analysis is to present a mechanism which aligns with the unusual observed seismicity of the 15 January eruption and the calculated forces involved, while also agreeing with the observed development of the eruption plume and magma properties.Based on the stacked P-wave seismograms (Figures 3b-3e), we attempt to decipher the four eruptive stages (A, B, C and D) in the context of conduit flow through magma-conducting conduits for each of the subevents.Because each of the seismic subevents exhibits similar waveforms, we focus our analysis on the E-1 subevent time series (Figure 4).The respective durations for Stages A to D are ∼3 s, 15 s, 3 s, and 3 s, with correlating the single-force directions of up/down/up/down, respectively (Figure 4b).Volcanic seismic sources usually give rise  Chouet & Matoza, 2013).Single-forces in the up-down direction rather than resolved force-couples are rare.In this regard, this HTHH eruption shares some striking similarities to the 2004 Mt.Asama activity (the most active volcano in Honshu Japan), which has well-resolved down/up/down single-forces in three stages for five episodic similar events that occurred over a ∼2.5 month period (Ohminato et al., 2006).Like HTHH, Mt.Asama also generated a strong air shock (Ohminato et al., 2006) indicative of a strong under-expanded volcanic jet with a vent magma pressure exceeding ambient atmospheric pressure.In the following, we use well-established ideas to discuss the four stages A∼D we have found.Although there is no doubt phreatomagmatic activity can play a role in generating local and regional seismic signals, its role on teleseismic signal generation is probably not very significant for the following two reasons.First, the earlier Jan 14 HTHH eruption generated no coherent global teleseismic signals, although strong phreatomagmatic activities and explosions occurred.Second, the 2004 eruption of Mt.Asama involved neither phreatic nor phreatomagmatic activity (GVP, 2004).
Mechanisms for producing upward forces include: viscous drag force along the conduit wall due to a steady state upward magma motion (Coppess et al., 2022;Ohminato et al., 2006), magma return flow at some barrier due to a uprising gas plug push (B.Chouet et al., 2003) (similar to the water-hammer effect), and constriction of flow (Morrissey & Chouet, 1997) akin to hydraulic shock (Kieffer, 1984;Morrissey & Chouet, 1997;St. Lawrence & Qamar, 1979;Streeter & Wylie, 1974;Yin, 2018).In the conventional 1-D water-hammer mechanism, water flows steadily through a pipe with an open valve.Rapid closure of this valve creates a surge of water hammer pressure P H , acting on the closure point which can sometimes burst a steel pipe.The effect can also be due to a piston of fluid impacting on an obstructive barrier (B.Chouet et al., 2003) such as a volcanic plug (Iverson et al., 2006) or flowing through a constriction (Morrissey & Chouet, 1997).We do not know whether such elements exist or not in the HTHH subsurface.However, our rudimentary magma-hammer model yields a remarkably consistent estimate of magma flux comparable to the vent discharge rate.Here, we sketch the basic water-hammer concept (necessarily simplified due to lack of subsurface information), modified to the parameters of a multi-phase magma.
Mathematically, the water-hammer pressure, P H , can be described by the well-known Joukowsky equation (Streeter & Wylie, 1974) for slow fluid flow (V << c): P H = ρcV, where ρ is the fluid density, c is the sound speed in the fluid, and V is the fluid flow speed in the pipe.Since magma is involved here rather than water, we apply the term magma hammer to describe this process.

Discussion of Forces in the Eruption Stages
Before Stage-A (time <15 s in the record), the stacked seismic record is characterized by fluctuating ground motion (Figure 4b).Based on infrasound wave speeds and satellite imagery (Yuen et al., 2022) the magma transport conduit had already connected to a surface vent and the plume had likely already been established.

Stage A: Magma Hammer and Estimates of Magma Flux and Discharge
Stage-A (E-1A) lasted about ∼3.2 s and the farfield ground motion was characterized by a distinct downward ground motion, meaning that the seismic source was likely to be an upward single-force (Figure S11c in Supporting Information S1).Full-wave seismic modeling (Text S1 in Supporting Information S1) reveals that this upward single-force source has a force magnitude ∼5.0 × 10 12 N (Text S1 in Supporting Information S1; Figure S3-S8).We argue that Stage A of E-1∼E-4 with similar ground displacements represents magma hammer breaching of a short-lived flow blockage of the magma conduit such as a caldera plug, generating an upward single force.

Magma Hammer and Piston Motion
Multiplying both sides of the Joukowsky equation  ( =  ) by the cross-sectional area, A, of the flow conduit, we get F H = P H A = ρcQ V = cQ m , where the hammer force, F H , can be directly estimated by the seismic modeling (∼5 × 10 12 N for Stage-A; Figure S11b in Supporting Information S1), Q V = VA is the magma volume flux (in m 3 /s), and Q m = ρVA is the magma mass flux (in kg/s), in the subsurface plumbing system.It is important to note that in this expression, magma density is that of the magmatic mixture of melt plus fluid (see Texts S2 and S3 in Supporting Information S1 for details).

Magma Equation of State and Flux Estimates
We now consider magma composition, degree of volatile saturation, and flow regime.Equilibrium crystallization of the last-erupted HTHH andesites was modeled using rhyolite-MELTS (Gualda & Ghiorso, 2015) to approximate the state of an andesitic magma as it undergoes closed-system degassing and ascends adiabatically to the surface from a reservoir depth of ∼7.5 km (P = 2 kbar).The state of the magma body was recorded at each P-T along an isentropic path and the volume fraction of the exsolved volatiles calculated (Text S2 in Supporting Information S1; Figure S12 in Supporting Information S1).The sonic velocity of a two-phase (melt + fluid) mixture was then computed for both bubbly flow (i.e., well mixed) and slug flow (i.e., gas pockets) regimes (see details in Text S3 in Supporting Information S1, Figures S13 and S14 in Supporting Information S1).We argue that slug flow is more likely, given the explosiveness of the eruption and short time scales involved although in a complex volcanic eruption it would seem reasonable to assume that both flow regimes may exist simultaneously at different locations and different times.Taking into account uncertainties, the sound wave speed in magma is estimated to be, c ∼ 1,000−2,000 m/s, and we can get a mass flux Q m ∼ 0.5 − 5 × 10 9 kg/s using the peak force in E-1A.Approximating the seismic wavelet in E-1A as a triangle, the average force should be halved.Alternatively, we can estimate the volcanic discharge into the atmosphere (Carey & Bursik, 2015).According to satellite images of the volcanic plume heights, the mass flow to support a 58 km volcanic plume from Yuen et al. (2022) is ∼1.4 × 10 9 kg/s.Our proposed magma-hammer model is consistent to within an order of magnitude or better with the force estimated from the seismic data, estimated vent discharge estimate into the atmosphere, physiochemical magma properties, and sonic velocities.

Stage B: Eruption of Volcanic Jet
Stage E-1B-a downward single force at the volcano-is the longest stage (15 s) in all four stages (Figure 4a).It is likely due to the reactive impulsive force associated with opening of the magmatic conduit and transient reestablishment of intense outflow for ∼15 s (Figure 4a).The eruption could have sent a series of pulses of volcanic product (ash, pyroclasts, crystals, lithics and seawater-derived steam) into the atmosphere based on observed maximum plume heights referred to earlier.This "upside-down rocket"-type eruption created a downward single-force on the order of ∼9 × 10 12 N (Figure S11b in Supporting Information S1).The magnitude of the force is comparable to that of the Mount St. Helens 1980 eruption (Kanamori & Given, 1982) of ∼10 13 N. Explosive magma fragmentation, which is a complex process, could decrease the magma viscosity and decouple the upward viscous drag force along the conduit wall to create a downward force (Coppess et al., 2022;Ohminato et al., 2006).In Stage-B, the net force first increases then decreases and both upward force and downward forces may coexist.If this is the case, seismic data and modeling show that the downward forces dominate toward the early part of the stage.

Stages C and D
Stage E-1C, initiated about 20 s after the onset of E-1A, exhibits an upward force of ∼4 × 10 12 N (Figure S11b in Supporting Information S1), and the downward force in E-1D is ∼3 × 10 12 N.We associate Stage-C with a second-phase closing/constricting of the conduit and the subsequent generation of an upward magma-hammer force or viscous drag upward force (details in Text S3 in Supporting Information S1).Stage E-1B exerted a sustained downward force, probably depressing both the less-fragmented magma in the conduit (Ohminato et al., 2006), the magma chamber and surrounding rock.At the end of E-1B, the force has diminished significantly and the possibility of the conduit elastically shrinking inward rapidly nearest the surface increases.Gravity can also assist this conduit closure as the second-order azimuthal variation in P-wave amplitudes suggests the conduit may be slightly tilted to the west with respect to the vertical (Text S1.3 in Supporting Information S1; Figure S9 in Supporting Information S1).The westward tilt is consistent with the post-eruption earthquake locations (Kintner et al., 2022).It is likely that Stage E-1C is related to the magma backflow impacting the overburden rock containing the constricted conduits, which creates an upward magma-hammer force and an upward viscous drag force due to the flow velocity.E-1D seems to be a reflected magma-hammer signal traveling down, which created a downward force as viscous drag or impact of the magma chamber.Stages E-1C and E-1D represent one cycle of oscillatory piston motion as magma motion travels through the conduit, from the reservoir to the surface.After E-1D (>45 s in Figure S4b in Supporting Information S1), we can see regular oscillating piston motions with a period T ∼ 9 s (Figure S11b in Supporting Information S1).Taking an estimated average traveling sonic velocity of c ∼ 1,000-2,600 m/s (Figure S13 in Supporting Information S1), we can determine the depth of the magma chamber at cT/4 ∼ 2.25 − 5.85 km using the formula given by Kieffer (1984), consistent with the geobarometric estimates (Brenna et al., 2022).

Episodic Seismic Events
Stages A-D for the E-1 event were repeated at least in E-2 and E-3 and likely also in E-4; E-4 is obfuscated by E-3, hindering our ability to clearly read the entire E-4 record (Figure 3).These large magnitude forces are short in duration and could produce the loud audible sounds as recorded by observers and infrasonic sensors around the world (Yuen et al., 2022).The repeated large magma forces can fracture the caldera basement and ballistically remove fractured rocks to excavate the caldera deeper and deeper, eventually forming a 600 m depression as shown in a recent bathymetry map (Cronin and Tonga Geological Services, 2022).The time intervals between the subevents (∼204.6,∼39.8, and ∼25.4 s) are long enough to allow for the conduit to be plugged again for additional magma hammers associated with geometric instabilities within the magma transport system.
The explosive HTHH eruption reveals a rapid but regulated nature of shallow magma transport dynamics manifested as episodic magma hammers.Examination and modeling of the seismic data suggest that HTHH produced the largest known magma-hammer force associated with a violent volcanic eruption.Future work will benefit from new 3D marine seismic imaging of the plumbing system to forecast future catastrophic eruption using the episodic nature of this system.

Conclusions
The 15-Jan-2022 HTHH eruption produced global seismic wavefields.Four episodic seismic events have been identified within a time window of ∼300 s in the intensifying phase of the volcanic activity.The four subevents have similar waveforms and durations.A westward tilted single force mechanism can fit well both the teleseismic waveforms and the regional surface waves.The magma hammer is one possible mechanism to cause the single force.The magma hammer also provides a way to estimate the conduit magma mass flux which is consistent with the vent discharge estimated from satellite imagery and magma properties.

Figure 1 .
Figure1.Timeline of events leading up to and including the onset of the 15 January 2022 eruption of Hunga Tonga-Hunga Ha'apai(HTHH).Eruption column heights (labeled by numbers) and eruptive history are reported by the Global Volcanism Program(GVP, 1988(GVP,  , 2009a(GVP,  , 2009b(GVP,  , 2015(GVP,  , 2022a(GVP,  , 2022b) )  and by our analysis of available satellite imagery (e.g.,Yuen et al., 2022).For days where no eruptive activity was noted, volcano symbols show no eruption.Days with eruptive activity but no substantial eruption column are marked by lava on the volcano symbols with Surtseyan eruptions, and eruption columns are to-scale for those days where substantial plume development was reported.Volcanic lightning data are also shown in black dashed line with lightning icons(Yuen et al., 2022).

Figure 2 .
Figure 2. Global seismic signals.(a) Distribution of global seismometers (blue triangles) used in this study and the Hunga Tonga-Hunga Ha'apai volcano (red dot); (b) Aligned seismograms based on the picked direct P-wave travel times showing four global subevents.The seismograms are vertical-component ground velocities, highpass filtered with a corner frequency of 0.01 Hz.The stations are aligned by the epicentral distance.

Figure 4 .
Figure 4. Eruption stages and magma plumbing system.(a) Volcano-reservoir schematics of states and forces (red arrow: magma hammer; black arrow: reaction force) corresponding to Stages A-D.The magma reservoir is centered at depth of 5-8 km based on geobarometry (Brenna et al., 2022).(b) Globally stacked ground velocity (vertical component) of event E-1.The amplitude of the waveform is not linearly proportional to the force history (Figure S11b in Supporting Information S1).Time zero corresponds to ∼04:14:45 15 January 2022.