Drumbeat LP “Aftershocks” to a Failed Explosive Eruption at Tungurahua Volcano, Ecuador

Highly periodic, repetitive long‐period (LP) earthquakes, known as “drumbeats,” have been observed at a range of volcanoes, typically during the ascent of degassed magma. Accelerating rates of drumbeats have been reported before explosions and potentially offer forecasts of future activity. However, the broader phenomenology of drumbeats is poorly understood. Here we describe an episode of over 900 LP earthquakes recorded in November 2015 at Tungurahua Volcano, Ecuador, that we believe are associated with a failed explosion. Rates of LP drumbeats accelerated for 10 hr, consistent with an inverse Omori's law. Before any explosion occurred, seismicity decreased following Omori's law, over a further 6 days. Despite earthquake rates decelerating, amplitudes, spectral peaks, Q values, and periodicity remain constant, suggesting there is little change in the source process with time. We argue that the decelerating seismicity is a result of progressive reduction of gas flux, unable to provide sufficient overpressure for explosion.


Introduction
Active arc volcanoes of andesitic-dacitic composition are often sources of rich seismic data. Signals at these volcanoes are often dominated by long-period earthquakes (LPs), commonly associated with processes occurring in and around the magma column. Understanding these signals could be key to improving our ability to forecast volcanic activity. LPs are characterized by frequencies between 0.5 and 5.0 Hz, emergent onsets, and missing clear S wave arrivals (Chouet et al., 1994). They often begin with a mixed frequency onset, followed by low-frequency coda that decays in amplitude with time. This characteristic shape has been modeled as a two-part process with an initial excitation trigger and subsequent resonance (Chouet, 1996). These are some of the features that have been used to distinguish different categories of volcano seismic events, attributed to different source processes (Chouet & Matoza, 2013) (supporting information Figure S1). Swarms of periodic, highly similar, repeating LPs occur in a phenomenon known as drumbeats. Drumbeat seismicity is commonly associated with degassed magma ascent; however, the broader phenomenology of drumbeats is still poorly established. Locating LPs is generally a very difficult process; however, with one or two stations, careful analysis of the waveforms and their frequency content can tell us about an evolving source mechanism.
Drumbeat earthquakes are best known from the dacite spine extrusion episode at Mount St. Helens between 2004and 2005. Iverson (2008 approximated long-term steady-state behavior and slowly changing drumbeat rates and amplitudes with frictional stick-slip at the conduit margins. However, drumbeat seismicity is known to display a variety of characteristics from many arc volcanoes. Drumbeat seismicity at Soufrière Hills Here we describe a 6 day sequence of accelerating and decelerating drumbeat LP earthquakes at Tungurahua during November 2015, associated with a "failed" explosive eruption. We use a Bayesian gamma point process model (Bell et al., 2017) to examine the acceleration of seismicity rate and the subsequent decelerating rate of seismicity. We find that the drumbeats both accelerate and decelerate according to a power law with an exponent value, p = 0.96 ± 0.51 and p = 0.97 ± 0.12, respectively. Despite prolonged decaying temporal rates of seismicity, the earthquakes show strong similarity with families persisting across the 6 day sequence and amplitudes unchanging. This suggests a slowing rather than a breakdown of the driving source mechanism following a failed eruption.
First, we introduce the activity and data recorded at Tungurahua during November 2015. We then present the seismic data, along with the statistical methods for analysis. We model the data using a Bayesian point process methodology, testing different rate models and estimating parameter posterior distributions. We analyze earthquake properties including waveform similarity, families, and Q factor values. We finally present a model for accelerating and decelerating drumbeats and discuss the implications this has for magma ascent dynamics at Tungurahua.

Tungurahua
Tungurahua is a 5,032 m high andesitic stratovolcano in the Central Cordillera of the Ecuadorian Andes (Hall et al., 1999). The most recent phase of activity occurred between 1999 and 2016 with notable sub-Plinian activity in 2006 . Unrest at Tungurahua was typically associated with high rates of LP seismicity. Between the major explosive episodes of 2014 and 2016, heightened seismicity accompanied deformation and repeating tilt cycles (Bell et al., 2017;Marsden et al., 2019;Neuberg et al., 2018). This study focuses on an episode of drumbeats during one such cycle in November 2015. The drumbeats persist for 6 days and did not culminate in any explosion. There was then a repose period of 3 months before the final explosions in February 2016. Drumbeat seismicity persisted for several weeks in April 2015 and was accompanied by small explosions and ash emissions (Bell et al., 2017), while in October and early November 2015, small pulses of drumbeat seismicity emerged and ceased over just a few hours or days and are as yet unstudied (Figure 1).

Monitoring Data
The Instituto Geofísico de la Escuela Politécnica Nacional (IGEPN) maintains a volcano monitoring network on Tungurahua. The network includes short-period and broadband seismometers, differential optical absorption spectroscopy (DOAS) gas flux stations, infrasound stations, tiltmeters, GPS, cameras, and acoustic flow monitors. From a seismic network of 11 stations, IGEPN maintains a catalog of detected, classified, and, where possible, located events. Over 90% of events were recorded at RETU, a short-period seismometer at elevation over 4,000 m, approximately 2,000 m from the crater rim. This proximity means the signal-to-noise ratio (SNR) is high and many small, shallow events are recorded. We manually picked 932 events from 25 to 30 November 2015 for this study, representing all detectable events at RETU. These events were only visible at the one station, and with emergent onsets and no clear S phases, locating the events was not possible. As the seismicity is only recorded at this uppermost station, we believe these LPs are associated with shallow processes in the top 2,000 m of the conduit . The similarity of the waveforms indicates that they are all closely colocated within a small depth range. Given this colocation, we use the maximum amplitudes of individual events as a relative comparison for magnitude. Six hundred fifteen of the manually picked events appear in the IGEPN catalog. However, there are only 20 events which are located and have estimated magnitudes, all of which are less than magnitude 1.5 and carry large uncertainties. The seismicity on 25 November is the first clearly identifiable sequence of LP events, as in the preceding days the signal at RETU is dominated by emission tremor.
Details of surface observations and ash column heights are extracted from daily reports produced by the Observatorio del Volcán Tungurahua (OVT) (https://www.igepn.edu.ec/) and used in conjunction with the seismic data for temporal analysis. Explosion counts and radial tilt measurements at station RETU are also collected from IGEPN catalogs.

Methods
The seismic data are initially processed using the ObsPy toolkit (Krischer et al., 2015). Thirty-second duration waveforms are sliced and bandpass filtered between 1 and 40 Hz. The maximum amplitude of each event is extracted. Fast Fourier transform (FFT) of each signal is calculated to generate a periodogram. We find the power spectral density (PSD) for frequencies sampled at an interval of 0.01 Hz and extract the maximum value as the fundamental peak frequency.
The Q factor for each event is calculated using an autoregressive moving average (ARMA) technique, adapted from Seismo-Volcanalysis software (Lesage, 2007). The Q factor is a nondimensional number that describes how quickly or slowly wave energy dissipates and is often strongly linked to the fundamental peak frequency. Autoregressive methods have been successfully used to analyze changing LP frequency contents (Kumagai & Chouet, 1999;Lokmer et al., 2008). The approach is similar to the commonly used Sompi method (Hori et al., 1989). A signal is composed of a number of individual harmonic decaying oscillations. Each component can be represented in complex frequency space and quantified by a peak frequency (f , Hz) and growth rate (g, s −1 ) (Kumazawa et al., 1990). We generated cumulative f-g diagrams for all filters between 2 and 30, and points that cluster around a pole at the spectral peak are used to calculate Q (Equation 1) (Cusano et al., 2008). In this automated adaptation of the ARMA method, a hierarchical clustering method is used to automatically select points in the complex frequency space (Eads, 2008).
We determine the maximum cross-correlation coefficient between 0 and 1 for all pairs of events in our catalog and use a threshold value to group events into families (Park et al., 2019;Waite et al., 2008;Yukutake et al., 2017). Following previous studies, including that of LP drumbeat seismicity in April 2015 at Tungurahua, the threshold is set at 0.7 (Bell et al., 2017;Petersen, 2007;Thelen et al., 2011).
We also calculate earthquake interevent times (IETs) and their periodicity to highlight times of pronounced drumbeat activity. Bell et al. (2017) defines periodicity as the ratio between the mean, , and standard deviation, , of IETs. Events randomly distributed in time, with an average rate, , will have a probability density function such that = 1. Clustered events have periodicities less than 1, whereas periodic events have periodicities greater than 1.
Finally, we considered models for the accelerating and decelerating components of the drumbeat episode. Previous studies of accelerating seismicity have modeled rates using power, exponential, and hyperbolic relationships Ignatieva et al., 2018). In accelerating and decelerating components we opted to model the event rates using an exponential relation (Equation 2) and a power law. For the decelerating event rates this is the modified Omori's law (Equation 3), and for the accelerating rates, an inverse Omori's law (Equation 4).
We model the drumbeat sequence as an inhomogeneous Gamma process . We define the point of maximum seismicity as, t 0 , separating the accelerating and decelerating components. We use a Bayesian approach with PyMC3 implementation (Salvatier et al., 2016). We use Markov chain Monte Carlo (MCMC) to sample the posterior distributions of model parameters. We run 5,000 iterations. We provide initial estimates for parameters p, t f , k, and . The prior distribution and rate parameters used are detailed in Table S1.

Results
Across the 6 day period from 25 to 30 November we see an initial increase in the rate of seismicity before a rapid deceleration (Figure 3). The peak in event rate occurs at 10:00 on 25 November (Figure 2). During the drumbeat episode the radial tilt increases and decreases in a range of 10 μrad.

Drumbeat Onset
The first 10 hr of the drumbeat sequence is markedly different from the activity observed thereafter. In the first 10 hr, the event rate increases, the individual event amplitudes increase slightly, and the seismicity becomes increasingly periodic (Figure 2). The point process modeling shows the accelerating rates of seismicity can be defined by a power law (Figure 3a). The best fitting exponent is p = 0.96 ± 0.51. The best fit value for t f is 0.53 ± 0.16. The t f value gives a "hindcast," indicating when an eruption may have occurred, if drumbeats continued to accelerate. The point process model finds t f and the possible eruption at approximately 12:43 p.m. This is nearly 3 hr after the peak rate of seismicity and onset of the deceleration.
Through this sequence, the Q factor shows little to no change. A few events before 10:00 on 25 November are included in families that persist through the episode, but there are no families that exclusively represent events prior to 10:00. This suggests that accelerating events before 10:00 are not generated by a distinct source that differs to that active later in the episode. There is no specific activity from the surface observations that can be pinned to the change in seismicity at 10:00. From 18:00 on 24 November to 18:00 on 25 November there were observations of energetic gas and steam emissions between 300 and 800 m. Notably, this was with no ash content, but for the majority of the day the summit remained obscured by cloud (IGEPN, 2015).

Drumbeat "Aftershocks"
During the afternoon of 25 November the earthquake rate decreases rapidly, while the periodicity remains greater than 1. Event amplitudes and the Q factor show no systematic change with time ( Figure 2). The average cross-correlation coefficient also remains stable through the sequence but begins to decrease after 12:00 on 29 November, and fewer events after this point are grouped into families. Unlike previous studies of event families where systematic shifts in families correspond to physical changes in activity (Bell et al., 2017;Green & Neuberg, 2006;Thelen et al., 2011), there is no chronological variation in the families identified here. The largest family contains 177 events, representing 20% of all picked events in the sequence, and persists across all 6 days, of both accelerating and decelerating seismicity. There is strong waveform similarity between families, as the master events also show high cross-correlation coefficients with one another ( Figure S2).
The decrease in seismicity rate is well modeled by a power law. The best fit parameters are calculated where p = 0.97 f ± 0.12 (Figure 3). The power law exponent is consistent with the Omori's law relationship where p = 1, and so the rate of earthquakes is inversely proportional to the period of time since the peak of activity. Previous studies examining tectonic aftershock sequences using modified Omori's law have found p values between 0.6 and 1.8 (Holschneider et al., 2012;Wiemer & Katsumata, 1999). Tectonic aftershock sequences are typically described as clustered, whereas seismicity here is quasiperiodic. The exponential model provided a less successful fit ( Figure S3 and Table S2).

Discussion
The surface observations, tilt data, and seismicity analyses allow us to build a conceptual model for the process generating this seismicity. The point process modeling reveals two distinct phases of accelerating and decelerating seismicity. The accelerating drumbeats follow a power law increase, according to the inverse Omori's law (Equation 4) with exponent p = 0.96 ± 0.51 (Table S2). Hindcast values for t f suggest that had seismicity continued to accelerate this may have culminated in an explosion nearly 3 hr later. This is comparable with the accelerating LP drumbeats prior to a large explosion in July 2013 at Tungurahua where the exponent is p = 1.05 . Five-minute real-time seismic amplitude (RSAM) was calculated for 00:00-10:00 on 25 November for comparison. Unlike July 2013, amplitudes of individual events in this study do not accelerate to the same degree, but through the majority of the drumbeat episode, they do generate RSAM values in the same range (∼50-150) ( Figure S4) . The decelerating phase closely mirrors the acceleration, as point process modeling fits the modified Omori's law (Equation 3) with an exponent, p = 0.97 ± 0.12.
Drumbeats that show short IETs and similar waveforms require a single, rapidly recharging, nondestructive source process (Neuberg et al., 2000). Families of events persist through the phases of accelerating and decelerating rates of seismicity. High periodicity and unchanging Q factors and amplitudes also persist throughout, suggesting the seismicity is likely generated by the same source location and process through the full episode.
The Q factor of individual waveforms can be used to determine the composition of resonating fluids generating LPs. Q values in the literature can vary from 3 to 300 (Molina et al., 2004). Although the Q values in November 2015 do not show a clear increase or decrease, they are consistently in the range 3-15, which fit well with a high gas fraction CO 2 or SO 2 water (Kumagai & Chouet, 2000). This is unlike the dusty ash-gas mixture examined in Molina et al. (2004) where Q systematically increases from 200 to 500. This also corresponds well with the daily records that on 25 November steam and gas emissions continue but notably without any ash.
A common component of models for LP drumbeats is a process of loading and release. However, the controlling mechanisms proposed have included stick slip processes (Iverson, 2008), magma failure (Green & Neuberg, 2006;Tuffen et al., 2003), and gas escape and depressurization (Bell et al., 2017). A plug model is illustrated in Bell et al. (2017), and plug sealing and destruction cycles during 2013-2014 at Tungurahua are well documented . Although these drumbeats follow a period of small explosions, there is no evidence that a plug was destroyed following April 2015 drumbeats (Bell et al., 2017). Instead, we consider the plug as an intermittently permeable barrier with fracture networks and transient pathways to allow gas, ash, and magma to ascend through it. The driving process is ascending magma in the conduit which is degassing, and there is a decelerating force that leads to longer loading and failure cycles at a plug in the upper conduit. Irrespective of the physics of the excitation mechanism in this episode, we can put constraints on the nature of the loading process.
We propose a model for the timeline of events in November 2015 (Figure 4). During 11-17 November initial magma ascent in the conduit, as detected in the inflation signal at RETU, caused depressurization and fragmentation leading to minor explosions and ash emissions (Cassidy et al., 2018). In the plug, there were persistent pathways for gas and ash escape, without increasing overpressure in the conduit. Observations of incandescent blocks on the flanks suggests magma was extruded. We assume this initial ascent was an aseismic process Salvage & Neuberg, 2016). Following the explosions, there was a further period of repose in surface activity (IGEPN, 2015), but the radial tilt still increases, indicative of ongoing magma ascent ( Figure S5). By 25 November, the radial tilt increase slows, and we observe increasing rates of drumbeats. This timing in the tilt cycle is the same as that of April 2015, though different to most in which SO 2 flux and seismicity peaks occur during minima in tilt inflation cycles (Bell et al., 2017). As the tilt cycles and episodes of seismic unrest are out of synchronization with one another, this suggests there was something unusual about the episodes of drumbeats in April and November 2015. Outgassing was ongoing and building pressure beneath the plug as the rate of release exceeds the rate at which it can freely escape. Each time the failure strength was surpassed, escaping gas interacts with this complex fracture network and generated LP events. Between 00:00 and 10:00 on 25 November the accelerating rate of LP drumbeats suggests there was an increase in the gas accumulation and hence the loading rate. This overpressure, however, was insufficient for an explosion. Seismicity persisted as outgassing continued but at a decreasing rate, and there was a decreased loading rate. The time to reach a failure strength increased, and so the repose period between events remained cyclic but increasing. The loading rate and seismicity decreased as a power law relationship until a state of equilibrium and repose. As the body was not ascending at this stage or under any further external pressures, we assume the plug was in a steady state. We suggest degassing continued at a reduced rate, but we cannot be certain for how long.
The next paroxysm in February 2016 included some of the most violent activity in the whole eruption, with pyroclastic flows and over 130 explosions in the first week of March (IGEPN, 2016). It is therefore significant that we see heightened seismicity and acceleration up to a failed explosion, before several months of quiescence. This provides valuable insight into forecasting potential explosions using seismicity at intermediate composition volcanoes. The power law fit during the accelerating seismicity indicated an eruption time at 12:43 on 25 November, by which time the seismicity rate had already decreased 30% from the peak rate.
Subsequent deformation analysis has used InSAR to examine inflation on the western flank during the first 3 weeks of November 2015, prior to the onset of drumbeat LP seismicity. The most likely source model identified was an inclined body of shallow, short-term, preeruptive magma storage (Hickey et al., 2020). This shallow magma storage may have contributed to processes in the conduit prior to drumbeat LPs.

Conclusions
The analysis of the November 2015 drumbeats contributes to the narrative of the final years of activity in the 1999-2016 eruption at Tungurahua. Although the sequence overall shows similarities with the drumbeats observed in July 2013 and April 2015, the driving mechanism and subsequent surface activity is likely different. Across accelerating and decelerating phases, the earthquake signals show minimal changes in their properties including Q factor and amplitudes, suggesting the source process is unchanging. Q values are consistent with LP excitation by a gas-rich mixture as opposed to magma or ash-gas mixtures.
We suggest that the observation can be explained by a model involving ascent of a degassing magma body.
We initially see acceleration of drumbeat LP events and an increasing rate of loading on the plug. Insufficient overpressure from gas flux causes an increasing repose time in loading-failure cycles at the plug and hence a decreasing rate of seismicity. We believe this is a "failed" explosive event, followed by several days of drumbeat LPs where magma in the conduit continues to degas.
This work contributes to further understanding the broader phenomenology of drumbeat seismicity and their source processes. This last instance of degassing is key to bringing the conduit into a state of equilibrium before several months of quiescence. Future work will look to extend these analyses to quasiperiodic seismicity throughout 2015 and into the final explosive eruptions during 2016. For volcanic hazard analysis, it is important to highlight that high rates of seismicity are not necessarily a precursor to substantial explosions. Careful analysis of the seismicity rate and individual event properties, alongside the surface observations, show these drumbeats are part of a complex process.