Edinburgh Research Explorer Statistical evidence of transitioning open-vent activity towards a paroxysmal period at Volcán Santiaguito (Guatemala) during 2014–2018

Long-term of oferuption, including a period of elevatedexplosiveactivity in2016. Here,we usea new catalogue of explosions to characterise changes in the eruptive regime during the study period. We identify four different phases of ac- tivity basedon changes inthefrequencyandmagnitudeofexplosions. Atthe two ends ofthespectrumofrepose timeswe ﬁ ndpairsofexplosions withnear-identicalseismic andacousticwaveforms,recordedwithin 1 – 10min of one another, and larger explosions with recurrence times on the order of days to weeks. The magnitude-frequency relationship for explosions at Santiaguito is well described by a power-law; we show that changes in b -value between eruptive regimes re ﬂ ect temporal and spatial changes in rupture mechanisms, likely controlled by variable magma properties. We also demonstrate that the distribution of inter-explosion repose times between and within phases is well represented by a Poissonian process. The Poissonian distribution describing repose times changes between and within phases as the source dynamics evolve. We ﬁ nd that changes in source properties restrict the extrapolation of explosive behaviour to within a given eruptive phase, limiting the potential for long-term assessments of anticipated eruptive behaviour at Santiaguito. © 2020 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http:// creativecommons.org/licenses/by/4.0/). Normal ’ explosions to at with the period of heightened activity erratic recurrence We analyse this characteristic across the eruptive phase transitions


Introduction
The Santiaguito lava dome complex, situated in the Western Highlands of Guatemala, has been continually active since 1922 (Harris et al., 2003). The complex consists of four domes, Caliente, La Mitad, El Monje and El Brujo, which formed along an East-West trending fracture at the southern foot of the collapse amphitheatre formed by the 1902 plinian eruption of Santa Maria Volcano (Rose, 1972). Since 1975 activity at Santiaguito has been centred at Caliente, mostly characterised by effusion of blocky lava flows and small-to-moderate gas-rich explosions (Harris et al., 2003). The Instituto Nacional de Sismologia, Vulcanologia, Meteorologia e Hidrologia (INSIVUMEH), responsible for monitoring the activity at Santiaguito, issues regular activity reports including information on the rates of occurrence of explosions and the height of the associated plumes. During 'normal' eruptive behaviour, as revealed by decades of monitoring, explosion plumes have been observed to reach heights between 0.5 and 1 km above the active vent (e.g. Global Due to its protracted eruption, which has lasted for nearly 100 years, and the unique vantage point offered by the summit of Santa Maria overlooking the active caliente vent, the explosions at Santiaguito have been the focus of many multi-disciplinary investigations. Previous studies have included the use of optical imagery, thermographic cameras (e.g. Sahetapy-Engel et al., 2004;Johnson et al., 2004;Sahetapy-Engel et al., 2008;De Angelis et al., 2016;Lamb et al., 2019), tiltmeter (e.g. Johnson et al., 2004;Johnson et al., 2014;Lavallée et al., 2015), infrasound (e.g. Sahetapy-Engel et al., 2008;Jones and Johnson, 2011;Scharff et al., 2014;De Angelis et al., 2016;Lamb et al., 2019), seismic (e.g. Johnson et al., 2004;Sahetapy-Engel et al., 2008;Sanderson et al., 2010;Johnson et al., 2014;Scharff et al., 2014;Lavallée et al., 2015;Lamb et al., 2019;Hornby et al., 2019), doppler radar (Scharff et al., 2014) and UV imaging (e.g. Holland et al., 2011;Esse et al., 2018). Bluth and Rose (2004) proposed that the primary mechanism of the small-to-moderate gas-and-ash explosions at Caliente is shearinduced fragmentation of dacitic magma near the margins of the conduit. This shearing is induced by the increased flow velocity brought on by closed system degassing as the ascending magma undergoes decompression, causing the magma viscosity to increase and prompting brittle deformation at the high strain rate experienced in the shallow conduit (Holland et al., 2011). The shearing events are thought to generate many small connected cracks (Rhodes et al., 2018), which are held open during release of the gas-and-ash mixture (cf. Kendrick et al., 2016). After the explosion the fracture pathways shut and begin to heal, resetting the system, leading to the cyclic nature of the explosive activity observed at Santiaguito (Holland et al., 2011). This mechanism is supported by the observations that 1) volcanic ash and gas are released along active fractures which can be seen on the dome surface at the onset of the explosions, 2) volcanic ash associated with these small explosions shows textural disequilibrium resulting from rapid frictional heating (i.e., unhomogenised melt schlieren containing fresh vesicles) associated with faulting events, and 3) more pronounced tilt signals and very long period (VLP) seismicity occur as a result of increased shear traction when ash is generated . The explosion source has been constrained to depths between 100 and 620 m below the vent through the comparison of acoustic and seismic signal onsets (Sahetapy-Engel et al., 2008), between depths of 100 and 600 m though the analysis of rock samples taken from the lava flow units (Scott et al., 2012), and at about 300 m depth below the vent via Mogi source modelling of VLP signals (Johnson et al., 2014). Yet, these models proposed for the small-to-moderate gas-and-ash explosions do not apply to the large explosions during the period of heightened activity of 2016, which were sufficiently powerful to excavate the lava dome . Protracted volcano monitoring provides us with datasets, which may be scrutinised to understand eruptive behaviour with the aim to find ways to predict recurrence in activity (e.g., Papale, 2018). The monitored datasets of open-vent volcanic systems are different to those of volcanoes ending a period of quiescence, where unrest can indicate an impending eruption. For open-vent systems such as Santiaguito, where activity has persisted for nearly 100 years, changes in the eruption style may be signalled in a more subtle way (Sparks, 2003). Therefore, understanding the day-to-day 'baseline' behaviour and associated trends and statistical attributes of geophysical signals during different phases of eruption is crucial in efforts to provide accurate forecasting of paroxysmal activity, which pose great hazards to the surrounding areas. In systems that produce many self-similar events which display a range of magnitudes, magnitude-frequency analysis can be used to determine the relationship between small and large events via the socalled seismic b-value. Magnitude-frequency analysis, most commonly used in seismic studies (e.g. Schorlemmer et al., 2005;Chen et al., 2006;Tsukakoshi and Shimazaki, 2008;Farrell et al., 2009;El-Isa and Eaton, 2014;Huang and Beroza, 2015), often refers to the Gutenberg-Richter (GR) relationship (Gutenberg and Richter, 1944) and Ishimoto-Iida's formula (Ishimoto and Iida, 1939), but has also been applied to explosive volcanic eruptions (e.g. Deligne et al., 2010;Nishimura et al., 2016;Sheldrake and Caricchi, 2017;Rougier et al., 2018). These studies have characterised the magnitude-frequency relationships for global volcanism as well as for individual volcanoes, and have shown how these distributions carry information on the explosion sources and eruption style. Repose time analysis has also become an important tool at open-vent systems (e.g. Connor et al., 2003;Watt et al., 2007;Lamb et al., 2014;Varley et al., 2018) in establishing relationships between source processes. Changes in the repose time between explosions have been used as an indicator of changes within volcanic systems, and differences in repose times between events has been used to show differences in the physics controlling eruptive processes (Varley et al., 2018).
Between November 2014 and December 2018, the University of Liverpool (UK) and the Karlsruhe Institute of Technology (Germany) deployed a permanent seismo-acoustic network around Santiaguito to continuously record the activity of the volcano. During this period multiple phases of activity were observed including effusive behaviour interspersed with small-to-moderate explosions, characteristic of the normal baseline behaviour as defined by Harris et al. (2003), and a period of heightened activity defined by large gas rich explosions during 2016 Wallace et al., 2020). In the initial 12 months following our installation at Santiaguito, activity consisted primarily of the extrusion of blocky lava flow from the summit of Caliente along with regular explosions emitting weak gas-and-ash plumes between 0.5 and 1.5 km above the vent (Global Volcanism Program, 2016b;Lamb et al., 2019;Wallace et al., 2020). During late 2015 and 2016 however, the activity was observed to shift in style and magnitude, as we noted a gradual decline in the number of explosions from its peak occurrence rate in early 2015, with explosions producing larger plumes Wallace et al., 2020). The explosions occurred with greater repose times, with large explosions occurring at a rate of less than once per day, and weak to moderate explosions occurring up to 4 times per day (Global Volcanism Program, 2017a). The large explosions generated larger proportions of pyroclasts as plumes rose up to 7 km above the vent (Fig. 1A). The largest explosions, which occurred in April-June 2016, excavated a large crater in the caliente dome structure (Fig. 1D;Global Volcanism Program, 2016a). During this period of heightened explosive activity, there was no extrusion of lava flows (Global Volcanism Program, 2016a). The explosions instead triggered pyroclastic density currents from column collapse and ejected large bombs up to 3 m in diameter to distances of 3 km (Fig. 1G;Global Volcanism Program, 2016a). During October 2016 the explosive regime came to an abrupt end, and the effusive regime resumed, accompanied with frequent smaller explosions at a rate of 25-35 per day (Global Volcanism Program, 2017a), similar to pre-2014 activity (Rose, 1987;Johnson et al., 2014). A lava dome emerged in the vent of Caliente in October 2016 (Fig. 1E;Global Volcanism Program, 2017a). During 2017 new lava continued to fill the excavated crater inside Caliente, leading to over-spill that caused block-and-ash flows (Global Volcanism Program, 2017b). Similarly to 2014/15, gas-and-ash plumes rose to heights of up to 800 m above the active vent, at a rate between 9 and 36 times per day (Global Volcanism Program, 2017b). By the end of 2017 the lava dome had appeared to be fully re-established its original (pre 2016) shape and INSIVUMEH reported that there was little change in activity between November 2017 and April 2018 with plumes rising approximately 500 m above the vent of Caliente (Global Volcanism Program, 2018a). During this time, between 15 and 21 explosions were commonly recorded per day. This level of activity decreased slightly until the end of the year, with explosions occurring between 11 and 15 times per day, producing plume heights of 500-800 m above the vent (Global Volcanism Program, 2018b).
Early work conducted on the data collected by the seismo-acoustic network can be found in De Angelis et al.  Lamb et al. (2019) analysed the seismic activity complemented by visual and thermal infrared observations to produce an early catalogue of 6101 explosions between 2014 and 2017. They presented a description of the types of volcanic processes and associated signals, and characterised the cumulative seismic energy which highlighted the increase in explosivity associated with the period of heightened activity recorded in 2016. A more recent study by Wallace et al. (2020) went further by constraining seismic energy, variable time delays between seismic and acoustic arrivals, thermal evolution and petrological changes associated with the increased explosivity. Characterisation of ash samples and ballistic bombs collected between 2014 and 2017 showed that changes in the chemical composition, mineralogy and groundmass texture throughout different eruptive phases occurred due to the fresh injection of a deepsourced, volatile-rich magma into the shallow magmatic mush leftover from protracted activity in the last decades, causing mingling and intensification of the explosive activity in 2016 (Wallace et al., 2020).
We have aimed to refine the identification of explosive events which led to the generation of a catalogue of explosions with improved completeness with respect to that used in Lamb et al. (2019). Algorithms used in the catalogue production constrained the occurrence of 18,895 explosions for the period between November 2014 and December 2018 (see Section 2.3). This dataset provides new insights into the volcanic and magmatic processes leading to shifts in eruptive style at Santiaguito. In this study we assess the seismic energy associated with individual explosions -a proxy for magnitude-and their variable occurrence rates; we exploit these measurements to investigate the magnitude-frequency relationship of explosions across different phases of eruptive behaviour with the aim of constraining four years of activity at Santiaguito.

Data acquisition
Between November 2014 and December 2018, we deployed and maintained a network of seismic and infrasound stations at the Santiaguito lava dome complex. The seismic network consisted of six Nanometrics Trillium T120 compact broadband seismometers (T = 120 s) and six Lennartz 3DLite short-period seismometers (T = 1 s). All instruments recorded with a sampling frequency of 100 Hz at a resolution of 24-bit. The twelve stations in the network (Fig. 2), were located between 810 and 7700 m from the active Caliente vent, and were deployed to achieve the best possible azimuthal coverage. Seven iTem prs100 infrasound sensors were also deployed, co-located with all broadband seismometers, as well as with the short-period seismometer at station LS01. Over the four years between 2014 and 2018 we conducted ten field campaigns in November 2014, April 2015, December 2015, January 2016, June 2016, February 2017, May 2017, January 2018, June 2018 and January 2019. Due to variable access to specific locations, and equipment malfunctioning, the stations in the network did not always operate simultaneously; however, with the exception of four data gaps, the network provided a quasi-continuous record of the activity over the four-year period (Fig. 2B).

Visual observations
Throughout the network deployment visual observations were routinely conducted by INSIVUMEH staff based at the local Santiaguito Volcano Observatory (OVSAN). These are documented in the Smithsonian Institution Global Volcanism Programme's weekly reports and monthly bulletins (published at: http://volcano.si.edu/volcano.cfm?vn=342030 last accessed 23/12/2019). Visual observations are frequently hindered by cloud coverage, and at night-time the local observatory is not manned. Furthermore, it is sometimes difficult to visually distinguish between passive degassing and small explosions. The observations made by INSIVUMEH were complemented by regular visits of the volcano of the Liverpool and KIT teams during the past 5 years. Despite the irregular nature of visual observations, they provide important information to link geophysical data to the ongoing volcanic activity. Visual observations were made throughout the recording period, making note of explosion frequency, plume heights, plume descriptions, lava dome appearance, and the presence of active lava flows.

Explosion catalogue
We designed and applied a detection and classification algorithm to extract and catalogue explosion signals from the continuous seismic data streams recorded across the network at Santiaguito between 2014 and 2018. The automatic algorithm consists of three main steps: 1) Pre-processing: Rapid signal detection at individual stations across the network based on signal to noise ratio (SNR) thresholding. To detect event waveforms in the continuous data streams, waveforms which had a SNR above 1.1 were extracted for the later processing steps. We fixed the noise level at a manually selected time void of signals. Although slower than a traditional STA/LTA, this ensured waveforms in periods of high noise, induced by non-explosive signals such as human activity, were also detected. To assist the detection of only explosions, the raw data streams were initially filtered between 0.1 and 10 Hz to remove noise in different frequency bands than explosions.
2) Processing: Producing a list of candidate explosion signals based on pre-selected waveform attributes. The detected waveforms were first compared to a preset group of frequency attributes. These attributes included the central and dominant frequencies, as well as the bandwidth at 50% of the dominant frequencies amplitude in the frequency spectrum. The envelope of waveforms which met the frequency criteria were then cross-correlated with an envelope of a template waveform, which consisted of a stack of manually detected waveforms. A cross-correlation coefficient threshold of 0.6 was set, with all waveforms correlating above the threshold labelled as candidate waveforms.
3) Post-processing: Using network association to determine if a candidate waveform can be catalogued as a true explosion. An explosion, originating within the network, should radiate seismic energy in all directions and be detected across the network by multiple stations. Noise however, could either be local to one station, or sweep across the network from one side to the other. The only noise which would be expected to behave in the same way as explosions is other volcanically generated events such as rockfalls, VT events or tremor. These events contain different frequency information, and therefore were removed during the processing steps. Candidate waveform detection times were compared across the network to determine if it was sufficiently detected to be classified as an explosion. Due to the changing station configuration, the criteria for sufficient network association to classify an event as an explosion changed accordingly. As the more reliable instrument, when multiple broadband stations were active, an explosion had to be detected at a minimum of 2 broadband stations. If only one broadband station was available, a detection had to be made by this station. In cases where no broadband station was active, an explosion had to be detected by 2 or more short-period stations. In the most extreme case when only one short-period station was active, the algorithm accepted all detections. However, to account for the change in criteria for a detection to be classified as an explosion, a 'trust' value was automatically assigned to catalogued events based on the number and type of stations which detected the event. If only one short period seismometer detected an event, a trust of 1 was assigned; if only one broadband or multiple short period seismometers detected an explosion, a trust value of 2 was assigned; a trust of 3 was assigned when one broadband seismometer and one or more short period seismometers detected an event; and a trust value of 4 was assigned when two or more broadband seismometers detected an explosion, with any number of short period seismometers. Following this post-processing step, an initial catalogue was produced, which was further refined through manual checks of events with a trust value of 1 and event clustering using cross-correlations. Events which did not cluster with other explosions, were manually checked.
The resulting catalogue contained 18,895 explosion signals, expanding the previous record  by a factor of approximately 3. The new catalogue includes a further 18 months of recording and separates explosions with smaller repose intervals, which were considered as one event in Lamb et al. (2019). For the range between November 2014 and May 2017, the two catalogues exhibit the same trends. Quality control on four weeks of our new catalogue of explosions was performed from January 2015, September 2015, February 2017 and September 2018 against the raw data streams and found the catalogue contained 90.5% of manually detectable events (610 out of 674), with fewer than 0.5% false detections in these periods. The four weeks contained high levels of noise, and with further checks on individual days throughout the recording period which maintained a comparable quality level, we believe the data checked in these weeks are representative of the most challenging conditions in the dataset. We checked for the influence of network configuration during a period of high-density station coverage in January 2015. We found that with a reduction of 50% of the available stations, and with the same criteria set for explosion detections, a drop of 15% in the detection rate is observed. Despite this drop we find that the trends in the data, and the results from the calculations made with this reduced data set vary minimally.
We calculated the size of all explosions in terms of their seismic radiated energy (SRE) and their associated energy magnitudes (M e ). The SRE and M e will be used to investigate the variability of the source dynamics for the different eruptive phases in the recording period, determine diagnostic features for these phases, and constrain relationships for the inter-explosion repose times.

Explosion energy calculation
The SRE was calculated for each event, defined as the elastic energy generated from an isotropic source at the surface of a homogeneous half space (Boatwright, 1980;Johnson and Aster, 2005). The equation is given as: where E s is the SRE, ρ earth is the rock density, C earth is the wave velocity, r is the radial distance from source to receiver, U is the ground velocity at the receiver, S is the site response and A is the attenuation. The integral is taken over a two-minute window which contains the full explosion waveform. Background noise included in the window is considered negligible to the calculation. Energy calculations vary from station to station due to unknown site effects, attenuation, and differences in the low frequency content between the short-period and broadband sensors. Therefore, relative station corrections were determined to obtain a consistent energy calculation of SRE from the network, so that all stations would obtain the same results irrespective of these differences. Station specific correction factors were calculated by determining the amplitude ratios for each station with LB03, which was used as the reference station. This is similar to the method of Lamb et al. (2019), who used LB01 as the reference station. LB03 was chosen as the new reference site to ensure consistency throughout the whole observation period. With variability between events, mean amplitude ratios were used as the relative station corrections for each station, taking into account re-deployments over the observation time span. Station correction factors used are shown in supplementary material, S1. To obtain robust energy values, SRE calculations were made using catalogued data from the station most active across the recording period to maintain consistency between calculations. The energy magnitude (M e ) for seismically observed volcanic explosions, based on ground motion velocity data, is a useful metric to assess the size of explosions at Santiaguito. Choy and Boatwright (1995) derived M e as: where Es is the seismically radiated energy, calculated from Eq. (1) (with constants later refined). We have calculated the energy magnitude for all events. As in all seismic catalogues, the smallest events may suffer from incomplete detection rates due to low SNR. The magnitude of completeness of a catalogue is defined as the minimum magnitude above which all events are reliably recorded. Using cumulative magnitude distributions, we calculated the magnitude of completeness for the catalogue to be 0.76 M e .

Evolution of explosive activity between 2014 and 2018
Eruptive activity, and in particular the nature and characteristics of explosions, evolved significantly at Santiaguito throughout the monitoring period; analysis of the new catalogue allows for tracking of the evolution of explosive activity during 2014-2018.

Explosion rates
The new catalogue provides a robust quantitative description of the rate at which explosions occurred during the four-year period. The catalogue shows highly variable activity, with explosions occurring over 500 times per week in January 2015, and b5 times a week in May 2016. Between these two dates, the rate of explosions decreased almost linearly. This decrease coincided with a shift from frequent small-tomoderate explosions to erratic and violent, ash-rich events (Global Volcanism Program, 2016a). Following return to a dominantly effusive regime in October 2016, our data indicate a rapid increase in the number of weekly explosions. By 2018 the rate of explosions became more consistent from week to week, with a steady level of detections recorded in the catalogue.

Secondary explosions
Over extended periods of time (~days to weeks), the explosions at Santiaguito often display an expected repose interval, statistically constrained at ca. 30-300 min; however, certain explosions are rapidly followed by a second, separate pulse of momentum-driven gasand-ash with identical acoustic signals within 10 min from the leading event (Fig. 3). We term these pairs as 'primary' and 'secondary' explosions. Secondary explosions commonly release less energy than the associated primary event, with 85% radiating b25% of the energy relative to the primary explosion (Fig. 3B). These paired explosions were commonly observed during 'normal' eruptive activity. Between 2014 and 2018, 4520 secondary explosions were found in the catalogue, making up 24% of all events. However, during the period of heightened explosivity in 2016, only 8% of the explosions were found to have repose times b10 min. We also observe no major difference between the seismic waveforms of primary explosions and isolated explosions that are not followed by a secondary explosion, as they display a high degree of similarity (Fig. 3A).

Energy trends
We used the weekly cumulative SRE to characterise different eruptive phases at Santiaguito (Fig. 4), calculated by summing the network SRE over consecutive 7-day periods. Although there may be bias introduced into the absolute values of SRE by the assumptions made in the site and path effects, our procedure ensures that relative changes are trustworthy. Changes in SRE, assisted by visual observations made during 2014-2018 indicate four distinct phases of eruptive behaviour. Phase 1 (November 2014-September 2015) is characterised by a high number of small-tomoderate explosions during a dominantly effusive eruption regime, with an increase in the cumulative weekly explosion energy observed. Towards the end of phase 1 explosions became less frequent but larger. During phase 2 (September 2015-October 2016), explosions are much less frequent and contain higher magnitude events which account for much of the weekly energy release. Phase 3 (October 2016-March 2017), which initiated when effusion of lava resumed at Caliente to fill the summit crater, is accompanied by an increase in SRE, caused by an high rates of occurrence of small-to-moderate explosions. Phase 4 (April 2017-December 2018) is characterised by continued effusion after a lava dome had become established within the crater. SRE remained low and slowly decreased due to a low but consistent number of small-to-moderate explosions. Although phase 4 contains a large gap in seismic data, visual observations do not indicate any additional change in behaviour during this time period.

Magnitude-frequency
Magnitude-frequency analysis is traditionally used to assess the ratio of small to large events within self-similar systems. Although most commonly used to analyse earthquakes, the GR relationship (Gutenberg and Richter, 1944) and Ishimoto-Iida's formula (Ishimoto and Iida, 1939), which link magnitude of earthquakes to their frequency of occurrence, have also been applied to explosive volcanic events (e.g. Deligne et al., 2010;Nishimura et al., 2016;Sheldrake and Caricchi, 2017;Rougier et al., 2018). Over the four-year recording period at Santiaguito, the energy magnitude of explosions varied over several orders of magnitude. We find a power-law dependency in our datasets between the magnitude of explosions and their rate of occurrence (Fig. 5), similar to the GR. In the logarithmic-logarithmic space (Fig. 5) the linear fit to a power law relationship follows: where f is the frequency of explosion occurrence, M e is the energy of an event, b is a parameter akin to the traditional b-value in the GR and a is a constant describing event productivity; the linear fit is performed only in the region above the magnitude of completeness of the catalogue. Over the entire four-year observation period the b-value is found to be 1.55 ± 0.06 (Fig. 5A). We investigated the variability of the magnitude-frequency relationship over the four phases of activity. We find that the b-value decreased from 1.84 ± 0.11 in phase 1 to 0.94 ± 0.06 in phase 2 (Fig. 5), reflecting the transition from small-tomoderate explosions in the effusive regime to large explosions in the explosive regime. In phase 3 the b-value increased to 2.28 ± 0.69 and remained high in phase 4 with a b-value of 2.40 ± 0.58 as the eruption regime produced fewer events during stable effusive behaviour. For all calculations only events over the magnitude of completeness are included, which we constrain as 0.73, 0.76, 0.90 and 0.77 M e for phases 1 to 4, respectively.

Repose times
Inter-explosion repose times indicate the occurrence and regularity of the explosions, and are commonly used to determine the statistical models which describe the probabilistic estimates which form the basis for eruption forecasting models (e.g. Connor et al., 2003;Watt et al., 2007;Lamb et al., 2014;Varley et al., 2018). We find that for the whole catalogue and for all phases of activity at Santiaguito repose times follow exponential distributions (Fig. 6). The linear slope of the repose time against the logarithm of occurrence for each repose time represents the rate parameter for the phase, that is the inverse of the mean repose time. The best fit to the repose times was made using an automatic piecewise regression function which tested the fit of two regression slopes at different breakpoints with a single regression line to find the best fitting slopes. The function compared the best fit with two regression lines with the best fitting single regression, if the regression of two slopes provided a significant improvement over a single slope, they were selected. We observe different rate parameters in each eruptive phase at Santiaguito, with phase 1 displaying two different rate parameters, indicated by a break in slope at a repose time of 0.17 days (i.e. 4 h). A break in slope is also observed over the whole catalogue, showing two dominant rate parameters over the four-year recording period (Fig. 6A). The rate parameters vary between 1.80 and 11.58, which occur in phases 4 and 1, respectively. The colorbar in Fig. 6a shows that the repose times in phases 1 and 3 have a time dependency, and do not exhibit a randomness of repose intervals, as observed in phases 2 and 4. The repose times in the overall catalogue also shows a time dependency.

Eruption phases and styles
Our observations agree with Harris et al. (2003) that the behaviour of Santiaguito is commonly characterised by effusion of blocky lava flows from Caliente lava dome, punctuated by small-to-moderate sized explosions reaching up to 1.5 km above the vent (Global Volcanism Program, 2016b), as seen in phases 1, 3 and 4. Previous studies at Santiaguito have stated that the explosions are regular in their occurrence with a repeatable source, generating explosions between 0.5 and 2 times per hour (e.g. Harris et al., 2003;Bluth and Rose, 2004;Sahetapy-Engel et al., 2004;Yamamoto et al., 2008). Over the four years of this study we see that the repose times are not so consistent, with event rates which varied between 500 times a week (3 times per hour) in January 2015 to b5 times a week (0.03 times per hour) in May 2016.
The combination of visual and seismic observations at Santiaguito revealed four eruptive phases; this builds on previous studies by Lamb et al. (2019) and Wallace et al. (2020), as we present for the first time a description of phase 4.
Eruptive phase 1 is characterised by the effusive eruption regime observed between November 2014 and September 2015, with high occurrence rates of small-to-moderate explosions. We attribute the progressive increase in energy observed in this phase, along with the decrease in explosion occurrence to the system preparing to transition into the explosive phase 2.
Phase 2 began in September 2015 with the emergence of larger magnitude explosions, which resulted in increased hazard to the surrounding population. These explosions were characterised by plumes which rose to heights of up to 7 km above the vent, and collapsed to generate pyroclastic density currents, while also ejecting pyroclasts of size up to 3 m diameter to distances of 3 km from the vent (Global Volcanism Program, 2016a). The large explosions excavated a deep crater in the Caliente lava dome. Wallace et al. (2020) split phase 2 into two phases, before and after March 2016, based on the componentry of collected volcanic ash which showed a progression from the dominance of dense brown clasts to porous and transparent glassy particles. For the purposes of this study, based on the lack of distinction in geophysical signals, these brief phases are simply combined in our analysis of explosion energy vs repose time (Figs. 4 and 6C).
Phase 3 saw the occurrence of effusive activity, interspersed with up to 200 small-to-moderate explosions per week, beginning in October 2016. SRE is observed to gradually increase during phase 3, which we associate with lava extrusion and dome growth throughout this phase (Fig. 4).
Phase 4 began after the lava dome had filled the excavated vent of Caliente in April 2017, and is characterised by attainment of a stable effusive regime. Throughout phase 4 the number and magnitude of explosions remained consistent at approximately 50 explosions per week with magnitudes up to 1.5 M e , which is reflected in the SRE trends. We note however, that statistically phase 3 and phase 4 are almost identical, as observed by their comparable b-values. We separate these phases through the observations of dome growth and upward trending weekly energy in phase 3, and continued effusion with an established dome and slowly decreasing weekly energy produced by a more consistent eruption rate in phase 4.
Phases 1-3 largely align with the phases outlined by Lamb et al. (2019), describing the same activity styles in each of the three phases. The differences between the phases here and those described by Lamb et al. (2019) are the boundaries between each phase. These differences are likely due to the methodology used to select the boundary times. In Lamb et al. (2019) the phase boundaries were chosen based on visual observations of the activity alone, whereas here, we use the SRE trends as a primary indicator, with visual observations used to corroborate the boundary choices.
From the observed trends in explosion rates and SRE, as well as the temporal distribution of repose times, we see that transitions between eruptive phases have occurred over different timescales, with gradual transitions occurring between phases 1 and 2 and phases 3 and 4, as well as a more abrupt transition between phases 2 and 3. The visual observations made throughout the four years again helped constrain these transitions and their timescales. The duration of the explosions, automatically calculated as the time between 2.5% and 90% of the cumulative seismic energy of the associated waveform, increased during the explosive behaviour in phase 2 to have a mean length of 25.8 s, compared to a mean of 18.9 s during the dominantly effusive regime of phases 1, 3 and 4. The duration of the explosions measured from the seismic waveforms matches closely with the visual observations of 30-60 s made by Bluth and Rose (2004) of the momentum-driven phase of gas release at Caliente. We speculate therefore, that it is likely that the seismic waveforms are a record of the vigorous gas venting phase while the fracture pathways remain open. The increased duration of the explosive events in phase 2 is possibly caused in-part by signal overprint from pyroclastic density current activity, and by larger gas overpressure developed over longer timescales in deeper parts of the magmatic column (cf. Wallace et al., 2020). Larger overpressures would enhance fragmentation efficiency (e.g. Kueppers et al., 2006), and deepen fragmentation, lengthening the travel distance in the conduit for gas and ash to erupt.
At the time of writing this paper, the activity at Santiaguito appears to be entering a new phase, with explosions becoming visibly larger, rising up to 1.3 km above the vent, occurring at increased rates of 35-40 explosions per day (Global Volcanism Program, 2019). This suggests that a transition into a fifth eruptive phase may have occurred since the end of our monitoring period in December 2018. Visual observations report that the size of the dome is currently larger than before the May   2014 dome collapse episode. From these visual observations, it is thought that the possibility of a dome collapse presents a significant risk, and such an episode could be hazardous should one occur at Santiaguito in the near future.

Source stability
We have noted through magnitude-frequency analysis that the explosions at Santiaguito obey a power-law relationship similar to the Gutenberg-Richter relationship used in earthquake seismology (Fig. 5). The power-law relationship depends on the b-value, which describes the magnitude of faulting events (i.e., ratio of small to large magnitude seismicity; Aki, 1967). In this paper we use the energy magnitude M e to describe the magnitude of explosion. We infer that the observed b-values reflect, in part, the state of the magma in the conduit system, which has controls on fragmentation. However, we argue that magma properties can vary both spatially and temporally, which makes it difficult to relate changes in bvalues to specific properties and thermo-kinetic conditions within the system (Roberts et al., 2015). Despite energy magnitudes spanning over 5 orders of magnitude, the power-law relationship indicates self-similarity between explosions at a mechanistic level (cf. Nishimura et al., 2016;Papale, 2018). Over the four-year period, Santiaguito expressed a b-value of 1.55 ± 0.06. However, across the four different eruptive phases we observed variations in the bvalue. The variations represent the evolution of spacing and size of events between phases; in other words, each phase may be considered as a discrete "experiment" with the sum of events defining its mechanistic character. The evolution from high b-value during phases 1, 3 & 4 (with frequent small events) to low b-value during phase 2 (with less frequent and larger events) suggests that the properties of the explosion source mechanism varied through time. Changes in source properties can involve differences in both the applied stress accumulation within the system as well as the scale and architecture of rupture, in part dictated by material properties.
The rupture and "strength" of magmas is controlled by several factors, including the fraction of heterogeneities (i.e., crystallinity and porosity), the viscosity of the silicate melt phase, and strain rate . As magmas are viscoelastic bodies (e.g., Dingwell and Webb, 1989), understanding magmatic fragmentation and resultant seismicity requires careful consideration of thermokinetic conditions (e.g., Papale, 1999;Zhang, 1999;Lavallée et al., 2012). The rupture of silicate melts results from an inability of the melt structure to relax an applied stress, provoking structural breakdown. The rate of structural relaxation of silicate melts is proportional to their viscosity, regulated by its chemistry and temperature (Dingwell and Webb, 1989); i.e. at low temperature, a melt's viscosity is higher and requires lower strain rate to rupture. So, understanding the rupture of silicate melts, or magmas (additionally hosting pores and crystals), requires knowledge of both viscosity and strain rates (e.g., Lavallée et al., 2008Lavallée et al., , 2013Cordonnier et al., 2009;Kendrick et al., 2013;Hornby et al., 2019). Material rupture results from the nucleation of micro-cracks, which propagate and coalesce in the build-up to system-size failure (Voight, 1989;Kilburn, 2003Kilburn, , 2012, and magma rupture ensues accordingly . Heterogeneities, commonly present in magmas, act as stress concentrators that focus the nucleation of micro-fractures (Sammis and Ashby, 1986); thus their presence lowers the strength of silicate melts (Vasseur et al., 2013;Cordonnier et al., 2012) and facilitates failure via characteristic acceleration in microseismic events that may be monitored to forecast rupture with increasing accuracy (Vasseur et al., 2015(Vasseur et al., , 2017. Zobin et al. (2014) speculated that a change in b-value associated with volcanic activity at different volcanoes may result from differences in magma crystallinity, which affected the viscosity of magmas and the stress required to failure. Whilst an increase in crystallinity may provoke changes in b-values of magma rupture, such a generalisation may be tenuous as the physico-chemical properties and stress conditions of magmas can vary widely between volcanic systems. As porosity generally has a greater impact on material strength than crystallinity (e.g., Coats et al., 2018), we anticipate that it would likely provide stronger controls on the development of material rupture and resultant seismic b-value, if other considerations (i.e., the viscosity of the melt phase and strain rate experienced) remained the same; this is supported by laboratory observations that single-phase melts rupture rapidly through localised fractures with large stress drops, whereas porous melts break slowly via multiple distributed small fractures with small associated stress drops (e.g., Vasseur et al., 2013). Likewise, a complementary study showed that the b-value resulting from magma rupture (under equivalent strain rates) generally increases with porosity (Vasseur et al., 2015). The development of porosity (Mueller et al., 2011) and pore pressure (e.g., Castro and Gardner, 2008) have previously been linked to changes in explosivity. Considering a single volcanic centre, as in our study, it is common to note a wide range of porosity in eruptive products (e.g., Lavallée et al., 2012;Mueller et al., 2011), yet only moderate changes in crystallinity within a given eruptive period (e.g., Bain et al., 2019), although magma viscosity may evolve regardless due to interstitial melt sensitivity to changes in temperature (e.g., Mastin, 2005;Blundy et al., 2006;Lavallée et al., 2015) and dissolved volatile content (e.g., Hess and Dingwell, 1996;Castro and Dingwell, 2009;Castro et al., 2005;Edmonds and Herd, 2007), and changes in ascent rate. Discharge rate and ascent rates can vary widely during volcanic eruptions; as a result, and crucial for the present seismic analysis, this implies that shear rates can vary by several orders of magnitude during explosions at lava dome eruptions (e.g., Quane and Russell, 2005). At Santiaguito, magma fragmentation has previously been linked to the rate of shear traction experienced in the magmatic column leading to explosions . Complementarily, laboratory testing has shown that shear rate plays a key role in the development of damage and associated b-value . In particular, laboratory experiments have shown that an increase in applied strain rate results in an increased localisation of fracture propagation during material rupture, which results in a decrease in b-value .
Thus, one needs to exert caution in inferring reasons for fluctuations in seismic b-values as many factors compete and in-situ magma properties and local thermo-kinetic conditions vary both spatially and temporally. Wallace et al. (2020) showed that changes in crystal textures, silica content and temperature (and thus viscosities) took place in the eruptive period studied here at Santiaguito. Hornby et al. (2019) demonstrated that magma porosity, temperature and applied strain rates were key controls on the tensile strength of Santiaguito magma, thus it is likely that some or all of these controlling factors have conspired to generate the b-values we resolved during the 4-year eruptive period studied at Santiaguito. Yet, it remains that small explosions which cause little-to-no damage of the lava dome result from small pressure releases, whilst large explosions on the other hand, release larger stresses in the system and tend to have ruptures with increased length-scales (Lavallee et al., 2013). If we return to our analogy that an eruptive phase represents an experiment, involving a material deformed within a bounding set of conditions (as identified by the seismic signatures and event spacing that define each phase) then the shifts in b-value between different eruptive phases demonstrate that explosions evolved from small, frequent events (high b-value) akin to the slow accumulation of damage during material deformation, to periods of larger scale, less frequent events which caused significant damage (low b-value), akin to wholesale failure during material deformation, driven by greater pressure accumulation. Thus, we must also stress that for meaningful interpretation, b-value must be considered over restricted and constrained time-periods, rather than considered as a stable, systemspecific value.

Secondary explosions
The data highlighted the occurrence of explosion duets. We define secondary explosions as explosions which occur within 10 min of an initial explosion, where the infrasound waveforms for both explosions are near identical yet separated by a period at a background level of activity. We observe that secondary explosions account for 24% of the explosion catalogue, making up a significant proportion of the explosive activity at Santiaguito. However, during the period of heightened explosivity in phase 2, the secondary explosions made up only 8% of explosions, predominantly following the smaller magnitude events during this phase, indicating that secondary explosions are a common feature of lava dome effusion at Caliente, occurring when the explosion fractures are more restricted. We also observe that over 85% of the secondary events radiate b25% of the energy compared to the primary event (Fig. 3B). We compared the infrasound signals of the primary and secondary explosions to investigate the relationship between the two events and show that the infrasound signals associated with primary and secondary explosions show a high degree of similarity ( Fig. 3C and D). T-tests on the time domain cross-correlations between all explosions infrasound waveforms show that the similarity between primary explosions and their secondary explosions is significantly higher than the crosscorrelation between any randomly chosen pair of explosions to a significance level of 0.1%. Example waveforms are shown in Fig. 4 C and D, where the waveforms of the 2 primary events show different shapes, yet the secondary events show high similarity to the primary events they follow. These observations indicate that the secondary events are somehow linked to the primary event, with a time dependence which cuts off at approximately 10 min. We speculate that the high correlation between primary and secondary explosions requires magmatic fragmentation under the same conditions; i.e. magma would likely fragment at the same depth, and the gas-and-ash products would be released via the same vent (geometry and size). We advance that this may be the case if the gas-and-ash are erupted from the same fracture pathways, as occasionally observed. Therefore, the repose interval may reflect the state of the fractures present in the lava dome. Due to the limitations in our observations however, we cannot validate these assertions.

Assessing future eruptive potential
The frequency of the inter-explosion repose times follows exponential distributions (Fig. 6), where the fit to the exponential distributions represents the rate parameter, which is the inverse of the mean return time. The rate parameters described by the exponential distributions are observed to change between eruptive phases as a result of the changes in source conditions and mechanisms between phases. Across the whole catalogue, two robust rate parameters are observed which relate to the mechanisms producing frequent low magnitude events during the effusive eruption regime, and infrequent large events in the explosive regime in phase 2 (Fig. 6A). The occurrence rate of different repose times is shown to correlate with time (colorbar - Fig. 6), which is a consequence of the explosion mechanisms and source conditions transitioning through time. During phases 1 and 3, there is also a transition from low to high repose duration through time, suggesting that the source behaviour is gradually transitioning (Fig. 6B,D). In contrast, phases 2 and 4 displays no correlation between repose duration and time, indicating that the source is more stable (Fig. 6C,E).
Above the cut-off magnitude (set at the magnitude of completeness), the stochastic explosion process can be represented by a Poisson distribution. Poissonian processes in explosion repose times are seen on both a global scale (De la Cruz-Reyna, 1991;Papale, 2018) and at specific volcanoes (e.g. at Volcán de Colima, Mexico; De la Cruz-Reyna, 1993). Marzocchi and Papale (2019) used the Poisson relationship for volcanic events of different sizes worldwide to determine the probability of events across varying magnitudes occurring within different time periods. However, on local scales at different volcanoes, inter-explosion repose times are commonly described by different statistical models such as Weibull and log-logistic, as well as Poissonian, which show variability in how explosions evolve within a system (Watt et al., 2007). Watt et al. (2007) showed that over time the statistical model describing the inter-explosion repose times can also change as the system evolves. On a global scale, only one exponential relationship is observed (Papale, 2018), however, we show that different phases are characterised by distinct frequency scaling of explosivity, which is likely caused by different source parameters and perhaps triggering mechanisms, which leads to the different rate parameters observed. Due to the nature of Poissonian processes, when an explosion has occurred, it is impossible to predict when the next event will occur, although probabilistic estimates can be given for the expected return time. Yet, probabilistic estimates of the expected return time are phase dependent (as each phase displayed different activity) and therefore change frequently on the timescale of several months. Other proxies may thus be necessary to enable the development of tools to ensure long-term assessments of eruptive behaviour at Santiaguito.
The maximum explosion magnitude expected within an eruptive phase can be estimated by extrapolating the linear fit of the magnitude-frequency distribution to a value of 1 event occurrence within the phase. For the entire catalogue, we obtain an estimate for the largest explosion to have an energy magnitude of 3.49, where the largest event recorded had an energy magnitude of 3.46. As with the overall catalogue, the estimates for each of the individual phases are overestimates. We calculated the largest events as 3.00, 3.92, 2.48, and 2.32 M e for phases 1 to 4, respectively, while the largest events observed had energy magnitudes of 2.82, 3.46, 2.11, 2.12 for phases 1 to 4, respectively. The level of caution in these estimates can easily be adapted; decreasing the event occurrence rate at which the estimation is taken increases the likelihood that the estimate of the largest possible event magnitude will be an overestimate. A caveat to this method is that the estimation assumes that the volcano will remain in the same eruptive regime, whereas the system has been shown to evolve rapidly, as observed between the dominantly explosive and effusive regimes in phase 2 and phase 3, respectively. Furthermore, the estimates made here are performed in hindsight; in real time, it may not be possible to determine which phase of activity is being exhibited, and for this reason we do not give errors in these calculations. Finally, changes in a volcano's eruptive behaviour affect the upper estimates of explosion magnitude, thus caution must be exerted when using it as a predictive tool.

Conclusions
Long-term seismic and acoustic monitoring at Santiaguito has revealed details on the changing nature of the explosivity at the active Caliente vent and revealed relationships between explosion energy and recurrences. Explosions occur at different intervals, ranging from a few minutes up to~6 days. On the shorter end of the scale, many explosions are followed within 10 min by secondary explosions (accounting for 24% of the total recorded explosions) with near identical acoustic signals to their primary explosion, and lower energy release. On the longer end of the scale, repose times between explosions however lead to contrasting signals and behaviour. Trends in the seismically radiated energy have provided an effective indicator of eruptive phase changes at Santiaguito, enabling the discrimination between one phase and the next. Magnitude-frequency analysis has shown that the b-value changes between eruptive phases. While the source mechanism of the explosions show a level of self-similarity, the changes in b-value suggest that the properties which control magma fracturing vary, thus local conditions cannot easily be reconciled. Yet, we infer that phases characterised by small, frequent events, resulting from small stress drops, have high b-values and more restricted damage, whilst phases characterised by large events, representing larger stress drops, resulted in lower b-values and more wholesale damage. Changes in the source properties between phases also influences the characteristic magnitudes and repose times, restricting extrapolation of behaviour to within a single phase and limiting the potential for long-term assessments of future trends in eruptive activity at Santiaguito.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.