Seismic Precursors to the Whakaari 2019 Phreatic Eruption detected in Five Other Volcanoes

Volcanic eruptions that occur without warning can be deadly in touristic and populated areas. Even with real-time geophysical monitoring, forecasting sudden eruptions is dicult because their precursors are hard to recognize and can vary between volcanoes. Here, we describe a general seismic precursor signal for gas-driven eruptions, identied through correlation analysis of 18 well-recorded eruptions in New Zealand, Alaska and Kamchatka. We show that the displacement seismic amplitude ratio, a ratio between high and medium frequency volcanic tremor, has a characteristic rise in the days prior to eruptions that likely indicates formation of a hydrothermal seal that enables rapid pressurization. Applying this model to the fatal 2019 eruption at Whakaari (New Zealand), we identify pressurization in the week before the eruption, and cascading seal failure in the 16 hours prior to the explosion. This method for identifying and proving generalizable eruption precursors can help improve short term forecasting systems.

Scolamacchia and Cronin, 2016); and/or active sealing by pressure applied to compressible clays (Kanakiya et al., 2021). Permeability barriers become even more important if gas input increases during magmatic unrest. Restriction of gas ow through the upper vent system leads to localized shallow pressurization, creating the conditions for an explosive eruption.
Seismic data analysis is standard practice at volcano observatories to track magmatic processes, volcano state and the possibility of future eruptions (Chouet, 1996 volcanoes (Gunawan et al., 2017), tremor prior to phreatic or phreatomagmatic eruptions are di cult to distinguish, or only identi ed in retrospect (Christenson and Wood, 1993;Christenson et al., 2010;Pardo et al., 2014;Dempsey et al., 2020;Yamada et al. 2021). This is partly due to the very shallow origin of these "small" explosions, with pressurization and triggering within a few tens to hundreds of metres of the surface (Montanaro et al., 2016;Jolly et al., 2020;Caudron et al., 2021). Signals are also dampened by formation of low-permeability seals (Christenson and Wood et al., 1993;de Moor et al., 2019). Quietening of a system can impart a false sense of safety if uids are being trapped/pressurised. A sudden pulse of additional magmatic gas may rupture these systems, but eruption initiation can also be any process that causes decompression, such as surface-propagating drying/cooling cracks, tectonic fracturing, sudden overburden loss due to landslides (Montanaro et al., 2016) or lake breakouts. Here, we analysed tremor time-series prior to 18 eruptions at six volcanoes: Whakaari, Ruapehu and Tongariro (New Zealand), Veniaminof and Pavlof (Alaska, USA), and Bezymiany (Russia) (Figure 1). We used feature engineering to extract latent patterns from the tremor and then systematically mined these for signi cant correlations across all volcanoes in the weeks prior to eruptions (see Methods). We identi ed characteristic peaks in a 2-day smoothed and normalized DSAR tremor (nDSAR; see Methods) that recurred with greater strength and frequency prior to eruptions at several volcanoes. By considering other data and analyses of the 2019 Whakaari eruption, we propose that this pattern indicates rapid sealing in the shallow hydrothermal system, which ultimately led to pressurization and a fatal phreatic explosion.

Correlations of pre-eruptive time-series features
In the four weeks leading up to an eruption, we observe potentially correlated patterns across at least three classes of derived tremor time series (see Methods). Prior to six major eruptions (VEI >2) at each of the studied volcanoes, there is a sustained elevation of nDSAR rate variance between 5 and 10 days prior to the event (Figure 1; see Figure S1 for equivalent time series on all eruptions). This behavior is not obvious in the raw data, but it indicates a change in the system state so that it is conducive to more voluminous sealing or pressurization with comparatively larger gas explosions.
A second archetypal pattern occurs in the nDSAR median, one week before eruptions ( Figure 2). Patterns of Whakaari-2019 and Veniaminof-2013 are especially similar. Both records exemplify a steady rise of nDSAR median to peak a few days prior to the eruption. The eruption itself often occurs after a day or two of decline. We quantify similarity in the pattern shape between pairs of eruptions by calculating the crosscorrelation coe cient, , for the 4-week records prior to eruption. For the Whakaari-2019 and Veniaminof-2013 eruption pair, = 0.71. The smoothed nDSAR pattern is broadly similar across all studied volcanoes with some variability in the timing, width and magnitude of the peak, which produces a range of values ( Figure 2a). The greatest similarity is seen between Whakaari, Veniaminof, and Ruapehu volcanoes, whereas Tongariro, Pavlof and Bezymiany, show a self-similar pattern of several decreasing cycles of the nDSAR median in the month before their eruptions.
The third pattern is observed as an increased strength of 75-min oscillations in nHF tremor, with similar timings prior to several eruptions as Whakaari 2016, 2019, Ruapehu 2009, and Bezymiany 2007a-b (See Figure S4). These eruptions were dominated by strong activity around two days prior with inverse RSAM exhibiting a linear decline ( Figure S5), indicative of cascading material failure (Voight, 1988;Chardot et al., 2015). In the next section, we show that correlation in this particular time series is spurious.

Statistical signi cance of possible eruption precursors
Demonstrating that a common pattern occurs prior to multiple events (i.e., is recurrent) is the rst step in establishing an eruption precursor. To be transferable, it should also occur prior to eruptions at multiple volcanoes. The nDSAR rate variance and nDSAR median are both recurrent and transferable (Figures 1 and 2). A third property of an eruption forecaster is that the pattern should also be rare (or absent) during non-eruptive unrest, or volcanic repose.
We tested whether values prior to eruptions were unusual during repose periods at the volcanoes (see Methods). If there is a statistically signi cant deviation, this indicates a higher degree of pattern recurrence prior to eruptive periods than expected by chance. We tested the nDSAR median pattern in the four weeks* prior to the Whakaari 2019 eruption (Figure 2b) for differentiability by convolving it with the entire Whakaari, Ruapehu, Tongariro and Veniaminof records** to generate distributions of ( Figure  3a; *see Figure S5 for equivalent analysis over three and ve weeks; **volcanoes were selected for the statistical analysis given their extended records, see Table S1). We found this pattern had high correlation prior to the four other eruptions at Whakaari, with values exceeding the distribution 90th percentile. Further, with a p-value < 0.001 (Table S3), these values are highly likely to belong to a different distribution than that of the inter-eruptive record. Eruptions at Ruapehu, Tongariro, and Veniaminof also rank highly in a percentile-sense, but the strength of dissimilarity is lower (p-values of 0.017, 0.066 and 0.072, respectively), possibly because there are fewer recorded eruptions.
A similar test was applied to the nDSAR rate variance pattern (Figure 1), again taking four weeks prior to the Whakaari 2019 eruption as the common measure of correlation (Figure 3b). Results suggest recurrency and differentiability of this pattern among the Whakaari eruptions, with all ranking above the 80th percentile (p-value of 0.006), but a much reduced degree of transferability to Ruapehu, Tongariro and Veniaminof eruptions (p-values between 0.094 and 0.149).
Finally, the 75-minute nHF harmonic prior to the Whakaari 2019 eruption ( Figure S4) was found to not be statistically signi cant (Figure 3c). Although strong activity is observed several days before six different eruptions (see Figure S4), this is not readily differentiable from other episodes that occur during repose.
Inferring eruptive processes in the Whakaari 2019 eruption Statistical analysis indicates an archetypal pattern of nDSAR median across multiple volcanoes ( Figure  2b), increasing to a peak ~2-4 days before an eruption. We use observations of the 2019 Whakaari eruption to explain this pattern in the context of the active hydrothermal system and aquifer within crater ll deposits above the vent. Speci cally, we identify transitions between ve regimes or state: (1) interaction and gas exchange between magmatic and geothermal systems, (2) pulsating gas release at the surface, (3) seal consolidation, (4) aquifer pressurization, and (5) Table S4).
A sustained RSAM harmonic signal lasting minutes to days (identi ed as "volcanic tremor"; Chouet, 1996), was detected at Whakaari between Nov 10 and 23 (regime 1, Figure 4c). This signal is considered to indicate interactions between the magma source (estimated to be 0.8-1.0 km deep at Whakaari; Jolly et al., 2018) and the groundwater (geothermal aquifer). At the surface (GeoNet, 2019a) geysering began and the crater lake level rose. On Nov 18, the Volcano Alert Level (VAL) was raised from 1 to 2, its highest noneruptive category (Potter et al., 2014).
Transition to a new regime is marked by an RSAM decrease, and concurrent increases in MF and HF after Nov 23. This was accompanied at the surface by more frequent gas emissions, water jetting and a stabilization of the lake level (GeoNet, 2019b), along with elevated SO 2 ux detected from satellite observations (Burton et al., 2021). Short, day-long cyclic oscillations are evident in both MF and HF bands (Figure 4c). Since the MF signal is stronger, the oscillations are associated with short-term rises in the unsmoothed DSAR ( Figure S7). Each cycle ends with a linear decline of inverse RSAM (Figure 4a), which suggests that cascading material failure is occurring (Voight, 1988).
Taken together, these observations imply that the Nov 23 to Dec 02 period was dominated by pulsatory gas uxing. Each cycle begins with building MF and DSAR that likely represents pressurization below a blocked conduit or weak seal. The cycle ends with cascading seal failure that renews gas ow paths, and generates surface jetting, accompanied by decrease in MF and DSAR. Sealing in the breccia-lled crater and conduit could be due to deformable clays, pore blockage by elemental sulphur, or rapid sulphate mineral precipitation along thin cracks. Seal failure and crack formation at Whakaari is likely transient (Kennedy et al., 2020). On Dec 02 there was an especially large gas pulse as seen in SO 2 ux (Burton et al., 2021), which coincides with a strong signal in the smoothed nDSAR rate variance (Figure 4a). This marked the end of the pulsatory regime.
After Dec 03, both the nMF and nHF medians (2-day normalized smoothed MF and HF tremor) decline for a few days, before a reversal and steady rise of the MF signal (regime 3, Figure 4b). Decoupling of MF and HF signals during this period is re ected in a concurrent increase of the nDSAR median (Figure 4a). There was no SO 2 ux (Burton et al., 2021) or elevated surface activity reported during this period and we hypothesize that it represents a more e cient seal consolidation and pressurization stage. This is consistent with prior association between short-term MF increase and pressurization in the pulsatory regime, but now with a larger disparity between MF and HF. A sub-surface seal will suppress uid release, and in-turn dampen surface process (decreases HF signal). Ongoing uid entry and pressurization below the seal causes stronger MF return, because HF is attenuated as signal transits from the subsurface. The mechanism of seal consolidation is uncertain but could be, e.g., determined by analysis of ballistics (cf., Kilgour et al., 2019).
The forth regime is characterized by a reversal of the nDSAR median, which declines steadily in the days before eruption (regime 4, Figure 4a). The decline mainly re ects a proportionally larger drop in MF (Figure 4b). This 'MF-quieting' likely shows that the reservoir/aquifer below the seal has reached an equilibrium pressure with the current state of deep uid recharge.
On Dec 08, ~16 hours before the eruption, a burst of activity occurred across all tremor bands (Figure 4c, waveforms in Figure S8), along with a pronounced increase in the 75-minute nHF harmonic (Figure 4a 2020) as a precursor and hypothesized to represent a sudden gas/ uid in ux. We show here, however, that pressurization and pre-conditioning of the system to eruption can be recognized much earlier.

Conclusions
Shallow-seated sudden explosive eruptions are notoriously di cult to forecast, especially in wet and noisy volcanoes. We performed a correlation analysis across time-series features derived from tremor data in the weeks preceding 18 sudden explosive eruptions across six volcanoes. We showed that there are common feature patterns, two of which can be precursors -because they are distinct from noneruptive repose. Monitoring for these observables could help to improve real-time eruption forecasting, and provide new tools and warning products, ultimately saving lives, especially at regularly visited touristic volcanoes.
Two precursors, nDSAR median and nDSAR rate variance, are especially similar across ve Whakaari eruptions, and also recognised prior to Ruapehu, Tongariro, and Veniaminof eruptions. We used the Whakaari 2019 eruption to relate these precursors with hydrothermal reservoir sealing and pressurization that primes the system for a phreatic eruption. Five distinct regimes track the onset to such shallow gas and steam-driven eruptions: (1) deep source uid input into a geothermal reservoir; (2) pulsatory gas uxing with weak temporary sealing over one or more weeks; (3) strong seal formation over the hydrothermal system several days before the event (4) leading to a critically pressured shallow aquifer (in near-equilibrium with deep gas pressures); before (5) external perturbation, and/or cascading material failure and crack formation over ~10-20 hours, leading to nal seal breach and explosive decompression.

Methods
Tremor data processing Continuous seismic waveform data were downloaded using ObsPy (Beyreuther et al., 2010) for six broadband stations close to the craters or vents of subject volcanoes (Figure 1; See Table S1). In total, 50 years of data were processed, ensuring that each eruption had at least one month of data prior (Table  S2).
The instrument response was removed from the daily waveform data and this was processed to obtain To isolate continuous signals, outlier detection is used to remove regional and volcano-tectonic (VT) earthquakes before the data streams are calculated. Events with especially long wave-trains, such as large magnitude subduction earthquakes, are removed by computing a two-window (20 minute) movingminimum of the data stream. For comparison between volcanoes, values in each data stream are transformed so they sit on a unit log-normal distribution (or a unit normal in log 10 -space). This is called zscore normalization and it has been shown to improve classi cation modeling (Singh and Singh, 2010). Thus, a normalized value of 1 indicates a datapoint at the mean, and a value of 100 indicates the signal is two standard deviations above the mean. Alternative normalization techniques based on station distance are detailed by Aki and Koyanagi (1981) and Thomson and West (2010).

Feature extraction
Each of the four data streams was subdivided into overlapping 48-hour windows, with each window advancing one data point (10 minutes) beyond the previous one. For each window, 524 timeseries features are extracted using the Python package tsfresh (Christ et al., 2018). Features included measures of distribution (parametric: mean, standard deviation; and non-parametric: quantiles, number of peaks), autocorrelation, frequency content, linearity and information content (entropy, energy, nonlinearity scores). Feature extraction yields 524 new time-series for each data stream. The effect of window length was explored by Dempsey et al. (2020) and shown not to affect resolution of short-term eruption precursors at Whakaari.
The feature referred to as nDSAR rate variance is labelled in the tsfresh python library to 'change_quantiles__f_agg_"var"__isabs_False__qh_0.6__ql_0.4'. It is calculated as the variance of paired differences (the "rate") for nDSAR values falling in the interquartile range 0.4-0.6. The feature referred to as 75-minute nHF harmonic has the corresponding tsfresh label 'fft_coe cient__coeff_38__attr_"real"' and is the real component of Fourier coe cient 38, which approximately corresponds to a frequency with period of 75 minutes.

Cross correlation analysis
For each feature time-series, e.g., nDSAR median, we analyzed four-week subsections prior to each of the 18 eruptions (See Table S2). We calculated standard correlation coe cients (Pearson method) between feature subsections from each eruption pair (see Figure 2a as reference). Once all the feature correlations are calculated (>40,000), we explored the matrices for high values between eruption pairs of eruptions and also clusters of eruptions that present strong correlations (see Figure 2a-b as reference).

Statistical signi cance tests
We conduct statistical signi cance analysis to test whether a pre-eruptive feature with high correlations is differentiable from non-eruptive repose periods. This veri es whether a result from data generated by testing or experimentation is not likely to occur randomly, but is instead attributable to a speci c cause.
To test statistical signi cance analysis, we rst select a possible 'signature', which is a four-week preeruptive feature time-series for one eruption. For example, one signature studied here is the nDSAR median prior to the 2019 Whakaari eruption (Figure 2b). The signature is then correlated with equivalent length records from each volcano (e.g., ~10 years for Whakaari), shifting forward one day at a time to calculate new correlations. This calculation is performed most e ciently as a convolution. Finally, we collect all correlation coe cients (e.g., ~2600 for the Whakaari record) and generate a histogram of their distribution and explore where the pre-eruption values fall (see Figure 2c-d; p-values are shown in Table   S3).
To quantify signi cance, we calculate p-values using a Kolmogorov-Smirnov test (KS test) (see Table   S2), which tests whether the pre-eruption correlation coe cients belong to the same distribution as the repose values. For example, a p-value for the archetype nDSAR median Whakaari-2019 (Figures 3a) correlated with the other four Whakaari eruptions generates high correlation coe cients that are very unlikely (p-value of 0.000148) to belong to the background distribution. We have tested a combination of archetype patterns and eruptions groups for different volcanoes (Table S3).
The statistical signi cance tests were performed over Whakaari, Ruapehu, Tongariro and Veniaminof record given. Records at Pavlof and Bezymiany are too short to perform conclusive statistical analysis (see Table S1).   Raw DSAR data are shown in light blue.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. SIrev20211116.docx