Stochastic Modeling of Explosive Eruptive Events at Galeras Volcano, Colombia

Abstract

A statistical analysis of explosive eruptive events can give important clues on the behavior of a volcano for both the time- and size-domains, producing crucial information for hazards assessment. In this paper, we analyze in these domains an up-to-date catalog of eruptive events at Galeras volcano, collating data from the Colombian Geological Survey and from the Smithsonian Institution. The dataset appears to be complete, stationary and consisting of independent events since 1820, for events of magnitude ≥2.6. In the time-domain, Inter-Event Times are fitted by various renewal models to describe the observed repose times. On the basis of the Akaike Information Criterion, the preferred model is the Lognormal, with a characteristic time scale of ∼1.6 years. However, a tendency for the events to cluster in time into “eruptive cycles” is observed. Therefore, we perform a cluster analysis, to objectively identify clusters of events: we find three plausible partitions into 6, 8 and 11 clusters of events with magnitude the 6-cluster partition being the preferred. The Inter-Event Times between cluster onsets (inter-cluster) and between events belonging to the same cluster (intra-cluster) are also modeled by renewal models. For inter-cluster data, the preferred model is the Brownian Passage Time, describing a periodical occurrence (mean return time ∼36 years) perturbed by a Gaussian noise. For the intra-cluster explosions, the preferred model is the Lognormal, with a characteristic time scale of ∼0.9 years. In the size-domain, we analyze only single events, due to the low number of clusters. Considering two independent parts of the catalog, we cannot reject the null hypothesis of the erupted mass being described by a power law, implying no characteristic eruption size. Finally, looking for time- and size-predictability, we find a significant inverse linear relationship between the logarithm of the erupted mass during a cycle and the time to the subsequent one. These results suggest that, presently, Galeras is still in the eruption cycle started in 2007; a new eruptive cycle may be expected in a few decades, unless the present cluster resumes to activity with magnitude ≥2.6.

1 Introduction

All these models assume that the repose intervals in a given catalog are a sample from a common parent distribution that governs a stationary (that is, the parameters of the parent distribution do not change over time) process generating the data.

However, a frequent characteristic of volcanic activity is its tendency to concentrate during temporal frames, that is, to occur in temporal clusters. This has been studied in several papers (e.g., Conway et al., 1998; Cronin et al., 2001; Bebbington, 2010; Bevilacqua et al., 2016) in different volcanic settings. In these models there is no specific objective identification of possible groups of events (or clusters) that are sometimes observed in the real datasets. In other words, clustering in these models may appear as a consequence of the nature of the parent distribution; for example, a process described by a Weibull distribution may exhibit a degree of clustering (in case the shape parameters is smaller than 1), as well as self-exciting models.

Another development in the modeling of repose time intervals is to include also the size of eruptions, parameterized by the erupted mass, erupted volume, or the eruption magnitude (Pyle, 2000) or by the VEI (Newhall and Self, 1982). For example, several studies have searched for a significant link between sizes of eruptions and their subsequent repose times, that is known as Time-Predictable Model (TPM), or between repose times and the sizes of their subsequent eruptions, that is known as Size-Predictable Model (SPM), both in the classical (De la Cruz-Reyna, 1991; Burt et al., 1994) or in a generalized form (e.g., Sandri et al., 2005; Marzocchi and Zaccarelli, 2006; Bebbington, 2008; Passarelli et al., 2010; Garcia-Aristizabal et al., 2012; Sandri et al., 2017). These models have the advantage of providing a simple conceptual model linking the repose time intervals to the process of loading and discharging of the magma storage system.

Finally, in analogy with the famous Gutenberg-Richter law in seismology (Gutenberg and Richter, 1944), when a catalog of eruptions from a volcano is available, one may also wonder how to model the distribution of eruption magnitude occurrences, in a so called frequency-size law.

A new up-to-date catalog for Galeras volcano has been recently collated by integrating official catalogs from the Colombian Geological Survey and data from the Smithsonian Institution. The main goal of this paper is to provide a statistical description of the main features of the eruptive activity at Galeras volcano (Colombia), and to better understand the processes generating the observed data. This goal is achieved by applying the methodologies described above.

We analyze the Galeras' new catalog in different perspectives: the temporal domain, the magnitude domain, and the time-magnitude domain. In general, we mostly use standard statistical analyses to describe the volcano's eruptive behavior; moreover, we adapt these analyses in order to obtain clues about eruption recurrence considering clustered processes.

We present the results of our analyses in Section 5, following the same logical order as in Section 4. Finally, in Section 6 we discuss and compare the results obtained in the three perspectives with the aim of trying to understand better the magma feeding system of Galeras volcano.

2 Galeras Volcano

Galeras volcano is an andesitic stratovolcano 4200 m high, located at 1.221^°N and 77.359^°W in southern Colombia, close to the Colombia-Ecuador border (Calvache, 1990; Hurtado Artunduaga and Cortés Jiménez, 1997; Cepeda, 2020). According to Calvache et al. (1997) Galeras volcano is located in the scar left of the former Urcunina stage collapse, which left an amphitheater that opens toward the west, the active cone being the current Galeras volcano stage. The eruptive activity at Galeras volcano, which began about 5.000 years ago (Calvache, 1990; SGC, 2015), is mainly characterized by relatively small-volume, vulcanian eruptions (Hurtado Artunduaga and Cortés Jiménez, 1997).

Historical records of Galeras volcano activity only began after the Spanish colonization. Stories of the chroniclers resulted in an almost complete record of historical eruptions (Calvache, 1990) since the 16th century. Most of the eruptions occurred in the500 years of historical record can be considered to be explosive of small to moderate in size (Bonadonna and Costa, 2013), with some of them accompanied by small pyroclastic flows (Espinoza, 1988; Hurtado Artunduaga and Cortés Jiménez, 1997). However, the risk for the population can be high, as almost 500,000 people inhabits in the proximity of the active crater (DANE, 2011) and more than 5,000 people live within the high hazard zone (UNGRD, 2019).

3 Data and Definition of a Complete Data Set

3.1 Collation of the Explosive Eruptive Events Database

Galeras is one of the most active volcanoes in Colombia (Calvache, 1990). A considerable number of eruptions at Galeras has been recorded and identified in the last 5,000 years. Although, as mentioned above, the majority of these events has been generally characterized by vulcanian activity, at least six major events have been recognized and analyzed (Calvache et al., 1997). The third version of hazard map of Galeras volcano considers a total of 45 events found in stratigraphic records (Hurtado Artunduaga and Cortés Jiménez, 1997). An additional catalog of historical records of the last 500 years registered by numerous authors (Calvache, 1990) was compiled and documented by Espinoza (1988).

The catalog of Explosive Eruptive Events (EEE) used in this work was produced integrating geological and historical information obtained from different sources (Calvache, 1990; Hurtado Artunduaga and Cortés Jiménez, 1997; Stix et al., 1997), including the Smithsonian Institution catalog (Global Volcanism Program, 2013) and data from the Colombian Geological Survey (formerly INGEOMINAS), who provide detailed reports of the last periods of activity (INGEOMINAS, 2005; INGEOMINAS, 2010). A total of at least 78 eruptive events were registered between 3120 BCE and 2010 CE. Furthermore, the catalog includes calculated data as eruption date, VEI, magnitude, and, for recent eruptions, the volume and column height.

Smithsonian Institution and Colombian Geological Survey were the primary sources of the data. Both offer a catalog of the events and provide mainly information about the eruption date and the VEI value. Other characteristics as mass, volume, and column height were taken also from complementary sources (Calvache, 1990; Hurtado Artunduaga and Cortés Jiménez, 1997; Stix et al., 1997).

However, the information available for some eruptions (in general the older ones) is derived from historical descriptions that do not include technical details or the exact date; in these cases, where possible, the missing information was inferred using relationships among parameters available in the literature (i.e., Newhall and Self, 1982). The final EEE catalog produced in this work is provided in the Supplementary Material (Section Supplementary Material).

In order to proceed with catalog analysis, we fill in the missing months and/or days with random dates, without altering the recorded succession of events. If only the day is missing, we pick a random day in that month; if also the month is missing, we pick a random day in the year. To check the influence of this choice, we also filled in the missing dates by using the mid point of the missing period with no difference in the obtained results.

3.2 Defining a Complete Dataset in Time and Size

To study the temporal pattern of EEEs occurrence it is necessary that the input catalog represents a complete record of events above a given size, avoiding bias in the analysis by possible under-reporting of events. Therefore, we first analyze the completeness of our dataset by looking for significant changes in the rate of events in time.

By considering both the whole data (Figure 1) and the subset of events with magnitude M 2.6 since 1820 (Figure 2), we notice that the increase in event rate since the 1820s may indicate a bias due to the under-recording of small magnitude events before then. For this reason, we consider as complete the catalog since 1820s only for events above or equal to magnitude 2.6, and these time and size constraints are used to set the dataset to be analyzed in this paper. In doing so, we realize we are discarding the largest portion of the most recent recorded events (approximately since 1990), which probably represents the most accurate record of what has happened at Galeras, complete for much lower magnitude. However, on the one hand we cannot use only the latter 30-years data to infer the long-term behavior of the volcano, as we would end up with too few events to conduct a sound statistical analysis, representative for a short period of time; on the other hand, as reported in Section 4.2, we use this last 30-years portion of the catalog to build a frequency-size distribution that is then compared with the results obtained using the longer-term, complete catalog.

FIGURE 1

FIGURE 2

4 Modelling Approach

4.1 Temporal Analysis

One of the aims of this paper is to study the temporal distribution of EEEs at Galeras volcano. In Figure 2C we plot the cumulative number of EEEs at Galeras in time; a visual inspection of this plot may indicate their possible temporal clustering, in which explosions exhibit a tendency to form groups of events close in time. Figure 2D shows the cumulative erupted mass (in kg) as events occur in the complete part of the catalog. A visual inspection of this plot suggests also a possible relationships between the repose times between subsequent events and the erupted mass in one of the two events.

Keeping in mind these observations, in this paper we analyze the catalog of eruptions of Galeras volcano from different perspectives; in practice, we perform a dual analysis by studying in parallel recurrence features of single EEEs (we informally call this analysis “all IETs together”), and groups of EEEs occurring close in time, that hereinafter we call clusters. Beyond a classical recurrence analysis of eruptions, we are interested in examining the behavior of the clusters.

4.1.1 Recurrence Analysis of Single EEEs (“all IETs Together”)

As a first step we perform a classical recurrence analysis considering the complete part of the catalog. With this aim, EEEs are treated as a stochastic point process in which individual occurrences are assumed to be independent random events in time, at times . We define IETs, the time passed between two consecutive EEEs, as , . Adopting this approach, and ignoring the event duration, implies assuming that the EEE onset time is the most physically meaningful information (Garcia-Aristizabal et al., 2012); therefore, our modeling aims to describe the waiting times of the physical processes governing renewed volcanic activity.

In order to fit our complete dataset with renewal models, the data need to be first checked for stationarity (i.e., generated by a stationary process) and independence (i.e., successive IETs must be independent). Thus, we check independence between successive IETs of events in the complete part of the catalog by rejecting a correlation between the duration of successive IETs as in Bebbington (2012). Further, we check stationarity in the IET-generating process (in the complete part of the catalog) by running the Kolmogorov-Smirnoff test to reject significant deviations from a uniform distribution of the onset times of events.

In the light of these tests, the complete part of the catalog, consisting of independent IETs generated by a stationary process, is described by considering alternative renewal models. The performance of the different statistical models is evaluated through the Akaike Information Criteria (AIC), trying to find a good balance between goodness-of-fit to the data and model simplicity (expressed by the number of free parameters). In particular, if our IETs can be taken as a sequence of random variables independent and distributed according to a function , then the original series of events is called a renewal process. To better understand the nature of these competing models (and in particular the preferred one based on the AIC), we calculate and plot the hazard rate, which for a history-dependent point process is the ratio of the probability density function to the survival function S(x): is the instantaneous event rate at time t conditional on survival until the time t, with , being the occurrence time of the last event before t (e.g., Bain, 1978).

The renewal models that we test are the Exponential, the Weibull, the LogLogistic, the Lognormal, the Gamma and the Brownian Passage Time (BPT). Except the Exponential distribution, characterized by one parameter, these statistical models are described by two parameters (see Table 1). Further, we also test a mixture of Exponentials (3 parameters) and a mixture of Weibulls (5 parameters), in analogy with previous works (e.g., Turner et al., 2008), trying to explain a possible bimodality in the set of the IET taken all together. In order to fit the two latter models, we implemented a Markov chain Monte Carlo (MCMC) approach; however, we check the consistency of the results using also the algorithms provided by Bebbington (2012). In fitting these models vs. the cataloged data, we also use the censoring information, i.e., no event has occurred since October 4, 2007, up to the end of 2019.

4.1.2 Recurrence Analysis Assuming Clustered Eruptions

Observing the Galeras data set and the value of the Coefficient of Variation of the IETs, equal to 1.65 (e.g., Marzocchi and Zaccarelli, 2006), it seems plausible that EEEs exhibit some degree of clustering in the time domain. For this reason, we hypothesize the possibility that, once a new event or a new eruptive cycle starts, the short-term behavior of the eruptive activity may follow different patterns during the evolution of the activity.

4.1.2.1 Cluster Analysis

The k-means algorithm used for cluster analysis (MacQueen, 1967) takes the input parameter (k), which represents the number of clusters, and partitions a set of n objects into k clusters so that the resulting intra-cluster similarity is high compared to the inter-cluster similarity.

Defining an adequate number of clusters k is a challenging issue (see e.g., Garcia-Aristizabal et al., 2017); in this work we use the popular Silhouette Method (Rousseeuw, 1987; Xu and Wunsch, 2005; Al- Zoub and al Rawi, 2008; Kaufman and Rousseeuw, 2008). Given a partition into k clusters, this method assigns a silhouette coefficient to every object in the dataset; this coefficient depends on i) the average distance between the object and the objects belonging to its cluster, and ii) the average distance between the object and the objects belonging to other clusters. For a satisfactory partitioning, the latter must be larger than the former, and in this case the coefficient takes positive values.

The silhouette method offers the advantage that it only depends on the actual partition of the objects, and not on the clustering algorithm that was used to obtain it. Its usefulness in the interpretation and validation of cluster analysis results has been discussed in a number of applications (e.g., Rousseeuw, 1987; Chen et al., 2002; Bolshakova and Azuaje, 2003; Garcia-Aristizabal et al., 2017).

4.2 Frequency-Size Analysis

In the previous paragraphs we presented the methods for analyzing the distribution of Galeras eruptions in the time domain considering both single and clustered eruptions. However, determining a frequency-magnitude (F-M) relationship for natural phenomena is another key issue for hazard and risk assessments since it allows one to estimate the return period of damaging events. As for other natural hazards, the frequency of explosive eruptions decreases as their size or magnitude increase; while for some phenomena (as e.g., earthquakes, extreme meteorological events, or flooding) there exist event catalogs populated with enough data for quantitative analyses, the record of eruptions from a single volcano is usually composed by a few events for which it is often difficult to properly define a F-M model. For this reason, the most sound F-M models proposed to describe the distribution of eruption magnitudes rely on regional to global analyses in which large eruptions from different volcanoes are merged and generalized in a single model (e.g., Simkin and Siebert, 2002 onwards). However, it is worth to note that while such aggregations may show an identifiable F-M distribution, it is not implied that this relationship holds for any particular volcano.

In this context, different F-M models for describing volcanic explosions have been proposed in literature, which are based, for example, assuming a exponential magnitude distribution (similarly to the Gutenberg-Richter model in seismology, as e.g., De la Cruz-Reyna, 1991; Mendoza-Rosas and la Cruz-Reyna, 2008), which is equivalent to assuming a power law distribution of eruption volumes (e.g., Papale, 2018) or mass, and it implies a so-called scale invariance in the eruption magnitude. Others authors have used extreme value theory (e.g., Mason et al., 2004; Coles and Sparks, 2006; Mendoza-Rosas and la Cruz-Reyna, 2008; Deligne et al., 2010; Furlan, 2010; Sobradelo et al., 2011).

In this study we try to fit our data with a power law describing the erupted mass in single events:where is the number of eruptions with erupted mass equal to m, and k and b are parameters of the law. By taking the logarithm of this equation, we get to the formulation that is similar to the Gutenberg-Richter law for earthquakes, i.e.:where M is the magnitude of an eruption, which is a linear function of (Pyle, 2000). If this holds, magnitudes obey to an exponential distribution (Aki, 1965; Marzocchi and Sandri, 2003):where , being the completeness magnitude in the catalog, is the magnitude bin (in our case 0.1 to 1), and ε is uniform random noise over mimicking the magnitude uncertainty.

Since we have already estimated the completeness magnitude in Section 2, we apply the scheme above for the data from 1820 and . However, since we also recognize a lower completeness magnitude if we reduce the catalog to the last 30 years (in particular, from 1990 with , see bottom panel in Figure 3) we repeat the above analysis also for this different dataset, which shares only one eruption (October 4, 2007) with the one from 1820 and . If we cannot reject the hypothesis of magnitudes following an exponential distribution for both datasets, we may compare the two b-values obtained and see whether they are consistent: in a scale-invariant process, they should have a unique value.

FIGURE 3

4.3 Time-Predictable and Size-Predictable Models

The recurrence analyses presented so far provide important clues regarding the stochastic behavior in time of eruption events above the minimum completeness magnitude of eruptions. Complementary to such analyses, it is possible to explore the predictive value of other data in the catalog, in particular the possible relationships between the size and the time interval between eruptive activity. Such kind of models are generally referred to as time-predictable and size-predictable (e.g., De la Cruz-Reyna, 1991; Burt et al., 1994; Passarelli et al., 2010; Sandri et al., 2017).

We test such models for both the individual eruptions and the clustered eruptions. Regarding in particular the clustered data, on the complete part of the catalog partitioned in clusters, we calculate the total erupted mass within each cluster (taken as an eruptive cycle); afterward, we analyze this information considering the time preceding or the time following the respective eruptive cycle (taken as the time between the end of a cluster and the onset of the following one), searching possible significant relationships typical of time- or size-predictable systems.

Here, we generalize the classical conceptual models at the basis of time- or size-predictable systems as in Passarelli et al. (2010), by treating clusters as eruptive cycles, and searching for relationships between the “size” of a magmatic cycle (parameterized by the total erupted mass in the cycle) and the time elapsed from/to the previous/following cycle. We explore both classical linear relationships (as in De la Cruz-Reyna, 1991; Burt et al., 1994), as well as more general ones, such as a power law (as in Passarelli et al., 2010). In this view, any significant relationship between the total mass erupted in the preceding cycle and the time to subsequent eruptive cycle can be taken as a generalization of time-predictability. Conversely, any significant relationship between the time from the preceding eruptive cycle and the total mass erupted in the subsequent cycle can be taken as a generalization of size-predictability.

5 Results

5.1 Temporal Analysis

5.11 Recurrence of Single Eruptive Events

We analyze the time series of the IETs from 1820s above magnitude 2.6, to exclude any significant dependency between successive IETs of events with M ≥ 2.6 (Figure 4). We reject any significant correlation between the duration of successive IETs, both in terms of linear correlation (Pearson coefficient = 0.08, p-value ) and rank correlation (Spearman coefficient ρ = 0.34, p-value ).

FIGURE 4

The stationarity Kolmogorov-Smirnoff test on events with M ≥ 2.6 since 1820 (Figure 2C) is not rejected at 1% significance level (value of test statistics is 0.25), and we see that our data are well within the 99% confidence level. Thus, we do not reject the hypothesis that our data represent a stationary time series.

In the light of these tests, we take the list of events with M ≥ 2.6 since 1820 as a complete dataset, composed by independent IETs generated by a stationary process, and we fit it with the most common and simple renewal models (Figure 5), listed in Table 1 (where we show the models tested, their maximum likelihood parameters, and their AIC value). According to the AIC, the best model that describes the IETs of all the eruptions in the complete catalog is the Lognormal distribution (in bold in Table 1).

FIGURE 5

The mixture of exponentials with the best AIC shows an intermediate degree of mixing between two different exponentials, and explains the data certainly better than the single exponential (see Table 1); however, its AIC is not better than the one achieved by the Lognormal distribution.

The mixture of Weibulls fit by our MCMC algorithm systematically points toward solutions with the p parameter close to one; such solutions imply the prevalence of one of the Weibull components, which is characterized by the same parameter values of the single Weibull reported in Table 1. On the other hand, the fitting procedure by Bebbington (2012) points to an alternative solution with a good AIC but characterized by one of the two Weibull distributions being close to a Dirac delta, i.e., very peaked and not useful to understand any process generating these data. The difficulty to fit this model probably arises because of a complex likelihood surface (characterized by 5 parameters) to be maximized, making it hard to find a proper solution to the problem with the few tens of data available; ultimately, this may be related to a difficulty in discriminating between two processes hardly distinguishable with the data at hand. Thus we conclude that the mixture of Weibulls is not helpful in our case, and that our best model when considering all the IET together is the Lognormal.

5.1.2 Recurrence Analysis Assuming Clustered Eruptions

5.1.2.1 Optimal Cluster Partitioning

The cluster analysis is performed in the time domain and considering the complete part of the catalog (M ≥ 2.6 data since 1820). We test different k values (where k is the number of possible clusters) and, considering the average silhouette value (Figure 6G), we identify at least three plausible clustering partitions: (Figures 6A,B), (Figures 6C,D) and (Figures 6E,F). In each of the three rows in Figure 6, the proposed k value was used as input parameter in the k-means algorithm for the identification of clusters; the resulting clusters of events are shown in the right panels (Figures 6B,D,E), in which onset times are plotted against the respective event magnitude.

FIGURE 6

These three possible cluster divisions exhibit high average silhouette values (vertical colored lines in Figures 6A,C,E) indicating relatively good partitions. In the case of 8 clusters however, one of the silhouette coefficients is negative (a negative value means that the event is more similar to the average event in another cluster than in its cluster). In the other two partitions (6 and 11), most of the clusters exhibit silhouette values above the respective average. However, in both cases there are poorly constrained clusters: 1 cluster in the 6-cluster partition and 2 clusters in the 11-cluster partition (Figures 6A–E) show a silhouette value not reaching the average value; furthermore, in the case of 11 clusters, 9 clusters consists of only one or two events having M 2.6.

Taking this latter limitation into account, we choose the simpler 6-cluster partition as the preferred one for analyzing the catalog of Galeras volcano events. Since the choice is ultimately subjective, in Section 6 we will also briefly discuss and compare the results for the alternative eight- and 11-cluster partitions.

Once the clusters are identified, the clusters start times are set using the onset time of the first event within each cluster. The start time of the 6 clusters identified in the Galeras data set (from 1820 and ) are summarized in Table 2, along with the number of events within each cluster. In the Supplementary Material we provide the same table for the eight- and 11-cluster partitions.

Stochastic Modeling of Inter-cluster and Intra-cluster Eruptive Behavior

Figure 7 shows the Cumulative distribution (CDF) of both the inter- and intra-cluster IETs, and the different models tested on the two datasets, in the case of 6-cluster partitioning, and considering the censoring data.

FIGURE 7

Table 3 shows a summary of the AIC values obtained for the different alternative models tested for describing the inter-cluster IETs (i.e., long-term process). In bold we highlight the preferred model, with lower AIC, that is the BPT.

Table 4 shows a summary of the AIC values obtained for the different competing models tested for describing the intra-cluster IETs (i.e., short-term process). In bold we highlight the preferred model, with lower AIC, that is the Lognormal, as in the “all IETs together” case. This similarity (in terms of model type, and consistency in the best-fitting parameter values) can be explained by the fact that the two datasets used for these two analyses are composed by the same data except the IETs between the pairs of events marking the end of a cluster and the beginning of the successive one. Given the low number of clusters (6), it means there are only five IETs less in the intra-cluster dataset with respect to the “all IETs together” one.

5.2 Frequency-Size Analysis

For each of the two subsets of possible complete datasets considered in this work (that is, subset 1 (S1): eruptions since 1820 with ; subset 2 (S2): eruptions since 1990 with ; see Figure 3) we apply the magnitude randomization procedure using the transformation shown in Eq. 4; in practice, we subtract to each magnitude value, and sum an ε randomly sampled from an uniform distribution in the interval [-0.05; 0.05] (assuming ). The Lilliefors test does not reject the null hypothesis of magnitudes being exponentially distributed for the S2 dataset (p-value ), where magnitudes are more precisely assessed. For the S1 we need to apply a larger noise to the estimated magnitudes in order to not reject the null hypothesis (with or larger, up to 1 unit of magnitude, p-value is always ). This is very reasonable, considering that these are older eruptions and the estimated magnitude is affected by larger uncertainty. By applying the maximum likelihood method we obtain consistent estimates of the b-value for the two datasets (Figure 3): and .

5.3 Time-Predictable and Size-Predictable Models

Figure 8 shows the plots of the total erupted mass against the time to the subsequent EEE (Figure 8A) or the time since the preceding EEE (8b) considering the data of each single event independently. In these two cases, the data are not fit by any significant TPM nor SPM model.

FIGURE 8

We perform the same analysis for the clustered EEEs; in this case, the erupted mass in a given cluster is plotted against the time from the end of that cluster to the start of the subsequent one (8c) or the time between the end of the previous cluster and the onset of the one erupting that mass(8 days). According to the plot in Figure 8C, there is some evidence (p-value of F-test: 0.0096) for a possible “inverse” time-predictable relationship in which eruption clusters characterized by large erupted mass are followed by shorter repose times (until the start of the following eruptive cycle). The fitted model shown in Figure 8C implies a power law relationship between time to the subsequent cluster and erupted mass:where (f-test with : is rejected with p-value ), and . According to this model, in case the current cluster does not resume to activity of magnitude 2.6 or above, given the present mass output from the 2007 eruption, we should expect about 30 years (from 2007) for the next cluster to start.

6 Discussion

A new catalog of eruption explosive events at Galeras was collated by using the datasets from the Colombian Geological Survey and data from the Smithsonian Institution.

We consider that the catalog is complete, stationary and composed by independent events since 1820, for events with . This catalog, taking it as a whole dataset of the realization of a homogeneous point process, is best fitted by a Lognormal model, implying that the logarithm of the IETs, taken all together, are normally distributed. The μ parameter of the Lognormal distribution represents the mean of the natural logarithm of the cataloged IETs, so that is their geometric mean; the latter can be taken as a characteristic timescale, which in this case has a value of approximately 1.6 years. However, being a geometrical mean, its value is strongly influenced by the shortest IETs.

Along with the CV value, some of the fitted distributions to the IETs all together have parameters pointing at a clustering of the EEE. For example, the parameter in the BPT model (reported in Table 1) is an indicator of the aperiodicity of the process (Garcia-Aristizabal et al., 2012) and here has a value much larger than 1, as does the ß-value of the Weibull model. These indicators support the visual observation of clustering in the EEE occurrence. Thus, we are able to identify different plausible cluster partitions. From these, the partition into 6 clusters seems to be the one providing the best trade off in the distance measures used for grouping the events and therefore we take this one as the reference partition for the discussion. Nevertheless, in both the main text and the Supplementary Material we complement the description presenting the respective analogue results obtained when using a different partition.

Considering the recurrence analysis of the eruptive events within an eruptive cycle, the preferred stochastic model describing the intra-cluster time intervals for this data set is the Lognormal (integrating the data from all the clusters), with a characteristic timescale parameterized by the μ parameter, translating into a geometrical mean of intra-cluster IETs of approximately 0.9 years. From Table 4 we see that a comparable fit is achieved by the LogLogistic distribution, with the same characteristic time scale (again representing the geometrical mean).

Conversely, considering the recurrence analysis of the eruptive cycles (i.e., long-term process), the preferred stochastic model describing the inter-cluster time intervals is the BPT, with a characteristic repose time (scale parameter) of approximately 36 years. The parameter α shown in Table 3 in the present case takes a value around 0.4, which indicates a slight tendency toward periodical behavior (see e.g., Garcia-Aristizabal et al., 2012). Although the BPT is the preferred model, all the other distributions except the Exponential fit the data similarly. Such similar values of the LogLikelihood are due also to the low number of clusters available in the record. In any case, we keep the BPT as the model best-explaining the inter-cluster onset times.

As mentioned in Section 4, in a history-dependent point process, the conditional probability that an event happens in a time interval , given an interval of years since the occurrence of the previous event (with t the current time of the assessment, and the time of the last event, so that x is the time from the last event), is the hazard rate or hazard function. The hazard function describes instantaneous failure rate, or in other words the conditional probability density of failure at time t given that no event occurred in the time interval x. An increasing hazard function at time t indicates that an event is more likely to occur in a given increment of time than it would be in the same increment of time in an earlier age. The hazard function describing the probability to experience the next event in the time interval given no event has occurred in x (that is again the time elapsed since the last event) can be calculated as follows:where T is the occurrence time of the next event, is the probability density function and is the cumulative distribution.

The implications of the different results in terms of renewal processes, found so far, can be analyzed inspecting the hazard rate function (see Figure 9A; in the Supplementary Material we provide the same plot for the 8- and 11-cluster partitions). The hazard rate of the Lognormal distribution fitted to all the IETs (without considering clustering of events) is a sharp increasing and then slowly decreasing function. It broadly means that the conditional probability of a new event is high shortly after an event, and slowly decreases as time progresses. This is a clear indication of a clustering effect. Given the current information that no event has occurred during the current censoring time x (approximately 12 years, since 4th October 2007), the probability of such an event in the next year can be calculated by Eq. 6 (Garcia-Aristizabal et al., 2012), with = 1 year. By using the best fitting model parameters from Table 1, we have a probability of about 6% in the year 2020. In Figure 10 (Box A) we also give an idea on how well this estimate is constrained: by sampling the Lognormal parameters within their uncertainties (from a bivariate Normal distribution and accounting for their covariance), we recalculate the probability of a new EEE in the year 2020. The 95% confidence interval for the probability of a new EEE in 2020 under this model ranges from 4 to 10% (Panel A1).

FIGURE 9

FIGURE 10

We can use the analysis of cluster to assess how likely a new event would belong to the ongoing cluster, or would mark the onset of a new cluster. The analysis on the 6-cluster partitioned data shows that the hazard rate function best describing the intra-cluster IETs distribution is again a Lognormal distribution, similar to the “all IETs together” one: the longer the time from the last events within a cluster, the lower the likelihood of having a new event. This well describes the gradual decline of activity during an eruptive cycle. Given the current repose time of about 12 years since the last event in the last cluster, the probability of a new M 2.6 in the year 2020 is about 9%, quite similar to the one obtained considering “all IETs together”. This is reasonable, considering the similarity between the two datasets (see Section 5.1.2). In Figure 10 (Box B) we see that the 95% confidence interval for the probability of a new EEE in 2020 belonging to the current cluster spans from 6 to 15% (Panel B1).

The preferred model describing the inter-cluster time intervals is a BPT, which can be seen as a delayed Poisson process: in our case, the hazard rate is very small for about 15 years after an eruptive cycle has begun (Figure 9A), and increases monotonically up to around 65 years, when it becomes flat as in the case of a homogeneous Poisson process, characterized by a constant hazard rate (e.g., Garcia-Aristizabal et al., 2012). The probability of a new eruptive cycle beginning becomes non negligible after a time passed since the start of the previous cycle of about 16 (or 20) years, i.e., when the survivor function is smaller than 99% (or 95%). Given the current elapsed time since the beginning of last cluster of about 12 years, there is a probability of about 0.1% of entering in a new eruptive cycle in the year 2020. In Figure 10 (Box C) we see that the uncertainty on the BPT parameters propagates into a 95% confidence interval on such value that spans from ∼0 to 4% (Panel C1).

A question we may try to answer is when to consider a cycle over, once Galeras is in an eruptive cycle, according to the clustered model for events. By using the fitted model parameters for the 6-cluster partition in the intra-cluster Lognormal model (Table 4), a repose of about 47 and 15 years yields respectively a 1 and 5% probability in the survivor function.

These results suggest that, overall, at present we are more likely still in the cluster that started in 2007: in other words, if the volcano re-enters in activity with or above during the current year (2020), according to this clustering model this EEE would be more likely associated to the ongoing eruptive cycle, rather than to the start of a new eruptive cycle. Speculatively, a slightly different judgment would be given in case no EEE occurs in the next 20 years approximately, i.e., until the year 2040: in such case, the probability in that year to experience an EEE overall from the “all IETs together” model would be about 3%, and similarly of an EEE belonging to the present cluster (4%), while the probability of an EEE opening a new cluster becomes about 6%. This is just an example to show how the results from the cluster model can be interpreted, keeping separate the intra- and inter-cluster models. In the right-column panels of Figure 10 we show the resulting probability for an event in year 2040, under different models, accounting for the uncertainty on the model parameters. The respective 95% confidence intervals for the three probabilities are (2–5), (3–8) and (3–18)% (Panels A2, B2 and C2, respectively).

Overall, when considering the “all IETs together” (as well as the intra-cluster) modeling, we obtain a best-fitting model that describes a decreasing hazard rate (Figure 9A) in time. This implies that, after a long time has elapsed since the last event of , the eruptive activity tends to vanish. In this view, the cluster model, through the description of inter-cluster time intervals, is useful to describe the recognized long-term behavior of Galeras, characterized by time-separated eruptive cycles at least in the last 200 years. As shown in Figure 9B, where we plot the ratio of the inter-to the intra-cluster hazard rate, the hazard rates for these two processes are comparable for a censoring time of approximately 30 years, a sort of threshold for discriminating which of the two processes is more likely to be the origin of the current censoring time.

The approach adopted here to purely describe the intra- and inter-cluster occurrences in the catalog has a remarkable limitation: both the distributions adopted for the two processes share the same time domain (that is the entire real axis both for IETs within and among clusters). This limitation does not allow us proposing a single merged model for eruption forecasting based on the clustering approach. In principle, it would be possible to overcome this limitation and create a single consistent model by truncating the IET distributions over separate time domains for the inter- and intra-cluster intervals. However, a major problem of such a model in the present case, characterized by few tens of data, would be to set the time boundaries necessary to truncate the two distributions in a meaningful and robust procedure. This is the reason why we prefer not to truncate them, and stick to a description of the catalog occurrences within the two identified regimes, separately.

The frequency-size analysis, although performed on a dataset restricted to the last few decades of activity, returns b-values that are consistent when considering two separate subsets of the catalog, spanning different periods of time and magnitude range. This indicates a common power law on the erupted mass, linking eruptions of different scale and energy in a scale-invariant relationship, and implying that there is no characteristic size in the eruptive process. We acknowledge that this is not a sufficient proof, but only a possible hypothesis not confuted by the data we have available at this stage.

In principle, we could apply the frequency-size analysis to the clustered data, to find a distribution for the cluster magnitude (i.e., the total erupted mass in an eruptive cycle), as in the time- and size-predictable analysis. However, we objectively identified only 6 clusters, and this number is too low to compute a b-value. However, this would represent an innovative analysis in eruption catalogs, and we recommend it as a possible future work for detailed volcanic records in which many clusters can be isolated.

Finally, we have found that the logarithm of the total erupted mass in an eruption cycle shows a significant inverse linear relationship with the time to the following eruption cycle. In other words, the higher the total erupted mass during an eruptive cycle, the shorter the time until the start of the successive eruption cycle. According to this model, if the current cluster will not experience a mass output much larger than the present output (from the 2007 eruption) we should expect about 18 years from now (or about 30 years from last EEE with magnitude ≥2.6 in 2007) for the next cluster to start. It is interesting to note that this time is similar to the one required by the inter-cluster model for a new cluster to begin (survivor larger than 5%), and larger than the probability of an intra-cluster EEE.

Many of the results found in this study on Galeras may be explained by a simple conceptual model of a volcanic system fed by two interconnected shallow and deep magma reservoirs, as already proposed for e.g., Soufriere Hills, Montserrat (Melnik and Costa, 2014), or Colima, Mexico (Massaro et al., 2019). The deep reservoir feeds the shallow one through a dyke-like structure, and the magma transfer from the deep to the shallow reservoir is strongly controlled by the pressure gradient between the two reservoirs. An eruptive cycle (originating events that determine the intra-cluster IETs) starts when the magma in the shallow reservoir overcomes a threshold in a state variable, for example the overpressure, which can be the result of renewed magma fed from the deeper reservoir. Meanwhile, the long-term recharging of the system (originating the inter-cluster time intervals), which can be primarily modulated by the pressure gradient between the two reservoirs, may take place in the background; magma flux from the deeper chamber is likely stimulated as the pressure gradient increases while an eruptive cycle goes on, revitalizing the system to start a new eruptive cycle. When looking at the all IETs together, we best explain the IETs by a Lognormal distribution characterized by a mean linked to the geometrical mean of the IETs, and by a thick tail that accounts for longer reposes, characteristic of the times between the last event in a cycle and the first event of the subsequent one. The erupted mass in each event does not show a preferred size. However, when considering clustering, we best explain the inter-cluster onset times by a BPT model. Such model is usually interpreted as a realization of a point process in which new activity will occur when a threshold is reached, the overpressure of the shallow reservoir in our case. In the time between two eruptive cycles, the shallow system is recharged by the deep one: this reloading is governed by two components, that are a constant-rate loading component (also called the drift term, for example a constant-rate recharging from the deep to the shallow reservoirs) and a random component (also termed Brownian motion) that may be due to random stress perturbations on the reservoir system (for example in terms of variations in the stress field due to regional seismicity or activity at other nearby volcanoes, or variations in the deep recharging process). The latter component may thus randomly disturb the volcanic system causing an aperiodicity in the mean repose time between eruptive cycles, although in the case of Galeras we have found that this aperiodicity is weak. The inverse TPM found for Galeras tells us that the eruptive cycles characterized by a large total erupted mass can create a higher pressure gradient between the deep and the shallow reservoirs, favoring the flux of magma feeding the shallow reservoir. Conversely, eruptive cycles with a low total erupted mass do not significantly perturb the pressure in the shallow chamber, hence the pressure gradient between the two reservoirs, inhibiting the “eruptability” of shallow magma.

7 Conclusion

We have analyzed a catalog of explosion events at Galeras volcano, with different perspectives, to extract the significant features of event occurrence when considering onset times, eruption magnitudes, and the two together.

We first identify as complete the part of the catalog from 1820 with eruption magnitudes 2.6 or larger.

On this complete catalog, the IETs are best described by a Lognormal model, with a characteristic IET of about 1.6 years. However, there is a clear tendency for explosion events to cluster in time, originating eruptive cycles. The identification of such cycles is inherently subjective: by adopting a set of reasonable criteria, we identify a preferred partitioning of the events from 1820 with magnitude 2.6 into 6 clusters. Under this partitioning, the distribution of events within a cycle is best explained by a Lognormal distribution, while onset times between clusters are best explained by a Brownian Passage Time model, that can be seen as a delayed Poisson process. The added value of the clustered model lies in its ability to describe the re-activation of volcanic activity, through the modeling of inter-cluster time intervals, with respect to the best-fitting model on all the IETs together: this latter model forecasts a gradual decline in the probability of a new event, after an initial peak, which would unrealistically lead to a vanishing of the explosive activity in case no event occurs for a prolonged period.

The eruption magnitudes of single eruptive events seem to be governed by a power law distribution, implying that there is no characteristic size for single events.

However, when considering the time between eruptive cycles and the volume erupted during each cycle, we observe an inverse time-predictability: the larger the volume erupted in the previous cycle, the shorter we will have to wait for a new cycle to begin.

These results suggest a simple conceptual model, consisting of a plumbing system characterized by two magma reservoirs located at different depths and connected by a dike system. In this model, the deep storage feeds magma into the shallow one through a constant-rate process that is randomly perturbed (for example by regional or local stress variations), triggering a new cycle onset when a threshold in the overpressure of the shallow reservoir is reached: a cluster of events in time occurs, consisting of many small events or few large ones (following, in general, a power-law distribution), relaxing in turn the shallow reservoir. If the total erupted volume in the cycle is large, it favors the uprizing of new magma from the deep reservoir to the shallow one by an increased pressure gradient between the two. In this way, the waiting time for a new cycle to begin is shorter.

According to the monitoring observations of the Servicio Geologico Colombiano (SGC), along with no significant eruptive activity at Galeras in about a decade, there has been also lack of geophysical signals that unambiguously indicate new magma feeding the system. Under this perspective, the last cluster seems to be over, and a new eruptive event could be classified as the first one of a new cluster. However, considering our probability estimates of having either a new eruption in the current cluster or the start of a new cluster next year, the preferred 6-clusters partition indicates a larger probability for a potential new event in the next year as belonging to the previous cluster rather than opening a new one, whereas when considering a 11-cluster partitioning, these two probabilities are very similar. It is worth noting that our analysis strictly focuses on the pure temporal occurrence of events (above a relatively high completeness magnitude), regardless of all the other observables that may further constrain the onset, or the end, of a cluster, such as the presence or absence of seismic activity, deformations, or geochemical anomalies. Including such observations to complement any clustering approach would be of invaluable support for analyzing recent and frequent activity, whereas it is less supportive for a long-term analyses, as the one presented in this article.

Funding

This research has been partially funded by project EXACT (Progetto 29 di Ricerca Libera INGV 2019 Delibera 53/2020).

References

• 1

  AkiK. (1965). Maximum likelihood estimate of b in the formula log (n) = a bm and its confidence limits. Bull. Earthquake Res. Inst. Univ. Tokyo43, 237–239.

  • Google Scholar

• 2

  Al-ZoubM.Al RawiM. (2008). An efficient approach for computing silhouette coefficients. J. Comput. Sci.43, 252–255. 10.3844/jcssp.2008.252.255

  • CrossRef
  • Google Scholar

• 3

  BainL. J. (1978). Statistical analysis of reliability and life-testing models, theory and methods. New York: Marcel Dekker.

  • Google Scholar

• 4

  BebbingtonM. S.CroninS. (2011). Spatio-temporal hazard estimation in the auckland volcanic field, New Zealand, with a new event-order model. Bull. Volcanol.73, 55–72. 10.1007/s00445-010-0403-6

  • CrossRef
  • Google Scholar

• 5

  BebbingtonM. S. (2008). Incorporating the eruptive history in a stochastic model for volcanic eruptions. J. Volcanol. Geoth. Res.175, 325–333. 10.1016/j.jvolgeores.2008.03.013

  • CrossRef
  • Google Scholar

• 6

  BebbingtonM. S.LaiC. D. (1996). On nonhomogeneous models for volcanic eruptions. Math. Geol.28, 585–600.

  • Google Scholar

• 7

  BebbingtonM. S. (2012). Models for temporal volcanic hazard. Statistics in Volcanology1, 1–24. 10.5038/2163-338X.1.1.

  • Google Scholar

• 8

  BebbingtonM. S. (2010). Trends and clustering in the onsets of volcanic eruptions. J. Geophys. Res.115, B01203. 10.1029/2009JB006581

  • CrossRef
  • Google Scholar

• 9

  BevilacquaA.FlandoliF.NeriA.IsaiaR.VitaleS. (2016). Temporal models for the episodic volcanism of campi flegrei caldera (Italy) with uncertainty quantification. J. Geophys. Res. Solid Earth121, 7821–7845. 10.1002/2016JB013171

  • CrossRef
  • Google Scholar

• 10

  BolshakovaN.AzuajeF. (2003). Cluster validation techniques for genome expression data. Signal Process.83, 825–833. 10.1016/S0165-1684(02)00475-9

  • CrossRef
  • Google Scholar

• 11

  BonadonnaC.CostaA. (2013). Plume height, volume, and classification of explosive volcanic eruptions based on the weibull function. Bull. Volcanol.75, 742. 10.1007/s00445-013-0742-1

  • CrossRef
  • Google Scholar

• 12

  BurtM. L.WadgeG.ScottW. A. (1994). Simple stochastic modelling of the eruption history of a basaltic volcano: nyamuragira, Zaire. Bull. Volcanol.56, 87–97. 10.1007/BF00304104.

  • Google Scholar

• 13

  CalvacheM. L.CortésG. P.WilliamsS. N. (1997). Stratigraphy and chronology of the Galeras volcanic complex, Colombia. J. Volcanol. Geotherm. Res.77, 5–19. 10.1016/S0377-0273(96)00083-2

  • CrossRef
  • Google Scholar

• 14

  CalvacheV. M. L. (1990). Geology and volcanology of the recent evolution of Galeras Volcano, Colombia. Master’s thesis. Baton Rouge: Louisiana State University.

  • Google Scholar

• 15

  CepedaH. (2020). Quinientos años de documentación histórica de actividad del volcán Galeras: escenarios eruptivos (Servicio Geológico Colombiaon).

  • Google Scholar

• 16

  ChenG.JaradatS. A.BanerjeeN.TanakaT. S.KoM. S. H.ZhangM. Q. (2002). Evaluation and comparison of clustering algorithms in analyzing es cell gene expression data. Stat. Sin.241–262.

  • Google Scholar

• 17

  ColesS.SparksR. (2006). “Extreme value methods for modeling historical series of large volcanic magnitudes,” in Statistics in Volcanology. Editors MaderH.ColesS.ConnorC.ConnorL. (London: Geological Society IAVCEI Publications), 47–56.

  • Google Scholar

• 18

  ConnorC. B.SparksR. S. J.MasonR. M.BonadonnaC. (2003). Exploring links between physical and probabilistic models of volcanic eruptions: the soufriere hills volcano, Montserrat. Geophys. Res. Lett.30, 1701. 10.1029/2003GL017384

  • CrossRef
  • Google Scholar

• 19

  ConwayF. M.ConnorC. B.HillB. E.ConditC. D.MullaneyK.HallC. M. (1998). Recurrence rates of basaltic volcanism in sp cluster, san francisco volcanic field, Arizona. Geology26, 655–658. 10.1130/0091-7613(1998)026 (0655:RROBVI) 2.3.CO;2

  • CrossRef
  • Google Scholar

• 20

  CroninS.BebbingtonM.LaiC. (2001). A probabilistic assessment of eruption recurrence on taveuni volcano, Fiji. Bull. Volcanol.63, 274–288. 10.1007/s004450100144

  • CrossRef
  • Google Scholar

• 21

  DANE (2011). Estimación y proyección de población nacional, departamental y municipal total por área 1985-2020. Tech. Rep.Bogotá: Departamento Administrativo Nacional de Estadística.

  • Google Scholar

• 22

  De la Cruz-ReynaS. (1991). Poisson-distributed patterns of explosive eruptive activity. Bull. Volcanol.54, 57–67.

  • Google Scholar

• 23

  DeligneN.ColesS.SparksR. (2010). Recurrence rates of large explosive volcanic eruptions. J. Geophys. Res.115, B06203. 10.1029/2009JB006554

  • CrossRef
  • Google Scholar

• 24

  EspinozaB. A. (1988). Actividad de1 volcán Galeras en épocas históricas. Tech. Rep.80 [in Spanish].

  • Google Scholar

• 25

  FurlanC. (2010). Extreme value methods for modelling historical series of large volcanic magnitudes. Stat. Model. Int. J.10, 113–132. 10.1177/1471082X0801000201

  • CrossRef
  • Google Scholar

• 26

  Garcia-AristizabalA.BucchignaniE.ManziM. P. (2017). Patterns in climate-related parameters as proxy for rainfall deficiency and aridity: application to Burkina Faso. ASCE-ASME J. Risk Uncertain. Eng. Syst. A.3, A4016001. 10.1061/AJRUA6.0000860

  • CrossRef
  • Google Scholar

• 27

  Garcia-AristizabalA.MarzocchiW.FujitaE. (2012). A Brownian model for recurrent volcanic eruptions: an application to Miyakejima volcano (Japan). Bull. Volcanol.74, 545–558. 10.1007/s00445-011-0542-4

  • CrossRef
  • Google Scholar

• 28

  Garcia-AristizabalA.SelvaJ.FujitaE. (2013). Integration of stochastic models for long-term eruption forecasting into a Bayesian event tree scheme: a basis method to estimate the probability of volcanic unrest. Bull. Volcanol.75, 1–13. 10.1007/s00445-013-0689-2

  • CrossRef
  • Google Scholar

• 29

  Global Volcanism Program [Dataset] (2013). Version. 4.9.1. “Galeras (351080),” in Volcanoes of the World. Editor VenzkeE. (Washington, DC: Smithsonian Institution).

  • Google Scholar

• 30

  GutenbergR.RichterC. (1944). Frequency of earthquakes in California. Bull. Seismol. Soc. Am.34, 185–188.

  • Google Scholar

• 31

  Hurtado ArtunduagaA. D.Cortés JiménezG. P. (1997). Third version of the hazard map of Galeras Volcano, Colombia. J. Volcanol. Geoth. Res.77, 89–100. 10.1016/s0377-0273(96)00088-1

  • CrossRef
  • Google Scholar

• 32

  INGEOMINAS (2010). Technical Reports. Resumen de actividad del volcan galeras del 29 de diciembre de 2009 al 4 de enero de 2010. Observatorio vulcanológico Y sismológico de pasto resumen. INGEOMINAS.

  • Google Scholar

• 33

  INGEOMINAS (2005). Technical Reports. Volcán Galeras Actividad histórica. INSTITUTO COLOMBIANO DE GEOLOGIA Y MINERIA INGEOMINAS.

  • Google Scholar

• 34

  KaufmanL.RousseeuwP. J. (2008). Finding Groups in Data: An Introduction to Cluster Analysis. John Wiley & Sons.

  • Google Scholar

• 35

  KleinF. (1982). Patterns of historical eruptions at Hawaiian volcanoes. J. Volcanol. Geoth. Res.12, 1–35. 10.1016/0377-0273(82)90002-6

  • CrossRef
  • Google Scholar

• 36

  LillieforsH. (1969). On the Kolmogorov-smirnov test for the exponential distribution with mean unknown. J. Am. Stat. Assoc.64, 387–389.

  • Google Scholar

• 37

  MacQueenJ. (1967). “Some methods for classification and analysis of multivariate observations,” in Proceedings of the fifth berkeley symposium on mathematical statistics and probability. Berkeley, California: University of California Press, Vol. 1, 281–297.

  • Google Scholar

• 38

  MarzocchiW.SandriL. (2003). A review and new insights on the estimation of the b -value and its uncertainty. Ann. Geophys.46, 1271–1282. 10.1002/2016JB013054

  • CrossRef
  • Google Scholar

• 39

  MarzocchiW.SandriL.HeuretA.FunicielloF. (2016). Where giant earthquakes may come. J. Geophys. Res. Solid Earth121, 7322–7336. 10.1002/2016JB013054

  • CrossRef
  • Google Scholar

• 40

  MarzocchiW.ZaccarelliL. (2006). A quantitative model for time-size distribution of eruptions. J. Geophys. Res.111, B04204. 10.1029/2005JB003709

  • CrossRef
  • Google Scholar

• 41

  MasonB.PyleD.OppenheimerC. (2004). The size and frequency of the largest explosive eruptions on earth. Bull. Volcanol.66, 735–748. 10.1007/s00445-004-0355-9

  • CrossRef
  • Google Scholar

• 42

  MassaroS.CostaA.SulpizioR.CoppolaD.CapraL. (2019). Cyclic activity of the fuego de colima volcano (Mexico): insights from satellite thermal data and nonlinear models. Solid Earth10, 1429–1450. 10.5194/se-10-1429-2019

  • CrossRef
  • Google Scholar

• 43

  MelnikO.CostaA. (2014). “Dual chamber-conduit models of non-linear dynamics behaviour at soufrière hills volcano,” in The eruption of Soufriere Hills volcano, Montserrat from 2000 to 2010. Editors WadgeG.RobertsonR.VoightB. (London: Geological Society of London), Chap. 3, Vol. 39 of 5, 61–69.

  • Google Scholar

• 44

  Mendoza-RosasA.la Cruz-ReynaS. D. (2008). A statistical method linking geological and historical eruption time series for volcanic hazard estimations: applications to active polygenetic volcanoes. J. Volcanol. Geoth. Res.176, 277–290. 10.1016/j.jvolgeores.2008.04.005.

  • Google Scholar

• 45

  MulargiaF.TintiS.BoschiE. (1985). A statistical analysis of flank eruptions on Etna volcano. J. Volcanol. Geoth. Res.23, 263–272.

  • Google Scholar

• 46

  NewhallC.SelfS. (1982). The Volcanic Explosivity Index (VEI): an estimate of explosive magnitude for historical volcanism. J. Geophys. Res.87, 1231–1238.

  • Google Scholar

• 47

  PapaleP. (2018). Global time-size distribution of volcanic eruptions on earth. Sci. Rep.8, 6838. 10.1038/s41598-018-25286-y

  • CrossRef
  • Google Scholar

• 48

  PassarelliL.SansòB.SandriL.MarzocchiW. (2010). Testing forecasts of a new bayesian time-predictable model of eruption occurrence. J. Volcanol. Geotherm. Res.198, 57–75. 10.1016/j.jvolgeores.2010.08.011

  • CrossRef
  • Google Scholar

• 49

  PyleD. (1998). Forecasting sizes and repose times of future extreme volcanic events. Geology26, 367–370. 10.1130/0091-7613(1998)026(0367:FSARTO) 2.3.CO;2

  • CrossRef
  • Google Scholar

• 50

  PyleD. (2000). “Sizes of volcanic eruptions,” in Encyclopedia of volcanoes. Editors SigurdssonH.HoughtonB.McNuttS.RymerH.StixJ. (London: Academic Press).

  • Google Scholar

• 51

  RousseeuwP. J. (1987). Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. J. Comput. Appl. Math.20, 53–65. 10.1016/0377-0427(87)90125-7

  • CrossRef
  • Google Scholar

• 52

  SandriL.AcocellaV.NewhallC. (2017). Searching for patterns in caldera unrest. Geochem. Geophys. Geosyst.18, 2748–2768. 10.1002/2017GC006870

  • CrossRef
  • Google Scholar

• 53

  SandriL.MarzocchiW.GasperiniP. (2005). Some insights on the occurrence of recent volcanic eruptions of Mount Etna volcano (Sicily, Italy). Geophys. J. Int.163, 1203–1218. 10.1111/j.1365-246X.2005.02757.x.

  • Google Scholar

• 54

  SGC (2015). Actualización del mapa de amenaza volcánica del volcán galeras—Colombia. Pasto: Servicio geológico Colombiano, Dirección de Geoamenazas.

  • Google Scholar

• 55

  SimkinT.SiebertL. (2002). GVP-3. Volcanoes of the world: an illustrated catalog of holocene volcanoes and their eruptions. Smithsonian Institution, Global Volcanism Program Digital Information SeriesAvailable at: http://www.volcano.si.edu/world/ (Accessed 2020).

  • Google Scholar

• 56

  SobradeloR.MartiJ.Mendoza-RosasA.GomezG. (2011). Volcanic hazard assessment for the canary islands (Spain) using extreme value theory. Nat. Haz. Earth Syst. Sci.11, 2741–2753. 10.5194/nhess-11-2741-2011

  • CrossRef
  • Google Scholar

• 57

  StixJ.CalvacheV. M. L.WilliamsS. N. (1997). Galeras volcano, Colombia interdisciplinary study of a decade volcano. J. Volcanol. Geotherm. Res.77, 1–4. 10.1016/s0377-0273(96)00082-0

  • CrossRef
  • Google Scholar

• 58

  TurnerM.CroninS.BebbingtonM.PlatzT. (2008). Developing probabilistic eruption forecasts for dormant volcanoes: a case study from mt taranaki, New Zealand. Bull. Volcanol.70, 507–515. 10.1007/s00445-007-0151-4

  • CrossRef
  • Google Scholar

• 59

  UNGRD (2019). Fortalecimiento de la gestión del riesgo de desastres en la zona de amenaza volcánica alta-zava del volcán galeras. (Bogotá, Colombia: Presidencia de la República).

  • Google Scholar

• 60

  VoightB. (1989). A relation to describe rate-dependent material failure. Science.243, 200–203. 10.1126/science.243.4888.200

  • CrossRef
  • Google Scholar

• 61

  WickmanF. E. (1966). “Markov models of repose-period patterns of volcanoes,” in Random process in geology. New York: Springer, 135–161.

  • Google Scholar

• 62

  XuR.WunschD. (2005). Survey of clustering algorithms. IEEE Trans. Neural Network16, 645–678. 10.1109/TNN.2005.845141

  • CrossRef
  • Google Scholar