Renewed Explosive Phreatomagmatic Activity at Poás Volcano, Costa Rica in April 2017

Salvage, Rebecca O.; Avard, Geoffroy; de Moor, J. Maarten; Pacheco, Javier F.; Brenes Marin, Jorge; Cascante, Monserrat; Muller, Cyril; Martinez Cruz, María

doi:10.3389/feart.2018.00160

ORIGINAL RESEARCH article

Front. Earth Sci., 16 October 2018

Sec. Volcanology

Volume 6 - 2018 | https://doi.org/10.3389/feart.2018.00160

This article is part of the Research Topic Towards Improved Forecasting of Volcanic Eruptions View all 21 articles

Renewed Explosive Phreatomagmatic Activity at Poás Volcano, Costa Rica in April 2017

$\r\nRebecca O. Salvage,*$ Rebecca O. Salvage^1,2^*

Geoffroy Avard¹

J. Maarten de Moor¹

Javier F. Pacheco¹

Jorge Brenes Marin¹

Monserrat Cascante^1,3

Cyril Muller¹

María Martinez Cruz¹

¹Observatorio Vulcanológico y Sismológico de Costa Rica, Universidad Nacional, Heredia, Costa Rica
²Department of Geoscience, University of Calgary, Calgary, AB, Canada
³Department of Earth Sciences, University of Oregon, Eugene, OR, United States

Phreatic and phreatomagmatic eruptions at volcanoes often present no short term precursory activity, making them a challenge to forecast. Poás volcano, Costa Rica, exhibits cyclic activity with phreatic and some phreatomagmatic eruptions separated by times of quiescence. The latest phreatomagmatic stage began in March 2017 with increases in crater lake temperatures, SO₂ flux, and the rate of seismicity, as well as accelerated ground inflation near the active crater. On 23 April 2017 at 04:12 UTC, a large phreatomagmatic eruption occurred at Poás, sending blocks up to 1 m in length to distances >1 km. Hindsight analysis revealed a precursory seismic sequence from 25 March to 22 April of similar seismic events (in terms of their frequency and waveform characteristics). Fourteen families of similar seismic events (containing ≥10 events per family) were identified during this precursory sequence, totaling over 1,300 events. An acceleration within the dominant family of LF (low frequency) waveforms was identified, suggesting that a forecast for the onset of the eruption may have been possible using the Failure Forecast Method (FFM). However, no confidence could be placed in the forecast generated, reiterating that not all accelerating trends are suitable for analysis using the FFM, in particular in conjunction with a least-squares linear regression. Our residual analysis further supports the concept that using a least-squares linear regression analysis is not appropriate with this dataset, and allows us to eliminate commonly used forecasting parameters for this scenario. However, the identification of different families of similar seismicity allows us to determine that magmatic fluid on its way to the surface initially became stalled beneath a chilled margin or hydrothermal seal, before catastrophically failing in a large phreatomagmatic eruption. Additionally, we note that 24 h prior to the large phreatomagmatic eruption, all LF families became inactive, which could have been falsely interpreted in real time as the waning of activity. Our results suggest that identifying families of seismicity offers unique opportunities to better understand ongoing processes at depth, and to challenge conventional forecasting techniques.

1. Introduction

Phreatomagmatic eruptions at volcanoes occur when magmatic fluid or gases migrating toward the surface interact explosively with ground or surface water. At the point of direct contact between the two, the magmatic fluid or gases, which are of a higher temperature than the static water, rapidly cools, and the static water abruptly heats and expands due to a rapid phase change from liquid to gas, generating an explosive eruption (Büttner et al., 2002). The products of these types of eruptions are usually fine ash deposits (as a result of the energetic fragmentation process of rising magma and the subsequent fragmentation of the surrounding edifice; Büttner et al., 1999), hydrothermal fluids, and large blocks in the vicinity of the erupting vent. Lahars are also common, depending on the amount of water at the surface and the local topography (Barberi et al., 1992), meaning these eruptions can pose a serious hazard to local populations.

Precursors to phreatic eruptions are often short-lived, and are often difficult to distinguish from normal background levels. The September 2007 eruption of Mt. Ruapehu, New Zealand, produced a steam column up to 15,000 feet, ballistic and surge flows, and lahars without any usable precursors indicative of an impending eruption (Christenson et al., 2010; Jolly et al., 2010). Similarly, the September 2014 eruption of Mt. Ontake, Japan, which left 64 people dead or missing, showed very little precursory activity with only small changes in the amplitude of recorded tremor, a small migration in high frequency seismicity and anomalous tilt measurements <10 min before the onset of eruption (Kato et al., 2015; Oikawa et al., 2016), and one very long period earthquake in the preceding 25 s (Maeda et al., 2015). Poás volcano, Costa Rica, on which this study is focused, showed very short-lived precursory fluctuations in lake gas composition (between 24 and 36 h) prior to a number of phreatic eruptions between April and June 2014 (de Moor et al., 2016). In other cases, precursors are identifiable but may not be able to be processed and analyzed in an appropriate time frame to give an alert, for example, the 1990 eruption of Kelut volcano, Java, which was preceded by spasmodic tremor in the hours prior to a number of phreatic explosions (Lesage, 1995). These types of eruptions are therefore particularly difficult to forecast in terms of timing and intensity, especially when lacking precursory activity, principally as groundwater reservoir locations and heat flow within the volcano are often poorly constrained. In addition, phreatic eruptions may not always involve the movement of magma toward the surface (which usually generates precursory signals) and may occur simply due to changes in the shallow hydrothermal system or variations in gas input (Rouwet et al., 2014; de Moor et al., 2016).

Here we report on the largest phreatomagmatic explosive eruption to occur at Poás volcano during the latest phase of unrest, which was first identified in real-time by the staff of the Observatorio Vulcanológico y Sismológico de Costa Rica Universidad Nacional (OVSICORI-UNA) on 26 March 2017, and lasted several months. A number of geophysical and geochemical parameters rapidly increased from this date onwards, suggesting a sudden renewal of activity, and a number of smaller phreatomagmatic eruptions were indeed identified in early April 2017, which became more energetic with time. As unrest was already elevated, further activity suggesting a much larger phreatomagmatic eruption was imminent was not identified in real time, but did occur on 23 April at 04:12 UTC. We show that although the similar seismicity identified at Poás volcano produced an accelerating trend, it was not possible to use it in conjunction with the Failure Forecast Method to generate a successful forecast. However, the techniques employed allowed the identification of far greater numbers of events compared to a traditional amplitude based detection algorithm or manual identification by OVSICORI-UNA personnel. The identification of families of similar seismicity also allowed us to gain a more detailed insight of ongoing processes at depth, including a better understanding of the number of active sources beneath the volcano.

2. Poás Volcano, Costa Rica

Poás volcano, Costa Rica (N 10.1968; W 84.2305, ~2,300 m a.s.l., Figure 1) is a large stratovolcanic complex located in the central valley of the country, 22 km north of the main international airport (SJO) and 35 km north-east of the capital city San José, where approximately one third of the countries' population live. Poás has three main volcanic centers: Botos cone, containing a cold, mildly acidic crater lake, which last erupted approximately 7,500 years ago (Prosser and Carr, 1987; Alvarado-Induni, 2005); the von Frantzius crater, an extinct volcanic cone (Casertano et al., 1983); and the principal crater, which is usually partially filled with the Laguna Caliente, a warm, very acidic crater lake (20–60°C, pH ≤ 1.8), which experiences frequent phreatic and more occasional phreatomagmatic eruptions (de Moor et al., 2016). The Laguna Caliente is the result of meteoric water interacting with hydrothermal-magmatic gases and rocks. Prior to its closure in April 2017, Poás was the second most visited national park in the country, with more than 400,000 visitors per year. Eruptions at this volcano consequently pose a great threat to those in this vicinity, and in particular to tourism, both in terms of health and economy.

FIGURE 1

Figure 1. (A) Topographical map showing the location of monitoring equipment before and after the 23 April 2017 large phreatomagmatic eruption at Poás (inset: Reference location of Costa Rica and Poás). The active crater (Laguna Caliente) is represented by the black diamond; and the dormant Botos crater is shown as a gray diamond to the south-east of the active crater. The Mirador outlook was the viewing station of the National Park for observation by tourists into the active crater, which was closed on 13 April 2017. Seismic and GPS stations are named as referenced in the text. CRPO, the Multi-GAS station (purple star) and the camera (green star) were all destroyed by eruptions in April 2017. (B) Photograph of Poás Laguna Caliente, 12 March 2017 from the Mirador before levels of unrest increased. Photo Credit: R. O. Salvage. (C) Example phreatomagmatic eruption at Poás crater, 19 April 2017 from the Mirador. Photo Credit: J. F. Pacheco.

Historical activity at Poás is reported as early as 1828 (strong degassing with sporadic phreatic eruptions; Casertano et al., 1983), with the first phreatomagmatic eruption reported on 25 January 1910, when a large steam and ash cloud reached 4 to 8 km above the summit (Martínez et al., 2000). The last period of strong phreatomagmatic activity occurred between 1953 and 1955 and saw the emission of both ash and larger bombs from the principal crater, followed by the complete drying up of the Laguna Caliente (Casertano et al., 1983). This eruptive period resulted in the creation of a small dome-like feature at the edge of the Laguna Caliente, growing up to 45 m high at its greatest, partly made up of primary pyroclastic deposits (Rowe et al., 1992). The crater lake returned by 1961 with sporadic geyser-like phreatic eruptions over the coming decades, sometimes with associated ash and rock emissions (Martínez et al., 2000). More recent activity occurred in 2005–2008 and 2011, 2014, and 2016 as reported by the OVSICORI-UNA (e.g., Global Volcanism Program, 2016). Intense fumarolic activity, in particular in and around the edges of Laguna Caliente, was reported between 1980 and 1985, 1986–1991, 1994, 1999, and 2005–2008. New fumaroles opened on the south-western inner crater walls between 1995 and 2000, and toward the east between 1999 and 2007, away from the Laguna Caliente (Martínez Cruz, 2008). Between 2009 and 2011 numerous phreatic eruptions occurred with column heights varying from 5 to 500 m. In the latter half of 2011, fumaroles around the dome reached temperatures up to 980°C. Further smaller phreatic eruptions occurred in April, September, and October 2012. On 1 and 2 May 2013, a series of phreatic eruptions occurred, which were preceded by a clear increase in seismicity from 29 April onward (≥200 low frequency events, above a background level of 50 events/day), as well as an increase in the temperature of active fumaroles on the dome from 102°C on 30 January to 380°C on 16 April (Fischer et al., 2015). In 2014, 46 phreatic eruptions with varying plume heights up to several hundred meters were detected using the seismic network (Avard et al., 2015). In early 2015, incandescence was observed around the dome, with fumaroles registering temperatures of 625°C, and mobile Multi-GAS measurements yielding high SO₂/CO₂ gas ratios from the acid lake. Renewed phreatic activity was observed at Poás between June and August 2016.

2.1. Monitoring Network

Systematic monitoring began at Poás volcano in 1978, conducted by the Universidad Nacional, presently known as the OVSICORI-UNA (Martínez Cruz, 2008). Prior to the volcanic crisis in April 2017, the volcano was monitored by 2 seismic stations in close proximity to the volcano, CRPO [Nanometrics Taurus digitizer with Trillium Compact sensor (sensitivity of 0.008–100 Hz)] at 550 m east of the dome, and VPVF [Guralp CMG-DM24-EAM digitizer with CMG-3ESPC sensor (sensitivity of 0.003–50 Hz)] 900 m north of the dome, and one seismic station (VPTE, Guralp CMG-DM24-EAM digitizer with CMG-3ESPC sensor) at 4.5 km from the volcano; 2 GPS (VPCR at 900 m north, and VPEV 1,350 m south-west of the dome); 1 stationary Multi-GAS station (300 m north of the dome); and a web cam (470 m north of the dome) (Figure 1A). The dome (which was subsequently destroyed during this eruptive episode) was located on the very southern edge of Laguna Caliente. Two additional seismic stations, VPLC and VPEM (both Quanterra Q330HRS digitizers with Nanometrics Trillium Compact sensors), were installed following the eruptions in April 2017. All seismic stations are broadband, 3-component sensors with 24/26 bit digitizers. Regular visits to the volcano by staff of OVSICORI-UNA allowed mobile DOAS measurements on foot along the western crater rim (remote sensing of SO₂ in the atmosphere based on UV absorption spectra); in-situ FLIR (thermal imaging); and sampling of lake water, fumaroles and eruptive products from within the principal crater containing Laguna Caliente.

2.2. Renewed Activity: April 2017

From January to March 2017, the number of seismic events identified at Poás volcano began to increase from an average of 30 events a day in January, to over 100 events a day from 12 February onward, and over 200 events a day from 26 March 2017 (Figure 3A). Events were identified manually from the incoming seismic record by OVSICORI-UNA personnel. The number of events per day was not considered unusual until the end of March 2017. No surface manifestations were observed until 1 April 2017 when a new water-rich degassing fumarole appeared 100 m west of the active fumarolic field on the edge of Laguna Caliente, that was active for ~24 h. On 7 April, a boiling (~90°C) and highly acidic spring (with a diameter of 4 m and the same pH and chemical composition as the lake, known locally as a “Borbollón") opened 200 m south of the dome. Due to high concentrations of detected SO₂, the administrators decided to evacuate the Poás National Park on 9 April and kept it closed for the rest of the day. On 13 April, at approximately 02:00 UTC, a new boiling vent opened on the dome in the south-eastern corner of Laguna Caliente, generating a small lahar which affected the village of Bajos del Toro, located 7.5 km west of the active crater; ~10% of the dome was destroyed; and ballistics up to 25 cm in length were reported on the western crater rim at distances of 400 m. The park closed on this day, and has remained closed up to the time of writing (September 2018). A phreatomagmatic eruption on 13 April at 21:45 UTC generated a plume of ash, water and magmatic gas to 500 m above the crater, destroying the Multi-GAS station and web cam in the process (stars, Figure 1A). On 14 April, another phreatomagmatic eruption occurred at 13:57 UTC, generating a plume that rose up to 4 km high and was visible from the capital city, San José (35 km away) as well as the closer provinces of Alajuela and Heredia, and destroyed a further ~80% of the dome. Ballistics up to 50 cm in length impacted the Mirador (national park tourist viewing platform) located 650 m from the main vent, at an elevation of 260 m above the crater. From 13-19 April, borbollón activity from the dome area, with sporadic energetic phreatomagmatic eruptions reaching elevations up to 500 m (Figure 1C) were observed, with erupted material progressively filling the lake around the active vent. On 21 April, this eruptive material reached the lake surface. The most energetic phreatomagmatic eruption of this eruptive episode occurred on 23 April at 04:12 UTC, launching ballistics up to 1.5 km away from the active vent. Spatter bombs up to 3 m in length were found 450 m away on the south-east crater rim. Approximately 80% of the forest within 60 m of the crater was completely destroyed, as well as a seismic station (CRPO, Figure 1A). Small, discontinuous eruptions continued until the end of May 2017, when the lake completely disappeared. On 6 June, Poás began a phase of passive ash emissions (i.e., non-explosive) and continuous magmatic degassing. Since the closest seismic station had been destroyed in the 23 April eruption, 2 new seismic stations (VPEM and VPLC, Figure 1A) were installed on the flanks of Poás in June 2017 to strengthen the monitoring network, in addition to two new web cameras and two new Multi-GAS monitoring stations.

2.2.1. Seismicity

The number of low frequency earthquakes (LFs, ≤ 5 Hz, Figures 2A,B), tremor and volcano-tectonic earthquakes (VTs, high frequency, 5–15Hz, Figures 2C,D) started to increase in January 2017. However, the values at the time were not seen as unusual based on Poás' previous eruptive episodes (Figures 3A,B). In hindsight, these events probably represent the renewal of activity at Poás. On 26 March, a total of 22 high frequency VT events registered over 11 h, although none of the events were locatable within reasonable errors. This was the first warning recognized in real time of renewed activity at Poás. The Real-time Seismic Amplitude Measurement (RSAM, Endo and Murray (1991)), was calculated by taking the average amplitude of seismic events every 10 min within two discrete frequency bands (1–5 Hz and 5.1–15 Hz) to analyse amplitudes related to both LF and VT events, respectively. Some authors use the term SSAM to describe the use of discrete frequency bands (e.g., Cornelius and Voight, 1994; Rogers and Stephens, 1995). We chose the more general RSAM whilst explicitly stating the frequency bands of computation. No significant change in the filtered high frequency RSAM 5.1–15 Hz (Figure 3F) was observed, suggesting that although there was an increase in VT event count at the end of March 2017, there was not a change in the energy released. A significant increase in VT seismicity (a three-fold increase in event count) was registered from 13 April onward (Figure 3B), which coincided with the onset of daily phreatic activity at Poás and an increase in high frequency seismic energy release (Figures 3B,F, light blue shaded areas).

FIGURE 2

Figure 2. Time series and frequency spectrum (normalized) of representative waveforms identified at Poás in March 2017. (A,B) Low frequency event, 29 March 2017 00:35:45, band-pass filtered between 1 and 15 Hz. (C,D) High frequency event, 26 March 2017 09:23:10, band-pass filtered between 1 and 15 Hz. (E,F) Very low frequency event, 27 March 2017 13:47:30, band-pass filtered between 0.5 and 3 Hz.

FIGURE 3

Figure 3. Time series of monitoring data at Poás from September 2016 to 30 April 2017. The green dashed line indicates the timing of the VT swarm on 26 March 2017; the light blue shaded areas indicate the timing of phreatic and phreatomagmatic eruptions from 12 April onward; the red dashed line indicates the timing of the large phreatomagmatic eruption on 23 April at 04:12 UTC. (A) Daily counts of LF (low frequency) seismicity. The daily counts in February and early March were not considered unusual. (B) Daily counts of VT (high frequency) seismicity. Very few seismic events were identified before March 2017. The VT swarm on 26 March is clearly visible (green dashed line). (C) GPS height change in meters (vertical component, station VPEV). A clear inflation of 4–5 cm can be observed from early March to 23 April. The large deflation signal in January 2017 is considered to be an outlier, as measurements returned to what was considered background levels the next day. (D) SO₂ flux, measured by walking a transect along the western crater rim. (E) 10 min filtered RSAM (1–5 Hz) recorded at station VPVF, 1 March to 30 April 2017. (F) 10 min filtered (5.1–15 Hz) RSAM recorded at station VPVF, 1 March to 30 April 2017.

On 28 March, a significant number of LFs and short duration tremor was registered, with more than 400 events a day following a systematic increase in LF events over the previous few days (Figure 3A). Significantly, on the same day, the first Very Low Frequency (VLF) earthquake was recorded for this eruptive phase (Figures 2E,F). LF events once again dropped below 300 events a day on 10 April, although this was still considered to be elevated. From the 12 April onward, although the number of LFs did not appear to increase until 21 April (Figure 3A), the energy release significantly increased (Figure 3E, shaded blue area). The maximum number of LFs between February and April 2017 was recorded on 22 April, with ≥750 events on a single day (Figure 3B, red dashed line). Since the end of April 2017, seismicity at Poás has gradually reduced to almost no registered activity by April 2018.

2.2.2. Deformation

From January to February 2017, GPS sites did not register any significant deformation. A clear inflation of the edifice began in mid March, which inflated by ~2 cm by the end of March, and by a further 3–4 cm by 23 April 2017. Following the large phreatomagmatic eruption on 23 April at 04:12 UTC, a significant deflation of 2 cm was registered over a 48 h period, probably related to the expulsion of mass during the eruption, after which the crater began to slowly inflate once more (Figure 3C). In addition, from February to April 2017, the distance between VPEV and VPCR (approximately 2,160 m) increased by 18 mm, after 2 years of relative stability. Similar to observations at other volcanoes (e.g., Soufrière Hills volcano, Montserrat; Voight et al., 1998; Kilauea, Hawaii; Tilling, 2008; Volcán Santiaguito, Guatemala Johnson et al., 2014), a clear correlation is observed between the deformation signals and earthquake activity (Figures 3A,C), in particular during February and early March 2017, where a sudden increase in both the count of low frequency events (Figure 3A) and high frequency events (Figure 3B) coincides with an increase in GPS signal (Figure 3C), associated with the inflation of the edifice.

2.2.3. Gas Measurements

The permanent Multi-GAS station located near Laguna Caliente registered a change in degassing behavior in late 2016. Prior to the end of November 2016 the SO₂/CO₂ ratio had been decreasing over a period of more than a year. In early December 2016 the SO₂/CO₂ gradually started increasing, which is considered a signal of increasing likelihood of eruption at Poás (de Moor et al., 2016). Dramatic changes were noted in gas composition immediately before the eruptive phase in April 2017. On 29 March, the H₂S/SO₂ ratio measured with a stationary Mutli-GAS station dropped from an average value of ~2.4 in March 2017 to ≤ 0.01 on 31 March. On the night of the 31 March, the SO₂/CO₂ ratio increased from 0.04 to 0.16, before the first visible manifestation at the surface (new water rich degassing fumarole on the edge of Laguna Caliente). The increase of the SO₂/CO₂ ratio continued to 0.44 on 1 April and from then on increased linearly until 12 April, where the average value was 7.4, until destruction of the station on 13 April. The SO₂ flux appeared very low at ~20 t/d on 28 March (at the lower limit of the detection threshold); at 180 t/d on 4 April; at ~440 t/d on 10 April; at 1,500 t/d on 13 April; and 2,200–3,000 t/d after 24 April (measured with a DOAS instrument, on a drone as it was considered too dangerous for personnel to enter the crater).

2.2.4. Petrology

An energetic geyser-like eruption at Poás in 2016 produced ash with ~6% juvenile material content (modal analysis on stereoscopic microscope, phi = 1–2). SEM analysis of this ash sample showed irregular angular to sub angular shards, as well as a few rounded shards. The 12 April 2017 eruption that destroyed 10% of the dome presented 9% juvenile material in the ash deposit. The juvenile component increased to ~30% on 14 April; to 63% on 20 April; and to 85% on 22 April, clearly showing the increasing influence of a magmatic component in the evolving eruption. SEM analysis of material collected from the 23 April eruption exhibited a sharp, glassy surface texture. After the 23 April 2017 eruption, the juvenile content of erupted ash remained constant at ~80%, until it increased further on 16 June to 96% (Cascante, 2017). Bombs from the April eruptions were highly vesicular and porphyritic in texture, with a matrix abundant in phenocrysts such as plagioclase (~30%), pyroxenes and a small number of altered olivines, suggesting an intermediate composition of basaltic-andesite to andesite (Martínez et al., 2017).

2.2.5. Lake and Boiling Springs Geochemistry

The temperature of Laguna Caliente showed a two-fold increase between late March and mid April 2017, increasing from 35°C on 28 March to 41°C on 4 April; 54°C on 10 April; and 64°C on 13 April. After this, it was not possible to visit the lake for direct sampling and measurements due to safety concerns. At the same time, the lake level increased by 1 m between 28 March and 11 April, and by another 0.5 m by 13 April, despite a lack of rain during this time, suggesting the increasing addition of material and/or fluid to the base of the lake from active fumeroles and vents. A boiling spring appeared on the crater floor 180 m to the south of the acid lake issuing acidic and salty hot waters (93°C) with chemical characteristics (pH, salinity, anion concentrations) similar to that of the Laguna Caliente brines, but over saturated in cristobalite, a high temperature silica polymorph. The pH of the lake remained stable passing from 1.44 on 28 March to 1.39 on 13 April. The SO₄ ²⁻/F⁻ ratios in the lake waters decreased during the first 2 weeks of April 2017, indicating input of halogen-rich fluids into the lake, supported by the increasing lake levels over this time.

Dissolved unreacted SO₂ increased in the acid lake from 5 ppm on 28 March to 400 ppm on 13 April. After the April 2017 eruptions, the ultra acidic crater lake of Poás shrunk rapidly until it disappeared in June 2017, allowing the subaqueous fumaroles to discharge directly into the atmosphere. Several ponds of molten sulfur and sulfur cones were observed at the dried bottom of the crater after June 2017. Some ponds contained bright yellow molten sulfur, but one contained molten pyrite-rich black sulfur, indicating boiling temperatures between 113°C and 116°C, depending on the type of impurities present. A new lake started to form by mid January 2018, due to the gradual slowing of magmatic activity throughout the second half of 2017, and the high levels of precipitation between August 2017 and January 2018 as a result of several tropical storms and low pressure atmospheric systems affecting the country. The new acid lake had a temperature of 60°C and a pH of 0.60 in January 2018. This lake dried out during March 2018, and was absent for approximately 1 month.

3. Similar Seismicity at Póas Volcano

Seismicity has remained a primary monitoring tool for the detection of volcanic unrest because it can be remotely analyzed in real-time, often by an automated system and it is considered a highly valuable tool for decision-makers. The characterization of seismicity in a volcanic environment is traditionally based upon the signals' time and frequency characteristics: different bands of frequency relate to different active source processes at depth, which can be distinguished from one another, although the frequency bands associated with each process may overlap (Lahr et al., 1994). Low frequency seismicity (LFs, Figures 2A,B) are characterized by frequencies between 0.5 and 5 Hz; show emergent P- and S- wave arrivals (Chouet and Matoza, 2013); and have been linked to resonance of seismic energy trapped at a solid-fluid interface within a crack (e.g., Chouet, 1988) or a volcanic conduit (e.g., Neuberg et al., 2000). The trigger mechanism of such seismic energy may be generated by a stick-slip motion along conduit walls as magma ascends (e.g., Iverson et al., 2006; Kendrick et al., 2014); the brittle failure of the rising magma itself due to an increase in strain and viscosity rates (e.g. Lavallée et al., 2008; Thomas and Neuberg, 2012); interactions between the magmatic and hydrothermal system (e.g. Nakano and Kumagai, 2005); or through the slow rupture of unconsolidated material on volcanic slopes (Bean et al., 2014). High frequency events (VTs, Figures 2C,D), which are characterized by frequencies >5 Hz, have clear, impulsive P- and S- arrivals; and are generated when magmatic processes create enough elastic strain to force the surrounding edifice into brittle failure (Arciniega-Ceballos et al., 2003). Many volcanic seismic events fall between these two end-member categories and are classified as “hybrid” events. Very low frequency events (VLFs) typically occupy the frequency range below 1 Hz (Figures 2E,F). VLFs are typically attributed to the coalescence or ascent of gas slugs within a volcanic conduit during migration (e.g., Ripepe et al., 2001; Chouet et al., 2008), or to magmatic gas release (Jolly et al., 2010), although trigger mechanisms are still widely debated. Site and path effects, as well as the type of sensor deployed, significantly influence the waveform shape and its frequency recorded at a seismic station, meaning that the same event may be classified differently at two different stations. In addition, it is possible that differences in the location of the seismic event may influence the frequency content recorded. Consequently, it is essential to take these effects into full consideration when trying to classify volcano seismicity.

The further classification of seismic events into “families,” which all have a similar waveform shape as well as the same frequency content, allows the depiction of temporal and spatial changes in the source mechanism and the source location on a much smaller scale (e.g., Thelen et al., 2011; Salvage and Neuberg, 2016). By definition, families of seismic events should be generated by the same source mechanism and at the same source location (estimated at between one quarter and one tenth of the wavelength; Geller and Mueller, 1980; Neuberg et al., 2006) in order for the detected waveforms to have the same recorded shape at the seismometer (Minakami et al., 1951), as long as site and path effects on the seismic wave are minimal. Families of similar seismicity were first identified at Usu volcano, Japan, during a dome building eruption in 1944 by Minakami et al. (1951) and during the 1955 eruption of Bezymianny, Kamchatka (Gorshkov, 1959), and have since been identified at a number of active volcanoes around the world, including Redoubt, Alaska (e.g., Buurman et al., 2013); Mt. St. Helens, USA (Thelen et al., 2011); Colima, Mexico (Arámbula-Mendoza et al., 2011); Merapi, Indonesia (Budi-Santoso and Lesage, 2016); and Soufrière Hills, Montserrat (e.g., Rowe et al., 2004; Salvage and Neuberg, 2016).

Waveform similarity in terms of shape and duration can be evaluated by a cross correlation procedure where identical signals will result in a maximum cross correlation coefficient of 1 or –1, dependent upon their relative polarity and signals with no correlation resulting in a cross correlation coefficient of 0. The choice of similarity threshold (above which events are considered similar) is important: if it is too low there is a risk of placing events that are not similar into the same family; if it is too high similar events can be missed (Salvage et al., 2017). We define a “family” of events as events which show similar waveform shape characteristics, defined by having a cross correlation coefficient >0.7. This is in agreement with Green and Neuberg (2006), Thelen et al. (2011), and Salvage and Neuberg (2016), who suggest a cross correlation coefficient threshold of 0.7 in andesitic volcanic environments, since this is significantly above the correlation coefficient that can be produced from random correlations between noise and a waveform (Salvage et al., 2017). Here we first identify families of seismicity using a simple amplitude ratio algorithm, in addition to high waveform similarity on a multiple station network (≥0.7 cross correlation coefficient) using REDPy (Repeating Earthquake Detector, Hotovec-Ellis and Jeffries, 2016). Secondly, we take the core event of each family identified in the previous stage (an average of the stack of all events in that family aligned to the point of maximum cross correlation), in addition to all events identified during the study period that were identified using an amplitude ratio detection algorithm, and use these as a template to find more events within the continuous waveform data that may have previously been hidden by noise or simply not detected using a trigger algorithm [EQcorrscan, Calum Chamberlain, Victoria University of Wellington (Chamberlain et al., 2017)]. This template matching technique identifies a repeating event when the absolute cross correlation coefficient sum for a given template exceeds a threshold of 0.7 on each station, across a minimum of three stations. If repeating events occur within 0.5 s of one another, the strongest detection within this time period is taken, i.e., the event which exceeds the threshold the most. Both of these algorithms are written in Python and are open source, available through GitHub.

REDPy uses a simple amplitude ratio algorithm to detect seismic events from the continuous record on a multiple station network and then determines whether any events are similar by identifying all events over a given cross correlation threshold (0.7 in our case). This method, however, fails to recognize events hidden by noise, or events that are too closely spaced to re-trigger the algorithm (events within ~7.5 s of one another in this case). The core event is defined using OPTICS (Ankerst et al., 1999), a sorting-based clustering algorithm. The core event corresponds to the event with highest “reachability” within each family, meaning it correlates highly with many other members of that family. We took the core event of each family identified by REDPy (as well as all events which triggered the detection algorithm but appeared to have no similar events during the study period) and used these as input for EQcorrscan, which relies upon template matching rather than a trigger algorithm. Each core event was 10 s long in order to ensure that the entire waveform and coda were included and filtered between 0.5 and 15 Hz. In order for an event to be identified as similar as one of the core events, we again required a cross correlation coefficient >0.7. By this combined methodology, we found that the number of events identified in the dominant family (family containing the greatest number of events) went from 268 (using REDPy) to 771 events (using REDPy and EQcorrscan combined). EQcorrscan can be computationally intensive, however significantly more events were identified using this method (Table 1). The high discrepancy in number of events detected by these two different methodologies is a consequence of their detection algorithms. REDPy uses an amplitude ratio algorithm that only detects events above a given (user designed) threshold so events occurring within sections of high noise for example are not likely to be detected, and it can only detect one event per window length (in this case ~7.5 s). EQcorrscan is a template matching algorithm which therefore can detect event during noisy periods as it is simply focused on finding matching waveform patterns. In addition, the minimum inter-event time (IET) is 0.5 s, meaning it can in theory detect 7 times more events than REDPy within the same time window.

TABLE 1

Table 1. Details of fourteen families of similar seismicity identified in this study.

One family of similar seismicity that contained >10 events was identified in hindsight at Poás using REDPy (containing a total of 278 events). Furthermore, 84 events were detected using the amplitude ratio detection algorithm, but were not classified into families as no similarity was detected between these events. Many more events (a total of 771 within this family) were identified from the continuous record using the template matching technique of EQcorrscan (Table 1). Fourteen families, each containing >10 events per family, were identified using the combined methodology of REDPy and EQCorrScan (Figure 4). We cross correlated the core event for each family (Figure 4A) with every other core event and found that none of the families were similar to one another (the maximum cross correlation coefficient determined was 0.49 between two families of waveforms), suggesting that no single event is likely to belong to more than one family. LF families contained many more events and were active for longer periods of time (Figure 4B, green) than VT families, although all families showed distinct temporal patterns, in particular in relation to the rate of events with time (Figures 4E,F). Significantly, LF families began much earlier in the sequence (around the end of March) compared to the VT families that were active at the time of the large phreatomagmatic event on 23 April. The dominant family of similar waveforms showed two phases of heightened activity (Figure 4C): (1) an increase in daily events from 28 March to 2 April 2017, followed by a steady decline; and (2) a smaller, apparently less significant increase in event rate in the days prior to the large phreatomagmatic eruption on 23 April. A similar pattern of two phases of activity was observed for all other low frequency families identified. VT families showed a single period of heightened activity, starting in early April (Figure 4D). Following the beginning of daily phreatic eruptions at Poás (12 April onwards), the dominant LF family registered a significant increase in the RMS (Root Mean Square) amplitude of events (e.g., Figure 5A, following the blue dashed line), which decreased again on 25 April, 2 days after the large phreatomagmatic eruption. Significantly, the dominant LF family also showed a clear decrease in FI in the hours prior to the large phreatomagmatic eruption on 23 April (Figure 5A). The Frequency Index (FI), developed by Buurman and West (2010), is a proxy for the spectral content of the waveform, and is based upon the ratio of energy in low and high frequency windows, with a base-ten logarithm in order to reduce the index to a single value. A negative FI means the waveform is dominated by low frequency energy; a positive FI demonstrates a majority of energy in the high frequency band; and a FI of zero suggests that the waveform has equal amounts of high and low frequency energy. Here, we define a low frequency window between 0.5 and 5 Hz, and a high frequency window between 5.1 and 15 Hz. A decrease in FI would therefore suggest that the events began to contain a higher proportion of low frequency energy, potentially related to the increase ability for fluid to move through the volcanic system with ease, as the fracture network (responsible for the generation of higher frequency events) is thoroughly developed. This decrease in FI occurs coincidently with an increase in RMS amplitude, suggesting that the lower frequency events are larger in amplitude. Within the dominant LF family, the IET appears to evolve slowly with time, as events become more sporadic (Figure 5A). The minimum IET detected within the dominant LF family was 1.18 s, although ≥95% of the IETs detected were greater than the event duration of 10 s. Furthermore, in the 12 h prior to the large phreatomagmatic event on 23 April, no events from the dominant LF family were identified (Figure 5A). Since the dominant VT family continues to register a small number of events during this period (Figure 5B), this could be interpreted as either a sudden shut off in the conditions necessary to generate these types of events, or the possibility of aseismic magma movement (e.g., Neuberg et al., 2006). VT families appear to show no significant patterns in RMS amplitudes of events, IETs or changes in FI over the investigated study period (Figure 5B).

FIGURE 4

Figure 4. (A) Core events (normalized) of each family identified by REDPy. VT events are shown in gray, LF events in green. (B) Duration of families, colored with respect to (A). (C) Daily count of dominant LF family (greatest number of events) corresponding to the first waveform in (A) identified by template matching technique (EQcorrscan). (D) Daily counts of dominant VT family (greatest number of events) identified by template matching technique (EQcorrscan). (E) Cumulative counts of earthquakes within each LF family. (F) Cumulative counts of earthquakes within VT families. The dashed black line indicates the timing of increased VT activity counted by OVSICORI-UNA personnel (26 March), the dashed blue line indicates the beginning of daily phreatic eruptions (12 April), and the dashed red line indicates the timing of the large phreatomagmatic eruption (23 April).

FIGURE 5

Figure 5. Example of family analysis conducted on the dominant LF family [30-03-2017 08:46, (A)] and the dominant VT family [12-04-2017 05:49, (B)]. The dashed black line indicates the timing of increased VT activity counted by OVSICORI-UNA personnel (26 March), the dashed blue line indicates the beginning of daily phreatic eruptions (12 April), and the dashed red line indicates the timing of the large phreatomagmatic eruption (23 April).

4. Hindsight Forecasting of 22 April 2017 Eruption

The ability to forecast the timing, intensity and type of volcanic activity is one of the key issues facing volcanologists today. Since the time series analysis of families at Poás in April 2017 identified an accelerating trend within the dominant LF family, the Failure Forecast Method (FFM) may be applicable for identifying a forecasted timing of eruption. The FFM is based on an empirical power-law relationship relating the acceleration of a precursor (d²Ω/dt²) to the rate of that precursor (dΩ/dt) (Voight, 1988) by:

\begin{array}{rcl} \frac{d^{2} Ω}{d t^{2}} = K {(\frac{d Ω}{d t})}^{α} & (1) \end{array}

where K and α are empirical constants. Ω can represent a number of different geophysical precursors, for example low frequency seismic event rate (Salvage and Neuberg, 2016), event rate of all recorded seismicity (Kilburn and Voight, 1998), or the amplitude of seismic events (Ortiz et al., 2003). The parameter α can range between 1 and 2 in volcanic environments (Voight, 1988), or may even evolve from 1 toward 2 as seismicity proceeds (Kilburn, 2003). An infinite dΩ/dt suggests an uncontrolled rate of change and here is associated with an impending eruption. The inverse form of dΩ/dt is linear if α=2, and therefore in this case the solution for the timing of failure is a linear regression of inverse rate against time, with the timing of failure relating to the point where the linear regression intersects the x-axis in graphical form (Voight, 1988). As a deterministic approach for forecasting the timing of volcanic eruptions, the FFM relies upon several assumptions, including that the acceleration can be described by a simple power law and that the time of the eruption is related to the time at which this power law reaches a singularity (Boué et al., 2016). The FFM has proved useful at accurately forecasting volcanic eruptions in both near real time (e.g., Cornelius and Voight, 1994) and in hindsight evaluation (e.g., Budi-Santoso et al., 2013; Salvage and Neuberg, 2016). However, Boué et al. (2016) have suggested from a study of data over a 13 year period at Volcán Colima, Mexico, and a 10 year period at Piton de la Fournaise, Reunion, that not all cases of accelerating seismicity are suitable for the analysis by the FFM, in particular if the acceleration does not follow a single power law increase. In fact, only 36% of eruptions at these volcanoes could be forecasted in real-time, although hindsight forecasting fared better, with ~50% of eruptions being successfully forecast when utilizing the entire precursory seismic sequence in the forecasting model.

The FFM was first developed as a tool for forecasting the timing of slope failure using accelerating material creep; a cause and consequence of one single active system generating failure (Fukuzono, 1985). However, a volcanic system is inherently more complex, and accelerating magma ascent could be detected at several positions in the magma plumbing system with different phase delays and amplitudes. It may be advantageous to use the FFM in collaboration with a single family of seismicity since a single family originates from the same source location and is produced by the same source mechanism (e.g., Geller and Mueller, 1980; Petersen, 2007; Thelen et al., 2011), which may produce a more accurate forecast, potentially a consequence of focussing upon a single active system at depth without interference from other sources of error (Salvage and Neuberg, 2016). In other cases, using only a single family in conjunction with the FFM produces a less useful forecast than when using all the identified seismicity (e.g., Tungurahua, Ecuador Bell et al., 2017).

Here, we set α=2, to allow for simple implementation, and because it has proved to be a good choice for the value of α at a number of volcanoes (e.g., Cornelius and Voight, 1994; Chardot et al., 2015; Salvage and Neuberg, 2016). We only consider data from 18 April onward in this analysis and consider the acceleration over a 96 h period from 18 April until 22 April. The event count is dependent upon events not occurring within 0.5 s of one another (a user defined minimum threshold in EQcorrscan), however, in some volcanic environments individual events, in particular low frequency seismicity, may merge into tremor as their IET decreases (Neuberg et al., 2000). Consequently, the event count recorded in our case will not include tremor episodes if individual events cannot be constrained, despite this signal being associated with precursory activity at some volcanoes e.g., White Island, New Zealand (Chardot et al., 2015). We use the R² value of the least-squares linear regression as an initial indication of the confidence of the forecast made. The closer the R² value is to one, the more confidence that can be placed in the forecast, although in this case we use it simply as an indicator of the fit of the regression before performing more detailed residual analysis. R² <0.65 are considered to represent a poor relationship between the observed data and the fitted FFM model (Barrett, 1974). Consequently, we define a successful forecast as one where the timing of the eruption is within 3 h of the known timing of the eruption, with an R² > 0.65.

Accelerations in the number of seismic events per hour prior to the large phreatomagmatic eruption on 23 April 2017 at 04:12 UTC (Figure 6) were only identified when using all families combined, and in the dominant LF family. No clear accelerations could be identified in VT families. We use the number of events per hour as an indicator of the activity, rather than events per day, as we consider this to be more useful if the process is to move into real time: the number of events occurring within a time period of <1 h would be difficult to process and understand quickly enough for decision makers; and the number of events per day may be too long between processing times to generate usable forecasts. The acceleration identified when using all families combined (Figure 6A) is subtle, and in fact a clear acceleration followed by a deceleration can be seen from 18 to 19 April, which contains higher event counts. Application of the FFM to the entire accelerating sequence produces a very poor forecast, with R² values much lower than what is deemed confident (R²≪0.65, Figure 6B), despite the forecasted time of eruption occurring within 5 h of the known timing of the eruption. Figure 6C shows the acceleration of seismic events per hour for the dominant LF family (30-03-2017 08:46). An acceleration can be identified from the 19 April onward (Figure 6C). The acceleration of this family appears to stop approximately 24 h before the large phreatomagmatic eruption, producing a period of quiescence where no repeating events are identified. Application of the FFM (Figure 6D) indicates a poor forecast, with an R² value of 0.13. The least-squares linear regression generates a forecast for 22 April at 06:48, approximately 22 h before the known timing of the eruption.

FIGURE 6

Figure 6. Hourly counts of seismicity and application of the FFM prior to a large phreatomagmatic eruption at Poás. Each data point represents either: the hourly count of seismic events in that family (A,C) or the inverse hourly event rate (B,D). The vertical red line represents the known timing of the phreatomagmatic eruption on 23 April at 04:12 UTC. Accelerating trends, and their corresponding least-squares linear regressions are shown by dotted lines. (A,B) represent the hourly counts from all families of seismicity identified during this time period; (C,D) represents data from the dominant LF family of events (30-03-207 08:46).

Both generated forecasts were extremely poor, with R² ≪0.65 and only one forecast was made within 5 h of the known eruption time, even though an observable acceleration in seismic event rate could be identified. Following Chardot et al. (2015), residuals were calculated through time in order to test the assumptions of using a least-squares linear regression with this data, and plotted as a histogram, where the number of bins is equal to 2n^1/3 (n is the number of samples), as according to the Rice rule, to verify the normality of the residual distribution (Figure 7). Residuals were defined as the difference between the observed hourly event count for each family, and the best fit model to this data. Using the residual analysis, the least-squares linear regression model can be deemed appropriate if: (1) the residuals do not follow a trend; (2) the residuals do not increase or decrease as a function of time; and (3) the residual distribution follows a Gaussian distribution. Our residual analysis (Figure 7) suggests that the error structure of the data is inconsistent with a least-squares linear regression model when α is equal to 2, since the distribution is not Gaussian. This has been suggested previously by Bell et al. (2011) for seismic event rates. In addition, all of our residuals appeared to follow a trend, which suggests that a more complex model is needed to define the data, and residuals increase as a function of time, indicating they do not exhibit equal variance. Consequently, we can conclude that the least-squares linear regression applied in this instance is not the most appropriate model for describing the accelerating behavior observed.

FIGURE 7

Figure 7. Residual analysis for least squares linear regressions presented in Figure 6. The y-axis of (A) and (C) depict the count. (A,B) Residual analysis for all families of seismicity identified. (C,D) Dominant LF family (30-03-2017 08:46). The vertical red line indicates the mean of the calculated residual values (A,C).

5. Discussion

5.1. Characteristics of Families Identified

The identification of families of similar seismicity at Poás (events with a cross correlation coefficient >0.7) suggests: (1) a stable source process (non-destructive); (2) that the trigger mechanism of such events must be able to be recharged quickly (in this case in <2 s since this was the minimum IET observed); and (3) that these conditions must occur consistently at the same location (e.g., Green and Neuberg, 2006; Petersen, 2007). In particular, the identification of a number of families that are simultaneously active suggests that a number of distinguishable sources must be active at the same time beneath Poás. Furthermore, the identification of both LF and VT families acting simultaneously suggests a diversity in the ongoing physical processes at depth, once path effects have been accounted for in the classification. As discussed earlier, the source mechanisms for LF and VT seismicity are often disputed, but here we suggest that both VT and LF occurrence is related to the movement of magmatic fluid and gases toward the surface, with VT families suggesting the generation of fractures and the potential opening of new magmatic pathways due to increased stresses at depth as a result of the presence of magmatic fluid (e.g., Lahr et al., 1994; Kilburn, 2003), and LF families potentially suggesting the movement of this fluid through these fractures as a result of an increased strain rate in the magma (e.g., Chouet, 1988; Neuberg et al., 2000; Thomas and Neuberg, 2012). Since the VT families contained fewer events than families containing low frequency events (Table 1), we suggest that the dominant ongoing process during this time (once a new fracture pathway had been opened, demonstrated by the swarm of the VTs in March 2017) was magmatic fluid movement, most likely, a mixture of magma and gases. This is supported by the significant increase in RMS amplitude of LF events (in particular after 12 April, Figure 5A) in comparison to VT events (Figure 5B), suggesting that the movement of fluid dominated the processes occurring at depth during April 2017. We note that this method does not allow for the characterization of volcanic tremor, but instead will either count individual events within a tremor episode if the IET is >0.5 s instead of a single tremor sequence, or miss the tremor episode entirely, affecting the final event count. Although there are instances where volcanic tremor has appeared as a precursor to eruptive episodes at some volcanoes (e.g., Chardot et al., 2015), here we are interested in families of repeating seismicity as a precursor, and not other signals. Further analysis is required to determine whether volcanic tremor played a significant role as a precursor to this eruptive event.

All identified families of low frequency events (including the dominant family) appear to return to fewer seismic counts immediately prior to the eruption after an acceleration (Figure 6). This could represent a decrease in magmatic flow rates due to a physical obstruction, the generation of a damage zone that is responsible for impeding fluid flow, or the intermittent advance of magmatic fluid through a fracture network. A sudden drop in seismic event rate could also represent a change in the system to more aseismic magma movement, as fluid pathways become fully open, meaning less seismic events are generated. In our case, since the seismicity in this instance can be classified as low frequency, either the physical conditions necessary for the generation of this family are gradually changing (e.g., Stephens and Chouet, 2001); a change in the fluid flow rate leading to the generation of systematically fewer seismic events is occurring; or the movement of magmatic fluid becomes aseismic. Hotovec et al. (2013) also noted a period of quiescence prior to explosive events at Redoubt volcano, Alaska during an eruptive period in 2009, although on the timescale of seconds, rather than hours. Rodgers et al. (2015) noted a systematic decrease in LF seismicity at the same time as an increase in VT seismicity, over a number of months prior to eruptive activity at Telica volcano, Nicaragua from 2010 to 2011. They attribute this changing seismicity to an increase in pressurization (indicated by the VT seismicity) and thus a sealing of an active hydrothermal system (indicated by the drop in LF seismicity), suggesting it may drive phreatic eruptions. At Poás this scenario appears unlikely since we see no increase in VT seismicity with a decrease in LF events, and the reappearance of LF seismicity at the time of the eruption suggests that this 12 h of quiescence does not mark the timing of a sealing of the hydrothermal system, as it is unlikely to occur on this timescale. A significant decrease in the cross correlation coefficient immediate prior to the period of quiescence suggests a small, but significant change in the conditions needed to generate the dominant LF events, which become less similar with time (Figure 5A). Petrological evidence suggests that changes in fluid flow rate in response to obstructions in the conduit may be a more plausible explanation. We interpret the increasing juvenile content of the ash sample over time as being related to a fresh batch of rising magma “cleaning” out a path to the surface. Initially the path is obstructed by older, altered material (such as a plug), which must first be ejected before fresh magmatic material can reach the surface. Consequently, earlier in the eruptive period, the ash content is dominated by interactions of the hydrothermal system with the rising magma, and the opening of a pathway to the surface. Seismic events within the dominant family are no longer recorded in the final hours prior to the eruption since a clear pathway to the surface has been generated, and therefore fluid movement becomes aseismic.

5.2. Forecasting Potential

At Soufrière Hills volcano, Montserrat, the use of similar seismicity in collaboration with the FFM appears to allow the generation of a more accurate forecast, since isolating a single system at depth avoids additional uncertainties introduced by averaging data over a number of different accelerating phenomena, and therefore reduces the misfit between the data and the forecast (Salvage and Neuberg, 2016). At Poás, the identification of similar seismicity and its use as a forecasting tool is promising, since our hindsight analysis highlights that identifying repeating seismicity allows a much more detailed interpretation of ongoing processes at depth, in addition to the identification of far greater numbers of events for analysis. Using the methodology described in this paper, it would have been possible to detect accelerations in families of similar seismicity in near real time, although the analysis of this data is unlikely to have produced confident forecasts. Firstly, R² values for the least-squares linear regressions identified were always considered to be unacceptable for a confident forecast (R²≪ 0.65). Secondly, a number of accelerations in these families of seismicity could have been identified earlier in April: some prior to eruptive events (e.g., 12 April), and others occurring with no associated surface manifestations. Thirdly, the identification of decelerating trends, or a sudden drop in the seismic event counts may have falsely suggested a decline in activity at Poás, and therefore reduced the likelihood of a confident forecast. Lastly, accelerations were not identified in all families. For example, no clear accelerations could be identified within the families of VT seismicity. These factors become particularly important when trying to forecast eruptive events in real-time, and therefore more research is required to determine how some or all of these issues can be accounted for in a forecasting model.

Previously, Bell et al. (2011) suggested that it is inappropriate to use the FFM with a least-squares linear regression, since this does not account for the correct error structure of earthquake count data. Our residual analysis for the least-squares linear regressions applied when α is set to 2 supports this, since the residual error structures are not Gaussian (Figure 7). They therefore suggest that using a Generalized Linear Model (GLM) and a Poissonian error distribution where α is equal to 1 may be more appropriate. However, using the FFM in conjunction with the GLM also violates certain assumptions which are associated with the volcanic environment. The GLM suggests that the system being modeled is memoryless and that past events do not affect the future. We agree with Hammer and Ohrnberger (2012) and suggest that this is not an appropriate assumption for seismicity occurring in volcanic settings. Calculation of α (e.g Boué et al., 2016) is an essential next step in the characterization of the accelerations observed at Poás. An alternative model which may be more appropriate is a maximum-likelihood methodology which utilizes observed data to determine a model result with the greatest probability to maximize the likelihood function (e.g., Bell et al., 2013). This methodology also does not require the binning of data into time windows, which is advantageous. Although there are many issues with using the FFM in conjunction with a least-squares linear regression analysis (setting α to 2), in particular in real time scenarios, other regression models violate different assumptions of the FFM, indicating the inherent complexity of this problem. Our research indicates that a least-squares linear regression is not an appropriate tool to use for forecasting this eruptive event at Poás, consequently eliminating one of the methods that is currently popularly used. We consider these “negative” results a contribution to our understanding of the scenarios that can be forecast using this methodology and hope to expand further on this analysis by implementing different regression models in the future to this data set.

5.3. Conceptual Model: Phreatomagmatic Explosive Eruptions at Poás

Investigations into the volcanic structure of Poás in the 1980s suggest that a high density cylindrical plug, approximately 1,000 m in radius, sits between 500 and 800 m beneath the crater floor (Rymer and Brown, 1986; Casertano et al., 1987), potentially connected to the surface by a (now solidified) vertical intrusion beneath the dome in the active crater (Fournier et al., 2004). Fischer et al. (2015) suggested that phreatic eruptions at Poás are therefore caused by a gas pressure build up beneath this sealed plug (which provides the surface crater lake with heat and volatiles and is likely a chilled margin resulting from the crystallization of an older magmatic body), which is catastrophically released through hydrofracturing. De Moor et al. (2016) showed that phreatic eruptions at Poás were accompanied by short-term increases in SO₂/CO₂ and higher SO₂ fluxes, again suggesting that high temperature magmatic gas injection drives phreatic eruptions. A similar mechanism for the onset of phreatic eruptions has been suggested for Mt. Ruapehu, New Zealand (Christenson et al., 2010), where sulfur within the system creates an impermeable plug within the volcanic conduit and consequently leads to the accumulation of gases beneath it, the elevation of pore pressures, and the sudden catastrophic release when critical pressures are reached.

The earliest identified indication of renewed activity at Poás in 2017 was a swarm of VT (high frequency) seismicity on 26 March (Figure 3B). High frequency seismicity is associated with the generation of new fractures when critical stresses are reached, that potentially allows fluid movement (e.g., Kilburn, 2003). This VT swarm is likely to represent the brittle fracturing of a previously sealed chilled margin or hydrothermal seal, either as a new batch of magmatic fluid beneath it begins to push upwards, or as a result of partitioning volatiles into residual melt during the crystallization process of the upper portion of the cooling magma body. It is unclear, however, as to whether the magma itself is forcing its way up toward the surface due to volatile decompression, or whether the fracture (generated as a result of increased pressurization from the nearby magmatic fluid and gases) creates a vaccuum as gas escapes out its top, allowing the magma beneath to be “sucked” upwards, as the magma below decompresses and migrates toward the surface. Increased fracturing of the seal may also result from the increased local strain rate around the intruding magma body, forcing the surrounding rock into failure (e.g., Fournier, 1999). In hindsight, the counts of LF events in February and early March 2017 (Figure 3A) may be an early indication of movement of magmatic fluid to beneath the seal allowing the build up of pressure, although this was not noted at the time as the event rate was not considered unusual.

In July 1980, an increase in VT seismicity was first observed before an increase in low frequency events at Poás. Although no large explosive eruption occurred associated with this seismicity, a significant increase in gas emissions was noted, and was believed to be caused by crystallization-induced degassing beneath a water saturated chilled margin caused by an intrusive episode (Casertano et al., 1987; Rowe et al., 1992; Rymer et al., 2000). In March 2017, the swarm of VTs was short lived and contained relatively few events (a total of 22 earthquakes) unlike the swarm in July 1980, that contained hundreds of events (Casertano et al., 1987). This may be evidence of the “Kaiser effect,” in which fracturing and seismicity cannot be generated unless the previous maximum stresses of the system are exceeded (Fredrich and Wong, 1986; Smith et al., 2009). Consequently, it would be expected that after the resealing of the chilled margin as volcanic activity wanes, greater stresses are needed in order to re-fracture that body in a new eruptive episode, even decades later. Tuffen et al. (2003) has suggested that the rehealing of small veins through which fluid and gases can travel in an andesitic volcanic environment may be able to occur on timescales of minutes to hours, meaning that Poás would have had the opportunity to re-seal itself since the last activity in 1980. Other explanations could include a more ductile environment at depth in 2017, and/or changes to the magnitude distribution of events between the two time periods.

An increase in pressurization beneath a seal at depth is supported by GPS measurements from early March 2017, which registered inflation of the edifice (Figure 3C), suggesting that the VT swarm in late March represents a time when pressurization at depth (induced by a seal) reached a critical level and allowed fracturing to occur. Hypocentre locations may help to define the depth of this seal, but this is beyond the scope of this paper. Prior to the 23 April 2017 eruption, only 2 GPS stations were recording, leading to large uncertainties in the location of a pressurization source at depth. Simply determining the intersection of deformation vectors suggests the pressurization source lies between 1 km and 6 km beneath the volcano. Casertano et al. (1987) and Rymer and Brown (1986) have previously suggested that a dense plug, which may act as a seal, sits between 500 and 800 m beneath the crater surface, in agreement with our observations for a pressurization source. The accelerated inflation of the crater in April 2017 suggests magma reaching shallower levels within the crust as it moves through a (now generated) fracture network toward the surface, likely behind the seal through which fluid migration is still impeded, but not impossible. The movement of magma to shallow levels is also supported by the increased SO₂ flux at the surface, increases in the SO₂/CO₂ ratio, increases in lake temperatures, as well as a significant increase in low frequency seismicity from 22 March onward (Figure 3). In particular, the significant change in the gas composition from 29 March 2017 onward suggests the injection of gas to shallower levels, probably facilitated by the developing fracture network on 26 March.

The identification of a number of families of seismicity indicates that the movement of fluid (LF families) and the generation of new fractures (VT families) occurred consistently in the same location and by the same source mechanisms, and occurred simultaneously during the precursory period. Therefore, it is likely that developing fracture networks occurred in a limited number of locations, and since low frequency events occurred more frequently and with greater amplitudes and energy release, that the movement of magmatic fluid was the principal ongoing process at this time. Similarly, from 1978 to 1990, counts of LF events were considerably higher than those for VT events, during a variety of phreatic and protoplasmatic activity at Poás(Martínez Cruz, 2008). Increased pressurization below the partially sealed chilled margin due to the build up of magmatic fluid and the degassing of this magmatic fluid in shallower reservoirs is likely to lead to the catastrophic failure of the seal, and therefore to phreatomagmatic eruption. As new magmatic fluid moved toward the surface through the newly fractured seal, it may have picked up some of the surrounding altered conduit material and breccia from the generation of fractures, meaning the first erupted material in early April contained lower percentages of juvenile material than later eruptions, as the conduit was not yet fully open. As the magmatic fluid reached very shallow levels, it came into contact with the active hydrothermal system at Poás including the crater lake at the surface, which is likely to further enhance explosive activity.

We therefore suggest that the large phreatomagmatic eruption at Poás on 23 April 2017 was the result of a fresh batch of magmatic fluid becoming initially stalled behind a sealed chilled margin or hydrothermal seal at approximately 1 km depth. Once pressures were critical, some fracturing of this seal occurred (VT swarm in March) allowing the movement of magmatic fluid, and in particular volatiles, to shallower levels (noted by the increase in low frequency seismicity, gas fluxes, inflation of the crater area and significant changes in the gas composition). As the magmatic fluid migrated through the system to shallower levels, rapid degassing occurred, and it picked up surrounding edifice material and material form the seal. When it came into contact with the active hydrothermal system and crater lake, an even more energetic explosive eruption was generated.

6. Conclusions

Phreatic and phreatomagmatic eruptions at volcanoes are poorly understood, difficult to forecast, and often pose a serious threat to local populations and tourists visiting volcanic areas. A large phreatomagmatic eruption occurred at Poás volcano, Costa Rica, on 23 April 2017 at 04:12 UTC, following approximately one month of unrest (rapid inflation, increased seismicity, increased gas fluxes and smaller phreatic eruptions), potentially as a result of pressure build up beneath a partially sealed, chilled margin or hydrothermal seal at depth. The sudden fracturing of this seal as critical pressures were reached triggered an explosive eruption, which was further enhanced by interaction of the magmatic fluid with the active hydrothermal system at Poás, and the crater lake at the surface.

Fourteen families of similar seismicity were identified in the days prior to this eruption, which included both LF and VT families. We suggest that the VT families are indicative of further fracturing of the volcanic edifice and the seal, creating new fluid pathways, and that the LF seismicity reflects the movement of magmatic fluid through these pathways toward the surface. Detailed analysis of these families suggests that only the dominant LF family produced an accelerating trend in the hours prior to the catastrophic eruption on 23 April, which we have used in collaboration with the FFM to forecast the timing of the event. When using all families of seismicity to forecast the timing of the eruption on 23 April, an accelerating trend was also identifiable, although it did not produce an accurate forecast. The dominant LF family of events showed an accelerating trend which produced a forecast approximately 22 h from the known timing of the eruptive event. Analysis of each LF family individually suggested a cessation of events in the 12 h prior to the large phreatomagmatic eruption, which was most pronounced within the dominant family of LF events. This deceleration could have been misinterpreted in real time to signify a slowing of the activity at Poás, but we suggest that at this time fluid movement through the system became aseismic. No confident forecasts were generated, despite an obvious acceleration in seismic event rate. The lack of acceptable forecasts may result from the use of a least-squares linear regression with the FFM, which based on residual analysis of this data is not an appropriate model to use. Our analysis allows us to eliminate some common parameters which are commonly used for forecasting volcanic eruptions, to search for other more complex models to explain the accelerating seismicity in this scenario. Furthermore, the use of similar seismicity rather than simply defining seismic events according to their frequency content alone allows a more detailed analysis of time series trends to be carried out. For example, this methodology allowed us to identify far greater numbers of events for analysis, and identified that a number of distinct specific sources were generating seismicity at depth, as well as significant changes in the frequency content of waveforms with time. Further investigation is required to determine whether all large phreatomagmatic eruptions at Poás are preceded by accelerating families of seismicity and consequently whether this can be successfully used as a forecasting tool for future events in real time, if more suitable regression analyses can be determined for use at this volcano.

Data Availability

All monitoring data for the period of interest, and all data generated in analysis for this manuscript is available upon request, without undue reservation, to any qualified researcher. Please contact Rebecca O. Salvage with requests.

Author Contributions

RS identified the seismic families, computed the family analysis and accelerations, calculated the forecasts using the FFM and is responsible for writing the manuscript. GA provided key information regarding the chronology of the eruption in April, as well as petrological analysis of ash samples and lake water samples. JdM conducted gas analysis, JP and JB were responsible for monitoring of seismicity and seismic counts. MC conducted SEM and microprobe analysis of ash and other collected samples. CM conducted deformation (GPS) analysis. MMC was responsible for the lake and boiling springs geochemistry analysis and SEM-EDS and ICP-MS analyses of some of the 2017 bombs of Poás. All authors provided critical feedback on the manuscript and the ideas which are formulated here, and approved the final version.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank the entire team at OVSICORI-UNA, but in particular the field technicians and volcanological personnel who tirelessly maintain the monitoring network for all the volcanoes in Costa Rica, including Poás. RS would like to particularly thank Lauriane Chardot for inviting her to contribute to this special issue; and Calum Chamberlain and Alicia Hotovec-Ellis for answering numerous questions related to their respective codes. JdM thanks the Deep Carbon Observatory (DCO), Istituto Nazionale di Geofisica e Vulcanologia (INGV), and the United States Geological Survey Volcano Disaster Assistance Program (USGS-VDAP) for contributions to gas monitoring efforts at Poás. We thank the United States Geological Survey Volcano Disaster Assistance Program (USGS-VDAP) for providing seismic equipment to replace damaged seismic stations following the eruption. Finally, we thank W. McCausland, Andy Bell, and an anonymous reviewer as well as Corentin Caudron, the editor, for helpful comments which greatly improved the manuscript.

References

Alvarado-Induni, G. E. (2005). Costa Rica: land of volcanoes. San José: EUNED.

Ankerst, M., Breunig, M. M., Kriegel, H.-P., and Sander, J. (1999). Optics: ordering points to identify the clustering structure. ACM Sigmod Rec. 28, 49–60. doi: 10.1145/304181.304187