## Summary

In the first part of this work, we make use of two non-parametric statistical pattern recognition algorithms and a multiple regression analysis to analyse seismic clusters occurring around Mount Etna, Italy. The aim is to determine if the onset of flank eruptions at Mount Etna is linked to variations in the regional seismicity at a timescale of few weeks. From the analysis, we find that the discrimination between clusters preceding flank eruptions and clusters not related in time to flank activity is mainly linked to the volume output of the previous flank eruption, in some cases together with the time elapsed from its end. Instead, we do not find any difference in the seismicity features characterizing different types of clusters, except for a very small contribution of the number of seismic events in the clusters. This result does not confirm the existence, suggested in the past, of a direct link between the regional state of stress at a timescale of few weeks and the occurrence of flank eruptions on Mount Etna volcano. On the contrary, the result suggests that a prominent role in the flank eruption occurrence is played by the re-charging of the feeding system. In the second part of this study we analyse the relationship between the magma volume erupted in an eruption and the interevent time following it, finding that a ‘time-predictable model’ satisfactorily describes the occurrence of eruptions at Mount Etna in the last decades. The latter analysis is carried out both on the flank eruption catalogue only, and on the complete catalogue of flank and summit eruptions, with comparable results.

## Introduction

Mount Etna volcano is one of the most extensively monitored volcanic systems of the world. Besides an almost continuous activity in the summit craters, every few years Mount Etna experiences episodes of flank activity. Due to its densely populated slopes, the Etna flank eruptions represent one of the most prominent source of volcanic damage in that area. The first step in estimating this kind of volcanic hazard is to provide a reliable model for the temporal occurrence of flank eruptions. Note that, in this work, we do not take into account the spatial density of flank eruptions, which has been previously studied, for example, by Wadge *et al.* (1994). Regarding the temporal aspect, previous studies provide some empirical constraints by studying the relationship between flank eruption occurrence and regional and local seismicity at a timescale of few weeks (Mulargia *et al.* 1991, 1992; Vinciguerra *et al.* 2001). In these studies, seismicity is assumed to be a proxy for the local and regional tectonics, and therefore is not related to magma motion beneath the volcano. Despite some differences in the results, all of these papers suggest a tectonic control on the occurrence of flank eruptions at Mount Etna, at least over a timescale of few weeks.

Alternatively, flank eruptions have been characterized studying the statistical distribution of the interevent times and sizes (Wadge & Guest 1981; Mulargia *et al.* 1987a; Gasperini *et al.* 1990), but no univocal conclusions were achieved. While these studies generally agree on a steady magma output in recent years (the periods studied are 1971–1981 by Wadge & Guest 1981; Mulargia *et al.* 1987a and 1978–1987 by Gasperini *et al.* 1990), they identify several, and different, change points in the interevent time distributions. Remarkably, they do not find any relationship between the size of a flank eruption and the subsequent interevent time (*cf.*Mulargia *et al.* 1987a).

The organization of this work is mainly determined by the results obtained; in particular we started to investigate on a specific hypothesis, and the results obtained led us to look also into other directions.

At first, we re-address the link between flank eruptions and regional/local seismicity by means of two non-parametric statistical pattern recognition (NSPR) algorithms based on very different approaches, and by means of a multivariate regression analysis. The statistical pattern recognition algorithms used have been previously tested on synthetic data that mimic the basic properties of our data set of seismic clusters (Sandri & Marzocchi 2004). Also, the original seismic and volcanic catalogues available for this study are much more reliable and complete than those used in past works (*cf.*Mulargia *et al.* 1991, 1992).

We analyse groups of seismic events that we call ‘seismic clusters’(see definition in the next section) occurring in a large area around Mount Etna. The goal of this part of the study is to check whether there is a recurrent pattern in this regional seismicity preceding a flank eruption at Mount Etna, by considering seismicity as a proxy for regional tectonics (Lay & Wallace 1995). Note that the low?\average magnitude of the earthquakes considered here (see next section) do not seem indicative of a direct causal relationship between earthquake occurrence and the triggering or promoting of flank eruptions, which has been demonstrated for much larger seismic events (e.g. Marzocchi 2002; Marzocchi *et al.* 2004).

Through multivariate regression analysis we study how the characteristics of the seismic clusters vary as the next flank eruption is approaching, without grouping the data into classes.

We test the validity of the patterns found by applying them to an independent data set.

The second part of the work is mainly motivated by the results obtained in the previous point, that show the prominent importance of the re-charge of the feeding system, apparently independent of the regional tectonic stress (at least over the timescale investigated). In particular, in the second part of the work we study in detail the relationship between interevent time (defined as the time between the onset of two successive eruptions) and size of eruptions, providing some empirical constraints for the mechanism that generates eruptions.

## The Data

We analyse two different catalogues: one for seismic events and one for flank eruptions of Mount Etna.

Regarding the seismic data, similarly to Mulargia *et al.* (1991), we take into account seismic events that occurred in a circle of radius 120 km, centred over the summit of Mount Etna (Fig. 1). The catalogue from which we extract the seismic events consists of the combination of different data sets. From 1981 to 1996 we make use of the ‘Catalogo Strumentale dei Terremoti Italiani dal 1981 al 1996 (CSTI)’(CSTI Working Group, 2001) where the instrumental local magnitudes have been re-valuated, according to Gasperini (2002), by the comparison with a data set of true and synthetic Wood?CAnderson estimates. We actually use a slightly improved version of CSTI where the hypocenter depths of events having large azimuthal gaps (>180°) are fixed to 10 km as well as the magnitude estimates associated with large *P*-wave residuals (>5 s.) are discarded. From 1960 to 1980, we essentially refer to the catalogue of the ‘Progetto Finalizzato Geodinamica (PFG)’(Postpischl 1985) but applying an empirical magnitude shift of −0.3 units from 1960 to mid-1976 and of −0.5 units from mid-976 to 1980 in order to make it homogeneous with the re-calibrated section described above (see Lolli & Gasperini 2003). Finally, for the earthquakes since 1997, we use the hypocentral locations reported by the INGV Seismic Bulletin (1997–2002), with magnitudes recomputed from duration and amplitude data according to Gasperini (2002).

Before analysing any seismic catalogue, it is necessary to establish its completeness for the magnitude threshold we are interested (see e.g. Albarello *et al.* 2001). We evaluate the completeness of the merged catalogue obtained by visually inspecting the cumulative number of earthquakes (Mulargia *et al.* 1987b) as a function of time for two different magnitude thresholds: *M _{c}* = 3.0 (Fig. 2a) and

*M*= 2.5 (Fig. 2b). The completeness can be assumed since 1974 and since 1983, respectively, although for the latter threshold it is actually possible to distinguish a slight intermediate change in the slope (Fig. 2b), roughly between 1983 and 1988. We decide not to consider this minor effect, because it is much less marked than the change in slope since 1983. In the following we will consider both data sets, as the former is longer, but with fewer events, while the latter is more detailed but only on the most recent period.

_{c}The seismic events extracted for the two data sets are grouped to obtain seismic clusters lasting at most 1 week. In particular, the first cluster consists of the first seismic event in the considered area and of all the seismic events occurred in the same day and in the following 6 days. The other clusters are formed, in the same way, starting from the first event following each cluster. We remark that we are in a volcanic area, without clear cluster sequences in the usual seismological meaning, and that could be defined using clustering algorithms such as the one proposed by Reasenberg (1994). The definition of our clusters is instead a simple grouping of earthquakes occurring in a time window of a selected length. The spatial distribution of such clusters will be considered in the analysis and discussed later.

We obtain, in total, 218 seismic clusters for the first data set (since 1974, *M _{c}* = 3.0, see Fig. 1a) and 389 for the second one (since 1983,

*M*= 2.5, see Fig. 1b).

_{c}Concerning the flank eruption data, we refer to Behncke *et al.* (2005) for the complete list of eruptions at Mount Etna. However, we first have to define what flank eruptions at Mount Etna are, and then extract them from this list. We borrow the definition of flank eruptions at Mount Etna from Gasperini *et al.* (1990), according to which flank eruptions are those linked to the opening of magma??filled fractures on the slopes of Mount Etna. In this respect, we define a maximum altitude (2900–3000 m asl) above which the rocks at Mount Etna are so incoherent and unconsolidated that cannot fracture. As a consequence, all the eruptions occurring below this height are considered flank eruptions, and are listed in Table 1. The time-series of the seismic clusters and flank eruptions in the Etnean area are displayed in Fig. 3 superimposed to the sequence of the clusters detected for both the above data sets. All but one the flank eruptions are preceded by at least a seismic cluster in the 100 days preceding the onset of the volcanic event. The only flank eruption occurring without precursory seismicity in a period of 100 days is the one starting on 1975 November 29.

## Searching for Recurrent Patterns in the Seismicity Preceding Flank Eruptions

The primary goal of this section is not to predict flank eruptions, which might be better forecast by tracking magma motion beneath the volcano by means of microseismicity. In this part of the paper we rather look for recurrent patterns in local and regional seismicity preceding flank eruptions. Considering such seismicity as a proxy for regional and local tectonics (Lay & Wallace 1995), we therefore try to provide some insights on the possible tectonic control of flank eruptions over a specified time interval.

For this purpose, we perform two multivariate statistical analyses by means of two different methods.

### Pattern recognition analysis

Pattern recognition (PR) is a powerful multivariate analysis technique allowing, in principle, the identification of possible repetitive schemes among the objects belonging to distinct classes. While usual data analysis takes into account only one variable of the process at a time, PR is able to extract information from any possible combination (linear or not) of variables that are suspected to have an influence on the process. Moreover, PR does not need the construction of a theoretical model, but is usually based on a sole basic hypothesis, that is, the assumption that the phenomenon under study is governed by a finite number of complex, but repetitive patterns of the variables.

These appealing properties led PR to be applied with success in several and diverse disciplines, which share the study of complex systems (e.g. Fukunaga 1990). For the same reasons, we believe that PR, and especially NSPR, might be a very promising tool also in Earth Sciences. Until now the only few remarkable efforts in this direction are M8 and CN algorithms (Keilis-Borok & Kossobokov 1990; Keilis-Borok *et al.* 1988), and applications to volcanology (Mulargia *et al.* 1991, 1992; Vinciguerra *et al.* 2001; Sandri 2004). With respect to these applications, our study represents an improvement because we make use of algorithms that have been previously tested to verify their performance on data sets with similar properties to the real ones (Sandri & Marzocchi 2004; Sandri *et al.* 2004).

Technically, the main goal of PR is to classify *objects*. Every object is represented by an array of qualitative or quantitative variables. The procedure of analysis consists of three different steps: the learning phase, the voting phase, and the control experiments. In the learning phase, a set of known and classified objects is analysed in order to recognize all the possible *patterns* that characterize each class, that is, the combinations of variables that allow to discriminate the objects belonging to different classes. This step turns out to be very useful because it leads to the recognition of the most important variables for the process under study. In the voting phase, the patterns identified during the learning phase are used to classify new objects, whose class is unknown to the algorithm. Finally, control experiments check the stability of the results by repeating the learning and voting phases using different values for the algorithm parameters.

We use two statistical NSPR algorithms named binary decision tree (BDT) and Fisher analysis (FIS). Both of these algorithms have been previously tested on synthetic data (Sandri & Marzocchi 2004). Such simulations have shown that both BDT and FIS provide reliable results when applied to data sets that mimic the basic features of our seismic cluster data set, that is, few objects, each represented by many features, often not normally distributed. For a more complete definition of the algorithms, see Appendices A and B, while, for a very detailed description, see Duda & Hart (1973) for FIS, and Rounds (1980), Mulargia *et al.* (1992) for BDT.

In this paper we perform the learning phase both on the totality of available data and on a randomly chosen subset (consisting of 80 per cent of the available data). In the former case we attempt to recognize, as a first step, all the possible patterns in our data set. In the latter case, we reserve part of the available data for the voting phase.

Before performing the learning phase, we first have to define the objects to be analysed and the classes involved in the problem. In our study, the objects are the seismic clusters. Any object is represented by a vector that contains all the measurements (the features) that we can associate to the object. For each seismic cluster we retrieved the following measurements:

- 1
maximum magnitude recorded in the cluster (MXM),

- 2
number of events in the cluster (NEV),

- 3
trimmed mean 40 per cent (i.e. the mean on the data between the 20th and the 80th percentiles) of the latitudes of the events in the cluster (LAT),

- 4
trimmed mean 40 per cent of the longitudes of the events in the cluster (LON),

- 5
season of occurrence of the cluster (SEA),

- 6
time elapsed (in days) from the previous cluster (TPC),

- 7
time elapsed (in days) from the end of the previous flank eruption (also called time-since-eruption, TSE),

- 8
maximum magnitude recorded in the previous cluster (MXP) and

- 9
volume output of the previous flank eruption in millions of

*m*^{3}(VOP).

The first five measurements (MXM, NEV, LAT, LON and SEA) are related to the cluster itself; in particular, LAT and LON describe the spatial distribution of the clusters. The latter four (TPC, TSE, MXP and VOP) are related to the recent volcanic and seismic history of the system. Due to the large differences between the maximum and minimum measurements in the catalogue for the TSE, for the time elapsed from the previous clusters (TPC) and for the volume output of the previous eruption (VOP), we decide to use the logarithm of these values. Thus, each object has the following components: MXM, NEV, LAT, LON, SEA, Log(TPC), Log(TSE), MXP and Log(VOP).

For the sake of clarity, from now on we call ‘time-to-eruption’ the time to the following flank eruption, and, as mentioned above, ‘time-since-eruption’ the time between the end of last flank eruption and the beginning of the cluster (i.e. TSE). We identify three classes of clusters:

where τ is the time window for defining precursory clusters (see below). Type ‘A’ clusters are precursors to flank eruptions, while type ‘B’ occur right after the end of a flank eruption. Type ‘C’ clusters are not temporally associated with flank eruptions (Fig. 4).If both inequalities (1) and (2) are verified, then the cluster is not considered in the analysis, as well as clusters occurring during flank eruptions.

Clearly, the attribution of an object to one of these types depends on the parameter τ. We arbitrarily define that a cluster has to be classified as ‘A’ if it occurs at most about 3 months before the following flank eruption, that is, τ = 100 days. This choice is related to the temporal scale of the tectonic interactions we are investigating. The timescales considered in our NSPR analysis (weeks to few months) are typical of the elastic response of the lithosphere to irreversible perturbations. An example of mechanism occurring at these temporal scales is the typical aftershock occurrence following the irreversible deformation due to the mainshock (Dieterich 1994). Hereinafter, we use the term ‘quasi-elastic’ to indicate the time, and spatial, scale used, being too short for the viscoelastic responses (from years to decades, e.g. Piersanti 1999; Kenner & Segall 2000).

With the selected value for τ, about 12 per cent of the clusters contained in our data set are of type ‘A’ In the control experiments, we will check the stability of our results to changes of the value of τ.

By means of such definitions (conditions 1–3) we have, in principle, a three-class problem. Since we are mainly interested to a possible regional tectonic stress pattern favouring flank activity at Mount Etna, we take into account only the precursory clusters (type ‘A’ from now on representing class 1) and the clusters occurring ‘away’ from flank eruptions (type ‘C’ from now on representing class 2). In fact, most of the activity immediately after the onset and during a flank eruption may be related to a stress re-distribution induced by the eruption itself.

### Results of the pattern recognition analysis

We define ‘repose’ time the time between the end of a flank eruption and the onset of the following one. As a flank eruption is approaching, the number of seismic clusters does not increase. In fact, as shown in Fig. 5, the frequency of the seismic clusters recorded in the data set is approximately independent on the time-to-eruption. In Fig. 5, the frequency is normalized, that is, for each value of tim-to-eruption *t* on the *x*-axis, the correspondent frequency *f*(*t*) is divided by the number of repose times, contained in our flank eruption data set, larger or equal to *t*. This is done to take into account the fact that a large time-to-eruption is possible only if the repose time is long enough, that is, at least longer than the time-to-eruption. If we do not normalize, we obviously obtain a much larger number of clusters for short times-to-the-eruption, because such short times are possible for any repose time, even for short ones. As we see in Fig. 5, for the case since 1983 (Fig. 5b) the normalized frequency of short times-to-eruption is even lower than the average frequency in the data set.

The results of the NSPR analysis, when all the data are used in the learning phase, are shown in Figs 6 and 7 (for BDT) and Figs 8 and 9 (for FIS). For both algorithms and in both cases (since 1974, *M _{c}* = 3.0; since 1983,

*M*= 2.5) the most important feature is the logarithm of the volume output [Log(VOP)] of the previous flank eruption. Also the logarithm of the time from the previous flank eruption [Log(TSE)] is important. In particular, the pattern found based on both features tells that, for a seismic cluster occurring around Mount Etna, the larger the VOP, the longer the TSE needed in order to consider the cluster as a precursory one (class 1 or type ‘A’. However, the most striking result is that no seismic feature is found with discriminating power between the two classes of seismic clusters.

_{c}The clusters of type ‘A’ correctly classified by BDT are 95 per cent and 78 per cent, respectively, for the two cases (since 1974 *M _{c}* = 3.0, and since 1983

*M*= 2.5), while those of type ‘C’ are 56 per cent and 64 per cent. The results obtained by FIS algorithms are 64 per cent, 67 per cent, 83 per cent and 86 per cent, respectively, showing that FIS performs worse than BDT on class ‘A’ objects, but better on type ‘C’

_{c}By using a randomly selected subset consisting of 80 per cent of the available data for the learning phase, and voting the remaining 20 per cent of data as new and independent objects, we obtain the results displayed in Table 2. In particular, we notice that:

- (1)
the two algorithms provide consistent results;

- (2)
the pattern found, based on Log(TSE) and Log(VOP), involves the same features as in the case when all the data are used for the learning step, but it is more stable;

- (3)
except for BDT in the case since 1974 (

*M*= 3.0), the fraction of learning objects correctly classified is comparable to the one obtained when classifying the voting data, that is, the risk of over?\fitting the data can be excluded._{c}

We perform some control experiments to check the stability of the results obtained. In particular, we repeat four times the NSPR analysis, each time considering one of the following variants: The control experiments (see Table 3) overall confirm the predominant importance of Log(VOP); furthermore, Log(TSE) and NEV show a recurrent discriminant ability. No other recurrent pattern is found.

- 1
a different value of the parameter τ, by using τ = 40 days;

- 2
a different cluster duration, by re-extracting clusters of one month (instead of 1 week), from the original seismic catalogues;

- 3
a different flank eruption catalogue, by excluding all the cluster following the largest eruption listed (the one starting in 1991, see Table 1). The reason why we chose not to consider these clusters is to check whether the pattern found (based on the logarithms of the TSE and of the volume output of the previous lateral event) can be due only to these clusters.

- 4
a different seismic catalogue, by including only the seismic events in a circle of 30 Km of radius centred on Mount Etna. The reason why we do this control experiment is to check whether the trimmed means of the coordinates are biased by the simultaneous occurrence of seismic events close and far away from the volcano.

Even though the results obtained in these four control cases are less stable than those obtained in the standard analysis, they confirm that no seismicity feature has a predominant discriminating ability in distinguishing between precursory clusters and clusters not related in time-to-flank activity. Remarkably, our result is different from those previously obtained by Mulargia *et al.* (1991, 1992) and Vinciguerra *et al.* (2001). All of these studies, in fact, recognized particular seismicity patterns preceding flank activity, linked, in some way, to the tectonic regime

### Multivariate regression fit

In this section we analyse the seismic cluster data set with a multivariate regression analysis (Draper & Smith 1998). The use of this statistical technique is particularly indicated in the present case where the TSE seems to be a relevant feature. In fact, in NSPR analysis the continuous variables [for instance, Log(TSE)] are necessarily grouped by the algorithms in order to establish the boundary between the two classes. This grouping does not allow distinguishing between different values of the variable inside the same class; for example, very similar times-from-eruption can belong to different classes if they fall on different sides of the threshold value, while, on the contrary, very different times-from-eruption can belong to the same class.

In this case, the use of multivariate regression analysis overcomes this possible drawback because it finds out how the time-to-eruption variable depends on the other variables considered without imposing any grouping. At the same time, the method allows checking the stability of the results obtained with the NSPR analysis, and to further investigate on the pattern found based on Log(VOP) and Log(TSE).

Before going through the details of the analysis, it is worth to make some cautionary remarks on this kind of analysis. In particular, when we do regression calculations on unplanned data (that is, data arising from continuing operations and not from a designed experiment), some potentially dangerous possibilities can arise. For example, a false effect (a bias) on a visible variable may be caused by an unmeasured latent variable. Another undesired effect is linked to the case in which the most effective variable is kept within quite a small range, that might lead to the misinterpretation that the variable has no significance. A third problem is that the use of unplanned data often causes large correlations between predictors; this makes it impossible to attribute a causal effect to one specific predictor.

From a qualitative graphical point of view, our multivariate regression analysis builds a linear model in which the dependent variable (*y*) is the logarithm of the time from the beginning of the cluster to the following flank eruption. The independent variables are the nine variables observed for each cluster. The first step of the analysis consists of selecting the relevant variables of these nine that account for the variability of the dependent variable. Here, we select the most relevant variables by using the procedure called ‘best subset search’ described by Garside (1971). The procedure first requires the fitting of every possible regression equation, which involves any combination of independent variables, and then selecting the case, which is the best with respect to one criterion. The criterion adopted here is based on the *R*^{2} coefficient, representing the ratio of the variance explained by the regressive model to the total variance. We report the case with the highest *R*^{2} coefficient obtained by using a single independent variable, two independent variables, three independent variables and so on (e.g. see Table 4). The ‘best’ regression is the one for which the further gain in *R*^{2} is significantly reduced. This method contains a degree of subjectivity (see discussion in Draper & Smith 1998), but straightforwardly identifies the parameters explaining the largest part of variance in the dependent variable.

### Results of the multivariate regression fit

Table 4 and the upper panels of Fig. 10 report the results obtained by using all the available data and, as dependent variable, the logarithm of the time to the next flank eruption. These results show that the ‘best’ regression is the one obtained when using three independent variables, in particular NEV, Log(TSE) (TSE in days) and Log(VOP) (VOP in millions of *m*^{3}), with coefficients:

*M*= 3.0, and for the case since 1983,

_{c}*M*= 2.5. The addition of any other feature improves very little the

_{c}*R*

^{2}coefficient. Note that the regression coefficient for NEV is one order of magnitude smaller than those for Log(TSE) and Log(VOP), and it has the largest percentage error. This means that these latter two parameters explain by far the largest part of the variance, as found in the NSPR analysis.

The results obtained on a randomly chosen subset of the data set available, consisting of only 80 per cent of the total amount of data, are shown in Table 5 and in the lower panels of Fig. 10. Here, again, the best subset of variables consists of Log(TSE), Log(VOP) and NEV.

In Table 6 the results of the multivariate regression analysis obtained in the same control experiments (as in the NSPR analysis) are displayed. In the majority of the experiments, these three parameters explain the largest part of variance in the dependent variable.

In any case, we can never explain more than 50–60 per cent of the variability in the data set (see Tables 4 and 5). This implies that the scatter of the data around the best multivariate regression line is too large to allow to use profitably the multivariate regression line as an efficient forecasting rule.

The multivariate regression fit confirm that, except for a minor influence of feature NEV, the discrimination between seismic clusters occurring just before or away from a flank eruptions is mainly due to variables related to the occurrence of the previous flank eruption. In particular, as mentioned above, the larger the volume erupted in a flank eruption, the longer time is needed before having another flank eruption. This is a necessary condition in ‘time-predictable’ systems, and is related to the re-charging of the feeding system. In the next section, we analyse this issue in depth, in order to check whether flank eruption occurrence at Mount Etna is a time-predictable process.

## Is the Occurrence of Flank Eruptions at Mount Etna Time ‘Predictable’

The pattern found by means of PR and multivariate regression analyses suggests us to check the existence of a time-predictable behaviour in the temporal occurrence of flank eruptions at Mount Etna. Actually, since no seismicity feature has been found as discriminant, we are excluding that the regional tectonic stress is the main feature responsible for the opening of magma-filled fractures along Mount Etna flanks. In this respect, if the re-charging of the magmatic system is the main factor controlling eruption occurrence, the distinction between flank and summit eruptions may be no longer useful. Because of this, in the following we verify the reliability of a time-predictable behaviour both on the flank eruption catalogue (Table 1) and on the joint catalogue of flank and summit eruptions (from now on, flank + summit, reported in Table 7). The latter catalogue has been taken from Behncke *et al.* (2005).

A very simple time-predictable model for volcanic eruptions (*cf.*Bacon 1982; Burt *et al.* 1994; Hill *et al.* 1998) assumes that eruptions always occur when the volume of magma in the storage system reaches a threshold value; if the magma input in the magma storage system is approximately constant (*cf.*Mulargia *et al.* 1987a), but the size of eruptions is a random variable (RV) following some kind of statistical distribution, then we have a time-predictable system with longer/shorter interevent times (defined as the time between the onset of two successive eruptions; Mulargia *et al.* 1987a; Burt *et al.* 1994) following large/small volume output eruptions. In this respect, the pattern that we have found with PR and multivariate regression analyses for flank eruptions at Mount Etna is compatible with a time-predictable behaviour.

In order to verify if the occurrence of flank and/or of flank + summit eruptions at Mount Etna really follows a time-predictable pattern, we need to check whether there is a significant linear relationship between the volume erupted during an eruption (flank or flank + summit) and the interevent time following it (the necessary condition for time-predictability to hold). In fact, for a time-predictable system, the time to the next eruption (next interevent time) is determined by the time required for the magma entering the storage system to reach the eruptive level. In this view (*cf.*Burt *et al.* 1994), we have:

*υ*is the volume erupted during the

_{i}*i*th eruption occurring at time

*t*, and

_{i}*r*is the interevent time between the onset of

_{i}*i*th and the onset of (

*i*+ 1)th eruptions (i.e.

*r*=

_{i}*t*−

_{i+1}*t*). If we cumulate eq. (6) over the

_{i}*N*eruptions occurred, the random noises cancel out. In this case, we obtain a linear plot with

*t*on the x-axis (where

_{i}*t*, as mentioned above, is the onset time of the

_{i}*i*th eruption), and, on the

*y*-axis, the cumulative volume

*υ*erupted during the previous and the present eruptive episodes (i.e.

_{i}*V*= Σ

_{i}

^{i}_{k=1}

*υ*,

_{k}*i*= 1, …,

*N*), as the one shown in Fig. 11 for flank eruptions and Fig. 12 for flank + summit. However, the regression analysis to identify a possible relationship between the interevent times and the erupted volumes must not be carried out on the cumulative curve

*t*versus

_{i}*V*. This is because a basic assumption of the regression analysis is that the data must be independent, and they are clearly not in case of cumulated data. In this respect, we must analyse the simple

_{i}*υ*and

_{i}*r*data. Furthermore, in order to reduce the too high leverage of some of the data (see Draper & Smith 1998), in this work we analyse the log

_{i}*υ*and log

_{i}*r*data, although we also display the plot

_{i}*t*vs

_{i}*V*because it represents the standard view when speaking about time-predictable systems. By means of the logarithmic transformation of variables, the ‘time-predictable’ model by Burt

_{i}*et al.*(1994) is generalized, because a possible significant linear relationship between log

*υ*and log

_{i}*r*would imply a power law between

_{i}*υ*and

_{i}*r*: If the exponent

_{i}*b*is equal to 1, we are in a classical (i.e. as modelled by Burt

*et al.*1994) ‘time-predictable’ system; in the other cases, the time derivation of eq. (7) implies non-stationary re-charging process.

For Mount Etna, we first check whether there have been changes in the eruptive behaviour. By means of the CHPT algorithm (Mulargia & Tinti 1985), we first analyse the flank eruption catalogue, in particular the time-series of:
The third time-series analysed is the one related to the slope of the best fit line in Fig. 13, representing the long?\term magma supply rate (Burt *et al.* 1994). So, we perform a regression analysis on the whole series of log *r _{i}* and log

*υ*(from 1971 to 2002, flank eruption catalogue) obtaining the regression results shown in Fig. 13 and in Table 8. We consider as independent variable the interevent times, their uncertainties being much smaller than those on the erupted volumes. As we see from Table 8, the regression model rejects the null hypothesis

_{i}*b*is the slope of the model, estimated from the data through the least-squares technique) at a 0.01 significance level, by means of the

*F*-test. This means that there is a significant linear relationship between the logarithm of the erupted magma volume and the logarithm of the following interevent time. Furthermore, the estimated slope of the linear relationship ( = 0.9 ± 0.3, see Table 8) is consistent with a classical time-predictable model as the one proposed by Burt

*et al.*(1994). The regression model explains only a part of the variability in the data (

*R*

^{2}= 0.44, see again Table 8).

- 1
the erupted magma volume

*υ*, identifying no change point at 0.05 significance level;_{i} - 2
the interevent time

*r*, identifying no change point at 0.05 significance level;_{i} - 3
the ratio of the erupted magma volume

*υ*to the following interevent time_{i}*r*, identifying no change point at a 0.05 significance level;_{i}

We repeat this analysis on the flank + summit eruption catalogue (Table 7). The results for the searching for change points are that:
So, we perform again a regression analysis on the whole series of log *r _{i}* and log

*υ*from the flank + summit catalogue (from 1971 to 2002) obtaining the regression results shown in Fig. 14 and in Table 9. As we see from Table 9, the regression model rejects the null hypothesis

_{i}*b*is the slope of the model, estimated from the data through the least squares technique) at a 0.01 significance level, by means of the

*F*-test. Again, this means that there is a significant linear relationship between the logarithm of the erupted magma volume and the logarithm of the following interevent time. Again, the estimated slope of the linear relationship ( = 1.1 ± 0.3, see Table 9) is consistent with the value of 1. The regression model explains a smaller part of the variability in the data (

*R*

^{2}= 0.29, see again Table 9).

- 1
in the erupted magma volume time-series

*υ*, we identify one change point at 0.05 significance level on 1995 July 30, but none at a 0.01 significance level;_{i} - 2
the interevent time

*r*, identifying one change point at 0.05 significance level on 1990 January 4, but none at a 0.01 significance level;_{i} - 3
the ratio of the erupted magma volume

*υ*to the following interevent time_{i}*r*, identifying no change point at a 0.05 significance level;_{i}

In order to check the stability of the results, we perform the same analyses by using the time elapsed since the onset of an eruption and the time in the middle of the onset and the end of the previous eruption. In this way, we account for the fact that the past eruptions have different durations, associating the erupted volume to the centre of the eruption interval (defined by the onset and the end of the eruption) rather than to its onset. The results obtained are the same than the ones previously reported.

The time-predictability analysis shows that the occurrence of eruptions at Mount Etna is governed by a time-predictable mechanism, since there is a significant linear relationship between the logarithm of the erupted magma volume and the logarithm of the following interevent time. Over the last 30 yr, this mechanism has been approximately steady, because we detect no change point in the ratio of the erupted magma volume to the following interevent time.

## Discussion and Final Remarks

In this work we first have applied statistical pattern recognition techniques to a data set of seismic clusters registered around Mount Etna, in order to check whether there is a statistically significant link between the regional seismicity and the occurrence of flank eruptions on this volcano at a ‘quasi-elastic’ timescale (i.e. ‘few weeks’.

The results of this analysis, confirmed by the control experiments, have shown that the only recurrent pattern found is linked to the volume output of the previous flank eruption. Sometimes this pattern is associated to the time elapsed from the end of the previous flank eruption, with a minor influence of the number of events registered in the seismic cluster. In particular, for a seismic cluster occurring around Mount Etna, the larger the volume output of the last eruption, the less likely a cluster is a precursory one. When the time elapsed from the end of the previous flank eruption is a discriminating feature, the larger the volume output of the last eruption, the longer we need to wait after the end of the eruption before considering the cluster a precursory one.

We have analysed the same seismic cluster data set with a multivariate regression analysis, confirming the recurrence of the pattern just mentioned.

The correct interpretation of these results deserve a careful discussion of the role of seismicity. We know that there are clear patterns of strongly localized seismicity in the volcano edifice related to its gravitational spreading (Borgia *et al.* 1992) and to magma movements in the plumbing system at depth (e.g. Patané *et al.* 2003). From our perspective, and according to Patané *et al.* (2003), this kind of seismicity is not a proxy for tectonics, but rather the response of the volcanic edifice to the forces exerted by gravitational spreading and magma motion, that does not alter the regional seismic rate and the regional seismicity features. The seismicity considered here as a proxy for regional tectonics has a broad spatial distribution and a larger magnitude average. Under this perspective, our results stand for a not significant link between variations in the regional state of stress and the occurrence of flank eruptions of Mount Etna volcano at a ‘quasi-elastic’ timescale.

At the same time, we have noted that the pattern found could be explained by giving a prominent importance to the process governing the re-charging of the feeding system. In order to make this point clearer, we have analysed the volcanic catalogue of the last decades. We have carried out two separate analysis, one considering only flank eruptions, and one considering every eruption (flank and/or summit). The latter type of analysis has been carried out because, assuming that the re-charging of the feeding system plays a predominant role, there is in principle no reason to distinguish between the different ways in which the feeding system reaches the surface and exhibits volcanic activity. By means of CHPT algorithm (Mulargia & Tinti 1985), we have found no significant change point in the ratio of the erupted magma volume to the following interevent time. Then, we have found that a time-predictable model is a reliable first-order mechanism for flank + summit eruptions mainly governed by the feeding of the system. In practice, the interevent time between two subsequent episodes is mainly related to the amount of magma erupted in the first one of the two: the larger the erupted volume, the longer the interevent time.

### Acknowledgments

We wish to thank Dr Marco Neri of Istituto Nazionale di Geofisica e Vulcanologia in Catania, Italy, for the information on the catalogue of Mount Etna flank eruptions and for keeping it updated as the new flank episodes occurred during the writing of the present paper. We also wish to thank Prof. Stefano Gresta for useful revision comments.

### Appendix

#### APPENDIX A: Binary Decision Tree

This method has been developed by Rounds (1980) and, slightly modified, successfully applied to volcanic data by Mulargia *et al.* (1992). It can be used only in the two-class problem and it was originally designed for hierarchically ordered data sets, even though tests on synthetic data have shown a very good behaviour also on different types of data sets.

Once the data have been collected, and objects and classes have been defined, BDT integrates feature selection and binary decision tree, according to the following steps:

- 1
Fix a level α for the decision rule. This level represents the risk we accept of a wrong attribution at each step. We use α = 0.1.

- 2
Compute the cumulative distribution in both classes, for each feature, taken one at a time, and identify the feature, and the relative threshold value, for which the statistical difference between the cumulative of the two classes is the largest. This means that the significance level of this statistical difference must be (a) lower than the level α and (b) minimum. The feature (if any) for which both (a) and (b) conditions are satisfied is the first-order feature, often called the

*root*of the pattern. On the basis of the root feature and its threshold value, each object is assigned to either one of two subsets, formed respectively by data with a value of the root feature lower/higher than the threshold. - 3
Identify the second-order features, and their thresholds, for which the statistical difference satisfies again (a) and (b) conditions. These features are found by re-analysing all the features as in (2) separately in the two subsets, and are at most two (i.e. one for each subset).

- 4
Repeat step (3) from each second-order feature, in order to identify progressively higher orders, as long as it is possible to find a feature for which the cumulative in the two classes are statistically different at a significance level lower than α. The progressive branching of the tree gives all the possible patterns. The procedure automatically terminates when no further branching is possible at the given level α.

Steps (2)–(4) are performed by means of a non-parametric Kolmogorov—Smirnov two sample statistics (Hollander and Wolf 1973). Note that the use of an *a priori* fixed level α reduces the possibility to obtain *over-fitting patterns*.

#### APPENDIX B: Fisherߣs Discriminant Analysis

This method (see e.g. Duda & Hart 1973) is based on the reduction of the *n*-dimensional space of the objects (where *n* is the number of variables describing the objects, that is, the dimension of the vectors) to an *L*—1 dimensional space (where *L* is the number of classes). In our two-class problem (*L* = 2), Fisherߣs method simply projects the objects onto a line. The basic idea, called Fisher criterion, is to project the objects onto the direction that maximizes the ratio of the dispersion between the classes to the dispersion within the classes. More rigorously, suppose we have *N* objects **x**, each represented by a vector consisting of *n* components *x _{k}*(

*k*= 1,…

*n*). Of these,

*N*

_{1}belong to class 1 and

*N*

_{2}to class 2. We linearly combine the components of

**x**, that is, the

*x*(

_{k}*k*= 1,…

*n*), in order to obtain a 1-D vector

**y**= (

*y*):

*w*

_{k}are the elements of an

*n*-dimensional vector that projects

**x**onto

**y**. In this way, we obtain

*N*objects

**y**= (

*y*), spread over the two classes.

The unknown in eq. (B1) is the projector, that is, the vector **w**. As mentioned above, we would like to choose the projection for whom the ratio of the dispersion between the classes to the dispersion within the classes is maximum. In order to do this, first we need to define some quantities.

We define **m**_{i} as the average vector of Class *i*:

**m**as the average of all the

**x**: Thus, the dispersion matrix within the class

*i*is given by: and the dispersion within all of the classes is The total dispersion matrix is given by It follows that The second addendum of the right side term of eq. (B7) is a dispersion matrix

*S*that gives an idea of the dispersion between the partial means

_{b}**m**

_{i}over the different classes and the total mean

**m**: In order to achieve the vector

*w*

^{*}that maximizes the ratio of the

*S*to the

_{b}*S*

_{w}, we need to project these matrixes onto the

**y**space and compute the

*w*

^{*}such that: Once the maximization has been carried out, Fisherߣs analysis projects the

**x**vectors onto the

**y**space, that is a line. Then, each object

**y**is assigned to the class

*i*whose mean

**m**

_{i}, projected onto the same line, is closest to

**y**.

In order to select the subset of relevant features among those available, we use a branch-and-bound technique. The basic concept in the selection of the optimal subset of features is to find, among all possible subset of the *n* features, the one leading to the lowest classification error *and* consisting of the smallest number of features. In such situation, we are confident that we are considering all of the important variables (otherwise the classification error would not be the lowest) and we are excluding the irrelevant ones (otherwise the number of features in the optimal subset would not be the smallest).

A simple, but very time consuming, way to find such an optimal subset consists of exploring the performance of the chosen NSPR algorithm on all the possible subsets of the *n* features. This becomes prohibitive as *n* increases, since we have to explore Σ^{n}_{k = 1} (^{n}_{k}) subsets. In order to avoid the application of the chosen NSPR algorithm to all the possible subsets of features, the branch-and-bound technique has been developed. This technique is applied iteratively *n* times; at each iteration *k*(*k* = 1, …, *n*), it allows the identification of the suboptimal subset consisting of *k* features, by applying the NSPR algorithm only on the ‘most promising’ subsets of *k* features. The suboptimal subset is then the one consisting of *k* features and leading to the lowest classification error.

The branch-and-bound method relies on a basic assumption, that is, it assumes that the noise introduced by irrelevant features does not deteriorate the signal given by the relevant features. In a previous study we have tested the validity of this assumption for algorithms BDT and FIS. Based on this assumption, when a certain subset of *k* features does not produce a good discrimination rule, the branch-and-bound method assumes that any other subset of *k* + *l* (*l* = 1, *n* − *k*) features containing those *k* features will not be the optimal one. In this way, a considerable portion of all the possible subsets is discarded *a priori*, thus saving computation time and effort.