Evolution of Seismicity During a Stalled Episode of Reawakening at Cayambe Volcano, Ecuador

Cayambe Volcano is an ice-capped, 5,790 m high, andesitic-dacitic volcanic complex, located on the equator in the Eastern Cordillera of the Ecuadorian Andes. An eruption at Cayambe would pose considerable hazards to surrounding communities and a nationally significant agricultural industry. Although the only historically documented eruption was in 1785, it remains persistently restless and long-period (LP) seismicity has been consistently observed at the volcano for over 10 years. However, the sparse monitoring network, and complex interactions between the magmatic, hydrothermal, glacial, and tectonic systems, make unrest at Cayambe challenging to interpret. In June 2016 a seismic “crisis” began at Cayambe, as rates of high frequency volcano-tectonic (VT) earthquakes increased to hundreds of events per day, leading to speculation about the possibility of a forthcoming eruption. The crisis began 2 months after the M w 7.8 Pedernales earthquake, which occurred on the coast, 200 km from Cayambe. Here we show that the 2016 seismicity at Cayambe resulted from four distinct source processes. Cross correlation, template matching, and spectral analysis isolate two source regions for VT earthquakes–tectonic events from a regional fault system and more varied VTs from beneath the volcanic cone. The temporal evolution of the LP seismicity, and mean Q value of 9.9, indicate that these events are most likely generated by flow of hydrothermal fluids. These observations are consistent with a model where a new pulse of magma ascent initially stresses regional tectonic faults, and subsequently drives elevated VT seismicity in the edifice. We draw comparisons from models of volcano-tectonic interactions, and speculate that static stress changes from the Pedernales earthquake put Cayambe volcano in an area of dilation, providing a mechanism for magma ascent. Our findings provide a better understanding of “background” seismicity at Cayambe allowing faster characterization of future crises, and a benchmark to measure changes driven by rapid glacial retreat.


Motivation
Seismic records at active volcanoes can be a rich and informative dataset for monitoring and understanding internal processes. The characteristics of a seismic signal can be indicative of the source process that generated it (Zobin, 2012;Chouet and Matoza, 2013;McNutt and Roman, 2015). If recorded across multiple stations then the locations of the earthquakes can be used to infer the source mechanism (Woods et al., 2018;Lehr et al., 2019). The rates at which earthquakes occur can also indicate whether a source mechanism is constant and stable or evolving. Accelerating sequences of earthquakes have been linked to material failure and can be significant in forecasting eruption (Voight, 1988;Boué et al., 2015;Bell et al., 2018). Whilst very long term, constant and sometimes periodic earthquake rates may suggest a stable or resonating process (Park et al., 2019). The emergence of new swarm seismicity including spasmodic bursts may suggest a change in the stress regime, and may be associated with fluid migration, or rock fracturing (Hill et al., 2002a;Hotovec-Ellis et al., 2018). Understanding and characterizing different types of volcano-seismic events allows us to attribute them to specific processes, and can be particularly important in hazard assessment and forecasting. This is particularly the case at volcanoes which have been dormant. When there is not an extensive historical record of seismicity at a volcano, recognising patterns may help to identify signs of "reawakening" .
Interpreting patterns of seismicity in terms of volcanic processes and evolving unrest can be challenging. Earthquakes in volcanic regions can originate from a number of different sources, and may be independent from processes that may drive new eruptive activity. Active fault systems, hydrothermal systems and permanent glacial cover can add another element of complexity to geophysical monitoring at volcanoes (Jónsdóttir et al., 2007;Cusano et al., 2008). Particularly within the tropics, small, high altitude glaciers can respond to small changes in climate on short timescales, and may be susceptible to significant impacts from El Niño and La Niña years (Ceballos et al., 2006;Basantes-Serrano et al., 2016). Whilst it is known that equatorial glacier retreat has been ongoing for 150 years, this has accelerated in the last 30 years and should be closely monitored at active volcanoes (Manciati et al., 2014). Glaciers on a volcano can act as a constant supply of fluid to alter and weaken the ediface (Carrasco-Núñez et al., 1993) and so increasing retreat may lead to an acceleration of these processes. In the event of a volcanic eruption, glacier melt waters can contribute to disastrous debris flows, such as those seen at Nevado del Ruiz in 1985 (Huggel et al., 2007) and Nevado del Huila in 2007 (Worni et al., 2012). Several large non-cohesive lahars involving glacial melt, have been mapped in the valleys surrounding Mt Rainier, where runouts extend over 100 km (Hoblitt et al., 1998).

Cayambe Volcano
This study focuses on Cayambe Volcano, in the north of the Eastern Cordillera of the Ecuadorian Andes, just 60 km north of the capital Quito. The complex is a notably large volcanic centre, with a peak elevation of 5,790 m. Above 4,800 m the volcano is covered in a 22 km 2 ice cap, which in places is as thick as 100 m (Guillier and Chatelain, 2006;Detienne et al., 2017). Tephrochronology studies have revealed previous explosive activity in the last 4,000 years, including pyroclastic flows, lava dome extrusions and collapse events (Samaniego et al., 1998). Hydrothermal alteration of flank deposits and thawing soils from glacial retreat act to destabilise the western slopes of Cayambe. Small perturbations in the system, for example, from a magmatic intrusion or seismic event, could lead to catastrophic failure and significant landslides (Hoblitt et al., 1998;Detienne et al., 2017). This would pose a threat to the town of Cayambe, population 20,000, that sits just 15 km west of the volcano (Samaniego et al., 1998). Mass debris flows and lahars descending the eastern flanks of Cayambe could also easily damage the Trans-Ecuadorian oil pipeline and major highways to the oilfields in northeastern Ecuador (Schuster et al., 1996). Whilst there has only been one historically documented eruption at Cayambe (1785-1786), it has been persistently restless in the last 20 years and mountaineers regularly report sulphur smells from the upper flanks. The 1785 eruption is thought to have been a Vulcanian type eruption ∼VEI 1-2 (Global Volcanism Program, 2013;Bernard and Samaniego, 2017), similar to the Cotopaxi eruption of 2015 (Bernard et al., 2016) and those seen at Tungurahua in 2014 (Hall et al., 2015).
Activity at Cayambe is monitored by the Instituto Geofísico de la Escuela Politécnica Nacional (IGEPN), the national geophysical monitoring agency of Ecuador. In 2016, IGEPN detected a seismic crisis at Cayambe. Swarms of highfrequency volcano-tectonic (VT) events in June and September were recorded by three seismometers on the flanks of Cayambe volcano. These swarms involved daily earthquake rates significantly above previously recorded baseline seismicity. These swarms included two events of magnitude 3.3 and 3.7 during November, the largest two VT events of the year at Cayambe (Global Volcanism Program, 2013;IGEPN, 2016), and raised speculation about the possibility of an imminent eruption. In this study we characterize the seismicity during 2016 at Cayambe and try and to quantify the differences between the events of the seismic crisis and typical baseline seismicity in the preceding years. Subtle changes in the character of the waveforms can be indicative of a changing source or path process. Previously, detailed retrospective studies of temporal and spatial patterns in seismicity, have revealed key characteristics in eruption styles at comparable andesitic volcanoes, including Volcán de Colima (Zobin et al., 2015), Redoubt Volcano (Buurman et al., 2013) and Soufrière Hills Volcano (Neuberg et al., 2006).
Also of note, on April 16, 2016, 2 months before the onset of seismic unrest at Cayambe, the M w 7.8 Pedernales thrust earthquake struck the coast of Ecuador, 200 km from Cayambe (Heidarzadeh et al., 2017;Nocquet et al., 2017;Yoshimoto et al., 2017). This event caused MMI intensity III to IV ground shaking around Cayambe, and was followed by seven M w > 6.0 events in the following 3 months (USGS, 2021). Tectonic earthquakes can generate both permanent static changes and temporary dynamic changes in the stress field (Manga and Brodsky, 2006). Static changes are a direct result of slip on the earthquake fault and crustal deformation and are more prevalent closer to the earthquake epicentre. Whilst dynamic stress field changes are related to passing teleseismic surface waves and can effect distances much further from the earthquake (Prejean et al., 2004). Stress changes from tectonic events are known to be able to trigger both eruptive activity and seismic activity at volcanoes (Hill and Prejean, 2015). However, distinguishing between causation and correlation is complex, particularly if there are limited records of historical earthquakes or eruptions (Prejean and Hill, 2018). Here we investigate the relationship between the Pedernales earthquake in April 2016 and unrest at Cayambe that commenced in June 2016.
The permanent ice cap on Cayambe presents a further potential complication for hazard assessments from seismic records. Studies further south in the Chilean Andes have shown icequakes to be much more common than previously thought (Lamb et al., 2020). These can be generated by basal slip on the glacier, slip at fractures within the glacier or the resonance of fluids within glacier cracks. Icequakes observed at Cotopaxi volcano, Ecuador and Katla volcano, Iceland, show waveforms that look very similar to long period (LP) earthquakes (Métaxian, 2003;Jónsdóttir et al., 2009). We refer to LP events as low frequency emergent signals, often associated with fluid movement, including magma emplacement and hydrothermal systems Woods et al., 2018). Therefore being able to clearly distinguish between the two phenomena is key to reducing uncertainty in hazard assessment. This is particularly the case when we cannot confidently constrain depths and locations to separate sources.
Here we introduce the seismic data records during the 2016 crisis. We then describe the template matching process used to expand the catalogue of picked events. We detail the different methods used to quantify sequences of significant events throughout the year. We then present an interpretation of the seismic swarms that occurred during 2016. Our results identify a stable source of long period seismicity, likely associated with a shallow hydrothermal system. We also demonstrate the two VT swarms to be separate populations, from both local tectonic and volcanic sources. We draw comparisons between the seismic observations and models of re-awakening volcanoes, and speculate the significance of the Pedernales event prior to activity at Cayambe in 2016. Finally we discuss the wider implications of a seismic crisis like this for hazard assessment in the future.

Monitoring Data
The volcanic complex at Cayambe is defined by three cones, and bound by a significant regional fault to the north east, the "La Sofía-Río Chingual" fault. Old Cayambe is a basic-andesitic stratovolcano that is more than 250,000 years old. It has mostly been eroded and the more recent Nevado Cayambe, a composite cone has been built in its place, and over the eastern flanks (Samaniego et al., 1998). There has been a permanent geophysical monitoring network on Cayambe volcano since 1988, maintained by IGEPN (Guillier and Chatelain, 2006). Initially a single station on the northern flank, the network is now comprised of three stations, with telemetry to IGEPN in Quito. The seismic network includes one, threecomponent broadband seismometer (Trillium Compact 120), CAYA, and two short period seismometers (L4C), CAYR and ANGU, each sampling at 100 Hz ( Figure 1). There was also a short temporary deployment of seismometers between 1997 and 1998 (Guillier et al., 1999;Guillier and Chatelain, 2006). Station CAYA is positioned between the Nevado Cayambe cone and the La Sofía-Río Chingual fault, whilst ANGU and CAYR are positioned close to the edge of the glacier line on the western flanks. For this study, we have access to 366 days of data from station CAYA in 2016, and data from January 2015 through to December 2017 at station CAYR to specifically examine the LP seismicity. Waveform data from station ANGU are unavailable during July-September 2016 and the station was upgraded from a short-period to a broadband instrument in December 2016. As such we are only able to use data from January to November in 2016.

Signal Processing
Seismic data are analysed and processed using the ObsPy toolkit (Krischer et al., 2015). Our initial analysis uses a catalogue of manually picked and categorized events, provided by IGEPN. Events in the catalogue are labeled broadly as either LP, VT, hybrid (HYB) or tectonic (TECT) events. We remove tectonic events, and cross reference the Cayambe catalogues, with catalogue data for nearby Reventador volcano and a tectonic events catalogue to ensure there are no falsely identified events. We filter the data between 1 and 40 Hz. After earthquakes are identified, we create a 30 s event waveform for analysis, slicing 5 s before the pick time and 25 s afterward. Station CAYR lies much closer to the summit of Nevado Cayambe (∼2.5 km) and exhibits a much higher signal to noise ratio (SNR). As a result, over 80% of picked events in the catalogue are identified at CAYR (Figure 2). The IGEPN catalogue also includes some located VT events. These are used to map and interpret some of the significant swarms in 2016.
First, using catalogue picks at station CAYR, we run a cross correlation analysis. Each event is cross correlated with every other event, generating a cross correlation coefficient (XCC) between 0 and 1 as well as and the shift, in seconds, required to align the two picks at the maximum coefficient. The coefficient can be considered a measure of similarity (Petersen, 2007;Waite et al., 2008;Thelen et al., 2011). We set a minimum coefficient threshold of 0.7 to group events into families (Bell et al., 2017;Park et al., 2019;Butcher et al., 2020). The event with the most coefficients greater than the threshold is selected as the master event of the first family and every event with a coefficient greater than the threshold is included in the family. This is repeated until the families are too small to be considered significant. We set a minimum requirement of 20 events to be considered a family. The master event from each family is then used in a template matching method to extend the catalogue (Gibbons et al., 2007;Park et al., 2019). The master event is cross correlated with the signal for each whole day, and when a cross correlation threshold is exceeded, an event is picked. We then consider patterns within and across the families identified from the cross correlation analysis. The master event of each family is cross correlated with every other master to determine which families might be closely related. Every event in a family is shifted in time to align with the point of maximum cross correlation and stacked. These stacked events are then also cross correlated.
For each event, we determine the Q factor of the signal. The Q factor has been used in volcano-seismology studies to quantify and group event types (Molina et al., 2004;Lokmer et al., 2008). Q describes how quickly energy in a signal is attenuated, either as a result of a source or path process, and can be used as a metric to identify patterns of change in series of earthquakes (Rossing and Fletcher, 2004). Q and the interpreted attenuation is often empirically related to resonance in a fluid filled crack model  (Kumagai and Chouet, 1999;Lipovsky and Dunham, 2015). This value can be calculated by a number of different methods and reported values for Q vary between 10 0 and 10 3 (Molina et al., 2004;Lipovsky and Dunham, 2015). This range is validated by a series of synthetic tests for fluid-filled cracks where the fluid compositions are varied (Kumagai and Chouet, 2000). In any one study, however, only one method is used to calculate the Q value for an event. These broad ranges in Q values can in turn lead to very different interpretations of the seismic event sequences. We apply three different approaches to calculate Q for each event and compare the output. We will use these Q values as a tool to distinguish earthquake events associated with different processes.
The decay method is used in a comparative study across Mt Etna, Colima and Campi Flegrei volcanoes (Elmore and Heald, 1985;Del Pezzo et al., 2013). The method relies on decay time, τ, defined as the time taken for the squared amplitude of the waveform to decay to less than 1 e of the squared maximum amplitude (Eq. 1). This calculates the Q decay value for the oscillations of a single frequency, and so assumes the waveform spectra is dominated by a single peak frequency, f 0 (Rossing and Fletcher, 2004). As such, this method works consistently for LP type events, and specifically single frequency tornillos or on infrasound data, but is less robust with broadband signals.
We also apply the bandwidth method, which is closely related to the decay method. It relies on identifying the bandwidth of the spectral peak frequency, δf . The bandwidth is the difference between the high and low corner frequencies of the power spectral density (PSD). Q band for this spectral peak is then defined as the ratio between the spectral peak and it's bandwidth (Elmore and Heald, 1985) (Eq. 2). This method is presented using smoothed spectra for the analysis of tornillo infrasound events at Cotopaxi Volcano in 2015 and 2016 (Johnson et al., 2018). In this study the bandwidth, δf is calculated for the single largest spectral peak. We do not apply smoothing to avoid imprinting any interpretation on the spectra.
Finally, we apply an auto-regressive moving average (ARMA) methodology to determine Q ARMA , as described in Lesage (2009). The ARMA approach assumes a signal is a superposition of many individual harmonic decaying oscillations, each defined by their own frequency, (f, Hz) and decay rate (g, s −1 ). We consider all AR filters from 2 to 50 and plot cumulative f-g diagrams (Lesage, 2008). The method is automated to identify points in the complex frequency space that cluster at a pole around the spectral peak (Butcher et al., 2020). This method is closely related to the Sompi method which uses eigen-decomposition to identify poles and calculate Q ARMA (Hori et al., 1989;Kumazawa et al., 1990) For each event, we also identify the peak spectral frequency from the periodogram of the signal, and the maximum absolute amplitude of every event. On single station records, we can then consider maximum amplitude as a proxy for magnitude and gather information about amplitude ratios between stations.
In order to identify clustered swarms of events, and locate any episodes of periodic activity we calculate the inter-event times (IETs) and periodicity (Bell et al., 2017). This is done for the whole catalogue initially and then within families and by event type. The periodicity is the ratio between the mean and standard deviation of the inter-event times. This is a useful tool to identify swarms and periodic episodes as it assumes events are randomly distributed and so have a ratio and periodicity equal to one. Clusters and swarms will have a periodicity much lower than one and periodic episodes of drumbeat-like, repeating seismicity will have periodicities much greater than one.

Located Events
Due to the sparse monitoring network on Cayambe, only 179 of the 11,990 VT events in the IGEPN catalogue are located, and even then with significant uncertainties. Despite this, the locations identify at least two distinct VT clusters throughout the year ( Figure 3). Manual inspection of the records show the events group broadly in four very distinct time intervals, with little overlap: June, July, September and November-December. We use these locations to aid our interpretations from the individual event analysis. The June events are widely spread over the north and eastern flanks, whilst the November events are much more concentrated near upper reaches of the Nevado Cayambe cone. The network geometry is such that the shallowest events are best captured by short period station CAYR on the highest flanks of the Nevado Cayambe, and is the most reliable station in the network.

Catalogue Overview
From the total 11,990 events in the catalogue in 2016, 9,848 events are recorded at station CAYR, of which 46% are labeled as LPs, 48% as VTs, and the remaining 6% as hybrids (HYB) and tectonic events (TECT). This is in stark contrast to 2015 where 98% of events recorded were labeled as LP, although there are significant data gaps. The cross correlation analysis generated 31 unique families of greater than 20 events ( Figure 4). Despite the catalogue containing almost an equal split of LP and VT events, only 84 LP events were grouped in two families representing more long term, background LP seismicity spanning the whole year. The remaining 29 families represent 1,545 VTs-seven of which highlight a short VT swarm in June and 21 families span more long term VT seismicity from September to December.
We use the master events from the 31 families to conduct template matching and identify any events not originally picked in the IGEPN catalogue. An additional 3,853 events are identified by this template matching method, extending the catalogue to 13,683 events in 2016. In particular the template matching identifies 364 new LP events-more than four times as many LPs that were in the original two LP families. The template matching identifies a further 1378 VT events in the June swarm and 2111 VT type events from September to December 2016. The template matching expanded the reach of some VT swarm families but no single family contains events in both June and November, showing a strong distinction between the source processes in each case. We also use the located VT events in the catalogue to verify the master events of the families. The located events in June and July correspond well with families 0, 6, 9, 10, and 23 (Supplementary Figure S1A). Whilst the master event of family 20 shows an XCC > 0.85 and is just 20 min prior to a located VT event (Supplementary Figure S1C). This suggests family 20 represents an isolated repeating source mechanism between the two large swarms in June and November. Figure 5 shows seismic event characteristics through 2016. All event types have a near exponential distribution of inter-event times (t) (Bell et al., 2017). The periodicity is calculated in fixed windows of 15 events, and so during swarm activity in June and December there are isolated incidents where seismicity is either clustered (periodicity <1) or cyclical (periodicity >1), but at generally low levels.
In general the LP events show a consistent peak frequency at 2.1 Hz. As expected this is lower than the VT events, but particularly in October-December when the VT events have a consistent peak frequency at 5 Hz. The June swarm VT events have uncharacteristically low spectral peak frequencies, less than 2 Hz and less than the background LPs. Thorough inspection of the June events show that they have typical spectral contents expected for VTs, between 1.0 and 15.0 Hz. This highlights how using a single peak frequency value as a defining characteristic can be misleading.
The Q values by all three methods return values across multiple orders of magnitude. A comparison for all three methods on the same family is shown in Figure 6. The Q ARMA values are calculated by an automated clustering algorithm and so the standard error from the clustered poles is used to plot the error bars in Figures 5E, 6F,L. This output from the three separate methods can lead to spurious interpretations and conclusions.

Persistent LP Seismicity
The catalogue shows a consistent level of baseline LP seismicity in 2015-2017. In 2016, LPs make up 46% of identified events at CAYR. However, from the cross correlation and family analysis, only two LP families emerge (families 7 and 26, Figure 4). The template matching expands the catalogue of LP picks (Figures Unlike some of the higher frequency events, these LP events are poorly recorded on stations ANGU and CAYA. In particular the high frequency response from ANGU impedes how well the LP waveform and spectra is recorded ( Figure 8C). Given the sensitivity of the cross correlation method, we expand the template matching search for LPs into 2015 and 2017 but only with station CAYR. Even beyond 2016, metrics such as the spectral peak frequency and Q ARMA show very little variation. The temporal analysis shows these LP families have a periodicity average at 1.1, suggesting an exponential inter-event time (t). Short episodes of clustered or cyclical seismicity are generally artifacts as a result of the chosen window length. The event rate is low compared to VT swarm event counts, just 1-2 events per day.
In addition to the repeating LP seismicity classified by cross correlation, some less consistent low-frequency events were identified in the data. These events are much longer (1-5 min in length) and have highly harmonic frequencies ( Figures  7D-G). Due to the varying lengths of these events, they are not classified into families successfully. They are identified throughout the year, although more are located in October and November. Single frequency events with long coda and distinct Hilbert transform envelopes are often referred to as tornillos (Narváez et al., 1997;Fazio et al., 2019). Some of the identified events match this description; however, many show interesting characteristics such strong integer harmonics ( Figures  7D,E). Sequences of these tornillo-type events can appear to merge and generate long episodes, which may be considered harmonic tremor (Figures 7F,G). Particularly in Figures 7E,G, when fundamental peak frequency increases through time, this may be an example of a short burst of gliding tremor (Almendros et al., 2014).

June 2016 Distal Swarm
The cross correlation analysis picks out seven distinct families in June 2016. The largest family found contains 229 events, all occurring between 8th and 17th June (Figures 6A-F). However, applying the template matching method identifies an additional 346 events and extends the duration of the family to 26th July. This family is characterized by short (<15 s), high frequency, VT events. They have a distinct signal envelope, clearly showing Pand S-phase arrivals. Three spectral bands emerge around 1.6, 4.5, and 7.9 Hz and are the peak frequency in 44% of events in the seven families which represent seismicity in June (Figures 6B, 9). When the master and stacked events of each family in the June swarm are cross correlated with one another, they show high correlation values (Figure 10A), suggesting a highly repeatable FIGURE 4 | Families extracted from cross correlation analysis of the initial 9,484 catalogue picks at CAYR. Individual events colored by their label, red shows LPs, blue shows VTs. Families numbered in order of size largest to smallest.
Frontiers in Earth Science | www.frontiersin.org June 2021 | Volume 9 | Article 680865 source and path. Q ARMA shows very small errors compared to some other swarms and families ( Figures 6L, 10). The June swarm events are clear on both CAYR and ANGU, although ANGU shows a particularly high frequency response. These events are harder to see at CAYA but, when stacked, do show an improved signal to noise ratio ( Figure 8B).

September 2016-Proximal Swarms
From September onwards, the families extracted from the cross correlation method persist longer than those in June. Some VT swarm families even contain a small number of events identified as LPs by the IGEPN, although in any family these amount to no more than 8% of the total events. These LP events were manually checked and verified to ensure the labeling is correct. These families represent VT events which are more LP-like in their signal properties. They have a more emergent onset, and less distinct P-and S-phase arrivals but a high peak frequency and more broadband spectra (Figure 9). Some families correspond to events which occur exclusively within a single day, whereas others represent more long term recurring events over several months ( Figure 4). Unlike the June families in which the families were all very cohesive and self similar with high cross correlation values, the November swarms show very little similarity to one another ( Figure 10B). The Q ARMA values show the events from September onwards have a typically higher Q ARMA value and peak frequency, than those in the June swarms (Figures 6, 11). Two significant located VTs occurred during November-a M3.3 event at 13:32 on 14th November, and a M3.6 event at 22:30 on 27th November. These are the two largest VT events at Cayambe during 2016. We manually pick events immediately before and after each of these M > 3 events to be sure we have an accurate picture of all seismicity which may be associated. We carefully model the rates of events following the M3.6 and M3.3 events; however, there is not a clear mainshock-aftershock pattern that can be modeled by a power (Modified Omori) or exponential relationship (Utsu et al., 1995;Bell et al., 2018). Family 1 groups a number of events around 14th November. However, it is unlikely that Family 1 exclusively represents this slip event as Family 1 spans October to December, indicative of a longer term process. For the slightly larger event on 27th   Figure  S2). We manually pick these high frequency spikes within the 5 min after 23:58:37UTC. However, station ANGU appears to have a defect and records recurring high frequency response between 20 and 30 Hz even during events known to be predominantly low frequency ( Figure 9C). As these events cannot be identified at stations CAYA or CAYR, we rule out any possibility that they are dynamically triggered events. As well as looking within the surface wave arrivals, we look for changes in average seismicity rate in the 24 h following the Pedernales event, for larger scale changes. We manually picked events as the catalogue only contained five event pick times, but visual observations of the file showed there were likely more. We picked at station CAYR, as this had the best signal to noise ratio. Overall there is a decreasing rate of seismicity through the 24 h period following the earthquake, from 30 events per hour, to less than 15.

The Meaning of Q
The Q values returned by all three methods span two orders of magnitude, which raises questions about the suitability of each approach. All methods were developed and tested for the purpose of low frequency seismicity, and specifically LP events generated by a resonating fluid filled crack (Kumagai and Chouet, 1999;Kumagai and Chouet, 2000). These methods are hinged on the spectral peak frequency of a monochromatic or harmonic signal (Lipovsky and Dunham, 2015). As such, the results of the Q analysis on broadband VT seismicity are not meaningful in a physical sense. Particularly in the case of Q ARMA -a metric based FIGURE 12 | Event "durations" defined by the τ parameter in the Q decay method. Events are separated into three broad populations identified by families: red marks persistent LP seismicity, blue marks the June VT swarm events, and black the VT swarm from September onwards.
Frontiers in Earth Science | www.frontiersin.org June 2021 | Volume 9 | Article 680865 on a single peak frequency in a broadband VT spectra, tells us very little. However, we can use Q decay as a quantifiable metric to describe the signal envelope. Determining a signal envelope that can be used to distinguish event types is often a subjective matter, as a certain degree of smoothing must be included. One more quantified approach is to define a relationship between the peak amplitude of a signal and a minimum "threshold" (Lees et al., 2008). The Q decay approach determines how quickly the signal amplitude decreases to 1 e of the maximum. This decay time, τ, shows a significant difference between the VT swarms in September and the persistent LP seismicity (Figure 12).
The Q ARMA method is also reliant on quantifying the growth rate (g) of a singular fundamental frequency (f). However, if there are insufficient points clustered at a pole at the peak spectral frequency, a Q ARMA value is not calculated. This reduces the number of erroneous Q values generated for broadband VTs. With this in mind, we use Q ARMA as the primary Q metric for interpreting this sequence of seismicity. For the analysis of LP seismicity this brings this study in line with previous efforts to monitor changes in Q across a series of unrest (Métaxian, 2003;Molina et al., 2004). And for VT events, if a Q ARMA value is calculated it is a function of the most dominant frequency with a broadband signal and used for comparison purposes, but not to interpret a physical process.

Volcano-Tectonic Swarms in 2016-2017
The temporal and spectral analyses do an excellent job to help quantify the different types of seismicity observed throughout one year at Cayambe. The two swarms in June and September onwards are related to two very distinct processes. The June swarms show clear separation between P-and S-phase arrivals, and have a distinct signal envelope in comparison to other families of VT events, where the signals are more emergent. Whilst there is significant uncertainty in the depth and magnitude estimates, the locations detailed in the IGEPN catalogue all show the June swarm events to be in the northern and eastern flanks of the volcano-quite separate from those later in the year. The high XCC of events in the June swarm, but dissimilarity to the September swarms suggests they have different source mechanisms. The events in June 2016 show a very different spectral signature to those later in the year, with much lower peak frequency and Q ARMA values. This may be the result of both source and path effects. Unconsolidated or hydrothermally altered flank deposits can act to extend the coda of failure events (Bean et al., 2014), which for proximal VT swarms may generate an increased Q value compared to a distal VTs. The disparity may also be dependent on depth and proximity relative to station CAYR. This June swarm was more likely associated with movement from the regional fault line along "La Sofía-Río Chingual." The Chingual fault is a major regional fault that extends 70 km from Cayambe, north through the Andes into Colombia (Ferrari and Tibaldi, 1992;Eguez et al., 2003;Alvarado et al., 2016) (Figure 1). As a named fault it is frequently cited in relation to earthquakes and volcanic systems in the region (Samaniego et al., 2005;Guillier and Chatelain, 2006;Yepes et al., 2016). Although it should be noted that the local tectonics in this region are complex and many smaller faults are not well constrained (Tibaldi et al., 1995). The last major tectonic earthquakes to occur in this region were in 1987 (M6.9 and M6.7), however several smaller thrust and right lateral strike slip events have occurred since (Kawakatsu and Cadena, 1991;Guillier and Chatelain, 2006;USGS, 2021). Unlike some of the swarm activity later in the year, these events are well recorded across all three permanent stations. Particularly at station CAYA, which lies closer to the Río Chingual and 5 km north of the Nevado Cayambe central crater.
The seismicity from September onwards, however, is more likely to originate from beneath the Nevado Cayambe central cone. The locations of the events are much closer to the central cone and to the south and west flanks. The events from the catalogue and subsequent template matching search show a more emergent signal envelope than those in June, with no clear separation between P-and S-phase arrivals. These events have higher peak frequencies overall and higher Q ARMA values. These VT events are quite LP-like in this respect. Whilst some families (15,17,22) represent events exclusive to just a few days, the majority span more than 3 months, indicative of a long term source process rather. Expansion of the template matching into 2017 shows a number of these families persist to the end of VT swarm activity in September 2017. However, unlike in June, the families are all quite dissimilar to one another suggesting there are multiple sources or the source is more heterogeneous or changing. The VT families after September 2016 are also less clear at station CAYA unless tens or hundreds of events are stacked, suggesting the origin is perhaps closer to the central cone and stations ANGU and CAYR. We believe these proximal VTs in late 2016 may be associated with progressive magma ascent and the formation of a proto-conduit (Fournier, 2007;White and McCausland, 2019). Whilst some families of events represent seismicity from a single day and perhaps correspond to a single rupture event (Supplementary Figure S1C), more long term persistent proximal VTs indicate an ongoing shallow process.

Origin of Long-Period Seismicity
The long period seismicity is informative as a baseline measure of day-to-day activity at Cayambe. Each event from both the catalogue and the template matching shows remarkable similarity throughout the year with a strong single peak frequency at 2.1 Hz and no strong additional harmonics. The LP events have typical durations between 10 and 15 s, and so cross correlation calculated in 10 s windows captures a thorough representation of the whole signal for these specific events. The LP seismicity is not well recorded on either ANGU or CAYA, and with no visible distinction between P-and Sphases, locating or magnitude estimating these events is not possible. Given these events are so well captured at only station CAYR, we suggest they are likely of a shallow source very close to CAYR within the Nevado Cayambe cone.
We consider all the possible source mechanisms for the persistent, repeating, LP seismicity in this episode of unrest at Cayambe. Low frequency, and in particular repeating low frequency seismicity have been described and modeled around the world. Minakami et al. (1951) first identified similar repeating events at Usu Volcano by manually measuring and matching peaks in the waveforms. These C-type earthquakes were attributed to a repeating source mechanism. Since these early efforts to categorise and interpret sequences of LP seismicity, a diverse range of LP seismicity has been described, and as such, there are many well established models for source mechanisms. Initial understanding of LP seismicity stemmed from studies of tremor and an idealised monochromatic or harmonic resonance process (Chouet, 1981). LP seismicity is inherently related to fluids and pressure changes and therefore their presence in the record is of interest for volcanic monitoring and forecasting. Generally, it is accepted that LPs are associated with the resonance of a fluid filled crack (Chouet, 1988). This has been extensively modeled by varying the crack media and investigate changes in the resultant waveform (Kumagai and Chouet, 2000). Further theoretical work has sought to model the initial excitation mechanism of the resonance process. For example, with the transportation of magma, flow-induced oscillations may excite a resonance mechanism (Julian, 1994). Whilst other studies have linked LP seismicity with the tensile hydraulic fracturing of rock (Bame and Fehler, 1986;Fournier, 2007;Lipovsky and Dunham, 2015). Due to network limitations at Cayambe, we cannot locate the LP events, but it is important to consider that shallow and deep LP sources may be connected in a volcanic plumbing system and triggered by pressure variations (Shapiro et al., 2017). It is important also to consider that the volcanic edifice is very heterogenous and LP seismicity may be the result of failure in unconsolidated materials rather than fluid presence (Bean et al., 2014). In this particular study we have a relatively small number of events that are captured clearly on just one station, and therefore can only make speculative suggestions around the source mechanism involved. However, if we consider the temporal patterns and rates of the repeating LP seismicity, then previous studies of LPs in volcanic systems have presented models of cyclical, periodic repeating events generated by magma emplacement and conduit processes (Petersen, 2007;Bell et al., 2018). Whereas episodes of LP seismicity from stick-slip glacial sources have shown swarm like behavior correlating with periods of snow fall (Allstadt and Malone, 2014). The LPs at Cayambe demonstrate neither of these extremes, and have an average periodicity ∼1, suggesting a constant average event rate and exponentially distributed interevent times. The LP event rate is also low compared to observations of seismicity at similar intermediate arc volcanoes. At Tungurahua volcano, hundreds of events per hour were recorded during activity in 2015 (Bell et al., 2017) and prior to explosions at Volcán de Colima, LP seismicity rate was in the range of hundreds of events per day . However, although neighboring and compositionally similar, these are two very active volcanoes. Rates of seismicity in the order of 10 0 -10 1 per day are more closely aligned with Puracé, Colombia (Arcila and Arcila, 1996) and Kawah Ijen, Indonesia  which like Cayambe, have not experienced recent explosive activity. Even in a scenario where LP seismicity is believed to be caused by glacial slip as opposed to volcanic processes, at Cotopaxi volcano, events were recorded at a rate of 60 events per hour (Métaxian, 2003). In the four instances where periodic LPs from the families are detected, the average rate at Cayambe is just one event per day. The LP event distributions in frequency and magnitude are also largely unaffected by the swarms in June and September, and show no change following the Pedernales M w 7.8 earthquake.
A study on Cotopaxi Volcano, Ecuador, used Q values to establish whether LP earthquakes were glacial or volcanic in origin (Métaxian, 2003). Using ARMA methods from the signal waveforms and estimating Q from the impedance contrast of ice and water, Q values 1-3 were identified as being glacial. However, in the records at Cayambe, the LP events generate Q ARMA values with a mean value 9.9. We use the same procedure (Aki et al., 1977) and consider a waterandesite contact with approximate andesite bulk density and V p values of 2,800 kg m −3 and 3,000 ms −1 respectively, where Q is approximately 9.5.
If we believe the LPs are associated with a hydrothermal system, then we may anticipate long term trends in the seismicity, and so we extend the template matching to identify LPs in 2015 and 2017 ( Figures 6G-L). The rate of seismicity is unaffected by continued VT swarm activity in 2017. Key identifying metrics including Q ARMA , peak frequency and high XCC values suggest these are from the same source as the 2016 LPs. Studies of Glacier 15α on Antisana volcano note that it is a good analogue for other glaciers in the Eastern Cordillera in Ecuador. This is due to their position on the equator, and direct exposure to the humid winds from the Amazon basin immediately to the East (Manciati et al., 2014). Unlike volcanic glaciers closer to the tropics, that are subject to extremes in climate such as El Niño and La Niña (Ceballos et al., 2006), the equatorial glaciers do not experience extremes of seasonality in precipitation, winds or temperature (Manciati et al., 2014;Basantes-Serrano et al., 2016). As such we might expect sub-glacial hydrothermal systems to be relatively consistent all year, rather than pronounced seasons of glacial unloading and meltwater recharge (Jónsdóttir et al., 2009). If the LP seismicity is related to hydrothermal interactions beneath the glacier, it is important to understand how these react-if at all-to external volcanic or tectonic seismicity.
Very low frequency, harmonic and tornillo-type seismicity was not easily identified and sorted by cross correlation. This is likely due to the duration of the events, and inconsistent peak frequency of these signals. The cross correlation method is applied over a fixed length, 10 s window for all events. This is suitable to capture similarities in VT and LP events that are up to 30 s in length for example, but does little to define events that last several minutes. In the harmonic seismicity, we see some examples of gliding frequencies and sequences of discrete events with different peak spectral frequencies ( Figures  7D-F). If a harmonic signal has only one very dominant spectral peak, this frequency can very significantly control the envelope and duration of the signal. This will limit the effectiveness of a simple, sliding window, cross correlation method, even with a very low detection threshold. As such these low frequency harmonic events should be carefully monitored, as they may be missed by automatic classification methods. Processes such as the harmonic product spectrum (HPS) (Roman, 2017) or as simple automated search for strongly peaked frequencies should be used to identify longer episodes of seismicity. Tornillo events and harmonic tremor at active volcanoes can be indicative of very shallow fluid processes and are often indicative of imminent activity (Hellweg, 2000;Hotovec et al., 2013;Arámbula-Mendoza et al., 2016). Highly monochromatic tornillo events at Galeras Volcano, Colombia, showed strong correlation between event duration and the time until eruption (Narváez et al., 1997). Whilst very long period (VLP) and LP events at Cotopaxi and Tungurahua, Ecuador were attributed to resonance of gas and ash filled cracks (Molina et al., 2004;Molina et al., 2008). These observations of tornillos are all quite different to the tornillo and harmonic seismicity documented here at Cayambe. The tornillo signals at Cotopaxi and Tungurahua were generally all very similar to one another, consistent or sequential in their dominant forcing frequency, and occurred during heightened activity or eruptions. The varied character of the harmonic seismicity at Cayambe, however, suggests this is not a stable or well established source process; without locations it is not possible to confirm whether these are from multiple sources. The first documented strong sulphurous smells near Picos Jarrín on the southwestern flanks of Cayambe, are in a Special Report on 1st December (IGEPN, 2016). However, without any quantified constraints on gas flux, we cannot be certain of the origin of this very heterogeneous harmonic seismicity. Due to the duration and pulsed like nature of some of the harmonic signals, they may also be classified as harmonic tremor (Arámbula-Mendoza et al., 2016;Bell et al., 2018). It is important to consider that anthropogenic sources such as nearby helicopters are known to generate harmonic tremor, and must be carefully identified and excluded from analysis -although this is often at much higher frequencies >10 Hz (Eibl et al., 2015).

Dynamic and Static Triggering
Large tectonic earthquakes cause stress changes in the crust, that may subsequently interact with volcanic systems. The stress change may be large and temporary (dynamic triggering) or small but permanent (static triggering) (Hill et al., 2002b;Manga and Brodsky, 2006). This relationship is complicated and causality versus coincidence is frequently debated (Lemarchand and Grasso, 2007). We use the seismicity rate as an indicator for dynamic triggering, given there is no eruption in this sequence. Following the Pedernales M w 7.8 earthquake, manually picked events at station CAYR do show a decreasing rate overall. It is significant that there is no dramatic increase in LP seismicity, as previous dynamic triggering studies have linked tectonic earthquake surface wave arrivals to instability in shallow magma storage (Lemarchand and Grasso, 2007;Lin, 2017). Although there are many more contributing factors, it has been shown that often shallow fluid storage does not necessarily make a volcanic system more susceptible to dynamically triggered activity (Prejean and Hill, 2018).
Given the first VT swarm begins at Cayambe just 7 weeks after the Pedernales earthquake, we speculate whether a small static stress change in the crust may have played a role in the resumption of activity. One study by Walter and Amelung (2007) presents a model for triggering of eruptions at arc volcanoes following megathrust earthquakes. Four case studies are presented from Sumatra, Kamchatka, Alaska and Chile, where a subduction earthquakes occurs within 200 km normal to a volcanic arc. Elastic numerical modeling of the volumetric strain demonstrates that each of these volcanic arcs would likely have been in a dilational area of volumetric expansion. This may in turn cause small static changes in the stress field and perhaps promote new eruptions. For example, a magma reservoir in an area of dilation can experience overpressure driven by magma mixing and exsolution of volatiles (Sparks et al., 1977). We believe the geometry of this model is paralleled at Pedernales and Cayambe. Further to this Béjar-Pizarro et al. (2018) model the stress change from the CMT focal mechanism of the Pedernales mainshock, and suggest that Cayambe lies in a region of "high level potential activation" where induced normal stress (Δσ) <0.1 bar.

Re-Awakening Dormant Volcanoes
The sudden increase in seismicity at Cayambe has implications for long term forecasting, particularly given there is no record for recent eruptions. One possible model for the re-awakening of dormant stratovolcanoes presents a 4-stage process . Different phases of characteristic seismicity can be attributed to each stage-from deep seismicity, to distal VT events, to vent clearing and then finally shallow repetitive seismicity ahead of an eruption. The study presents this model exemplified on volcanoes across timescales from days to years. Stage 1 including deep seismicity is associated with lower crustal magma storage. As is the case on many volcanoes, the limited three station network on Cayambe, and in turn the uncertainty in depth estimates mean it is not possible to constrain with confidence any VT seismicity between 10 and 40 km depth. However, it is possible that the variety of seismicity recorded at Cayambe in 2016 and 2017 may be indicative of some stage 2 and 3 re-awakening, described in this model.
The VT swarm during June 2016 may be considered distal seismicity and is consistent with the stressing of a regional fault system, described in stage 2 of White and McCausland (2019). We speculate that the M w 7.8 Pedernales earthquake generated a small static change in the stress field at Cayambe (Figure 13). This may have promoted a new pulse of magma to ascend (Walter and Amelung, 2007), or perhaps pressurized local aquifers and allowed slip along a pre-stressed fault (White and McCausland, 2016;Coulon et al., 2017).
The subsequent VT swarm from September 2016 onwards recorded events located within the Nevado Cayambe cone, and may be considered proximal VTs. If the distal VTs were associated with a new pulse of magma, then the proximal VTs may document the progression of ascent and possibly the formation of a proto-conduit. Repeated long-term families of VTs may indicate repeated rupture and sealing events of the brittle-ductile transition (Fournier, 2007), whilst isolated families may be indicative of individual rock fracturing episodes to accommodate ascent of magma and gases (Hill et al., 2002b). VT events that could not be classed into families may also be the spasmodic and irregular VT seismicity described in White and McCausland (2019) and associated with rupturing the brittle-ductile seal. We do not see accompanying low frequency seismicity at this stage, although this is possibly due to limitations in the network. The model for new activity at Cayambe is not wholly reliant on a new pulse of magma. It is also possible that a dilational regime in the crust prompted further exsolution of volatiles in pre-existing magma storage (Sparks et al., 1977). This is an option explored in a similar model in Roman and Cashman (2018).
The next stages of the onset are eruption are vent clearing and shallow repeating seismicity, as magma reaches the shallowest point in the conduit and eventually the surface. Evidence from the seismic record, however, does not suggest Cayambe volcano reached this stage in 2016-2017. The LP seismicity in particular is not typical or indicative of these late stages of reawakening. Empirical calculations of acoustic impedance between magma-rock contacts would suggest Q values for these types of LPs should be an order of magnitude higher (Kumagai and Chouet, 2000;Métaxian, 2003). We also might expect these events to have larger amplitudes and more prominent cyclical trends as magma ascends and interacts with the conduit margins (Iverson, 2008;Bell et al., 2018;Butcher et al., 2020). Whilst we may speculate some harmonic tremor and isolated tornillos are coincident with records of sulphurous smells up on the flanks of Cayambe (IGEPN, 2016), they are not characteristic of rapid magma ascent or gas venting processes in the latter part of stage 3 and stage 4 of the model (Hotovec et al., 2013;Arámbula-Mendoza et al., 2016). During the crisis, IGEPN reported deformation data but neither GPS or tiltmeter readings suggested any change beyond the margin of any error, ±7 mm (IGEPN, 2016). This further supports the idea that whilst there may have been magma ascent, it was not particularly rapid or shallow enough to cause significant flank deformation. The VT swarms show a decreasing rate until the end of 2017 (Figure 2). There has since been no significant increase in the rate of LP seismicity, indicative of further magma ascent and a transition to stage 3 . Whilst it is possible that exsolved gases in an ascending magma body will act to increase volumetric expansion, exsolution of water will act to stall ascent as the viscosity of magma increases (Dixon et al., 1995). It is possible that exsolution of water played a role, or simply the static stress field change was not enough for a magma pulse to ascend to the surface at Cayambe.
We suggest the persistent LP baseline seismicity recorded during 2016 and 2017 is of hydrothermal origin ( Figure 13). As such any increase in LP seismicity with a significantly different envelope, Q ARMA or spectral signature; coupled with an increased in volcanic VT type events should be carefully monitored.
Harmonic tremor and tornillo events should also be carefully monitored, coupled with DOAS measurements to identify significant changes in degassing. Periodically running cross correlation analyses and updating master event signals is important to trace subtle changes in seismicity and incorporate more recent events into the monitoring and interpretations. The 4-stage process  presents observations of the model in action across time scales from 3 to 4 3months at Mt St Helens (Thelen et al., 2011) to nearly 4 years at Mt Sinabung . It is therefore not unusual to observe this slow transition across a period of several years. The models for reawakening are not necessarily linear and if the process of reawakening did stall in 2017, then any resumed activity may include stage 2 distal VTs as FIGURE 13 | Schematic diagram of seismicity at Cayambe in 2016-2017, adapted from the crustal stress model presented in Walter and Amelung (2007). Map inset shows CMT focal mechanism for the M w 7.8 Pedernales earthquake (Dziewonski et al., 1981;Ekström et al., 2012) and relative location to Cayambe. Vertical scale on right, only an indicator of magmatic processes and not applicable for subduction depths.

Implications for the Future
Low latitude glaciers are rapidly shrinking (Thompson et al., 2011). Evidence from analogue glaciers on Antisana suggests that whilst glacial retreat has been ongoing for over 150 years, this has accelerated in the last 30 years, as a result of anthropogenic climate change (Manciati et al., 2014). Studies in Iceland have shown that unloading effect from ice retreat can impact stress conditions within a volcano (Sigmundsson et al., 2010), and cause long term changes to failure conditions in shallow magma storage. Rapid unloading may then promote further magma ascent and impact pre-stressed faults in the volcanic ediface. Station ANGU, although with it's defective instrument response, is installed on the very edge of the western reaches of the summit glacier. Carefully monitoring for any changes here may be particularly important. Even very localised slip may have further implications for secondary hazards from the glacier itself. Furthermore, in terms of hazard assessment, there are immediate hazards present from ice retreat and rapid glacier melt, particularly during an eruption (Tuffen, 2010). It is estimated around 4.2 km 2 of ice was lost during the 1985 eruption of Nevado del Ruiz, contributing to devastating secondary lahars (Ceballos et al., 2006). In the short term, increased meltwater also has the potential to interact further with hydrothermal systems, increasing pore pressure. If we consider Cayambe Volcano to be re-awakening, where the next stage includes further magma ascent and shallow storage then additional meltwaters could lead to potentially very dangerous phreatic activity. Using the existing network, close monitoring of the LP seismicity could be crucial. Due to the unique positioning of glaciers in the Ecuadorian Andes, Cayambe is not subject to seasonal extremes in climate and so the record of seismicity from the hydrothermal system is consistent all year round. Slight changes in the Q ARMA value, or spectral signature could then be indicative of a change to what appears to be a very stable hydrothermal system in balance through 2016 and 2017. All of these recommendations for monitoring signal properties, can only be reliably determined from one station. An extended network of seismometers on Cayambe would allow for greater accuracy in locating events, to underpin interpretations from the event characteristics.

CONCLUSION
This study identifies several unique subsets of volcanic seismicity, observed at Cayambe during the seismic crisis of 2016. We explore the potential for re-awakening at a volcano that has been seismically restless for 20 years but last erupted in 1785. The seismic crisis at Cayambe volcano occurred just 3 months after the M w 7.8 earthquake, 200 km west on the Ecuador coastline. We draw comparisons from models of dynamic and static triggering of volcanic eruptions to the seismicity documented in 2016 and 2017. Should the seismicity have been generated by a new pulse of magma beneath Cayambe volcano, then the process of reawakening is likely ongoing, although slowed or stalled currently. Given the potentially vulnerable state of the glacier, the unknowns around eruption dynamics and significant populations living within the immediate vicinity, Cayambe remains a dangerous volcano. We believe that using the metrics presented in this study, however, it is possible to carefully monitor the seismicity on the existing permanent network, and quickly identify when earthquake signals are evolving.

DATA AVAILABILITY STATEMENT
Seismic data are curated by the IGEPN and available on request (from https://www.igepn.edu.ec/solicitud-de-datos).

AUTHOR CONTRIBUTIONS
All authors conceptualised the research project. MR and SH provided geophysical monitoring data from Cayambe volcano (curated by IGEPN). AB supervised the research. SB lead the data analysis. All authors contributed to writing the article.

SB was funded by Edinburgh Earth and Environment Doctoral
Training Partnership (E3 DTP) Grant NE/L002558/1. AB was supported by NERC GCRF award NE/S009000/1.