Volcanic hazard assessment for the Canary Islands (Spain) using extreme value theory

The Canary Islands are an active volcanic region densely populated and visited by several millions of tourists every year. Nearly twenty eruptions have been reported through written chronicles in the last 600 yr, suggesting that the probability of a new eruption in the near future is far from zero. This shows the importance of assessing and monitoring the volcanic hazard of the region in order to reduce and manage its potential volcanic risk, and ultimately contribute to the design of appropriate preparedness plans. Hence, the probabilistic analysis of the volcanic eruption time series for the Canary Islands is an essential step for the assessment of volcanic hazard and risk in the area. Such a series describes complex processes involving different types of eruptions over different time scales. Here we propose a statistical method for calculating the probabilities of future eruptions which is most appropriate given the nature of the documented historical eruptive data. We first characterize the eruptions by their magnitudes, and then carry out a preliminary analysis of the data to establish the requirements for the statistical method. Past studies in eruptive time series used conventional statistics and treated the series as an homogeneous process. In this paper, we will use a method that accounts for the time-dependence of the series and includes rare or extreme events, in the form of few data of large eruptions, since these data require special methods of analysis. Hence, we will use a statistical method from extreme value theory. In particular, we will apply a non-homogeneous Correspondence to: R. Sobradelo (rsobradelo@ija.csic.es) Poisson process to the historical eruptive data of the Canary Islands to estimate the probability of having at least one volcanic event of a magnitude greater than one in the upcoming years. This is done in three steps: First, we analyze the historical eruptive series to assess independence and homogeneity of the process. Second, we perform a Weibull analysis of the distribution of repose time between successive eruptions. Third, we analyze the non-homogeneous Poisson process with a generalized Pareto distribution as the intensity function.


Introduction
The Canary Islands are one of the major oceanic island groups of the world and have a long magmatic history, which began at the bottom of the ocean more than 40 million years ago (Araña and Ortiz, 1991). The Canary Islands are an active volcanic region where all islands except for La Gomera show Holocene volcanic activity. Historical volcanism (last 600 yr) has been reported on the islands of La Palma (1430,1585,1646,1677,1712,1949,1971), Tenerife (1492Tenerife ( , 1704Tenerife ( , 1706Tenerife ( , 1798Tenerife ( , 1909 and Lanzarote (1730Lanzarote ( -1736Lanzarote ( , 1824. In all cases, they were eruptions of basaltic magmas characterized by emission of lava flows and construction of scoria cones. The Canary Islands are a populated ultra-peripheral Spanish region and one of the most popular touristic destinations in Europe. The presence of recurrent historical volcanism in this region is a good reason to undertake volcanic hazard assessment in order to guarantee the safety of its inhabitants Published by Copernicus Publications on behalf of the European Geosciences Union. and of its numerous visitors. Volcanic hazard is the probability of any particular area to be affected by a destructive volcanic event within a given period of time (UNESCO, 1972). As for any active volcanic region, volcanic hazard assessment on the Canary Islands requires knowing how volcanism has behaved in the past and determining the recurrence of volcanic eruptions. The first can be approached by detailed physical volcanology studies of past eruptions, in particular of those for which there exist historical chronicles (Table 1). The recurrence or eruption frequency needs to be based on historical records and precise dating of older events. Unfortunately, this is not an easy task as the reconstruction of the geological record of this volcanic region is far from accurate, lacking systematic dating of recent eruptions, and the historical records are imprecise and lack detail in some cases.
Despite these limitations, we can still analyse the volcanic hazard using the available historical data, covering the last 600 yr since the Spaniards occupied the archipelago. A set of fifteen relatively well documented eruptions form the historical volcanism of the Canary Islands (Table 2). A few other eruptions have also been reported in historical times but their age and location is imprecise and do not constitute reliable information. Most of the historical eruptions in the Canary Islands have been short lived (from few weeks to few months) basaltic, strombolian to violent strombolian eruptions, which have generated scoria cones of different sizes and lava flows of different extent (Romero, 1991). All the eruptions occurred in historical time, which goes from 1402 to present, have typically been separated a few tens of years but in some cases some have occurred in a very narrow period of time (e.g. Arafo, Fasnia, Siete Fuentes in Tenerife), or have lasted for some years (Timanfaya eruption in Lanzarote, [1730][1731][1732][1733][1734][1735][1736]. Studies of volcanic time series have been done using stochastic principles to study eruption patterns on specific volcanoes or volcanic groups (Wickman, 1976;Reyment, 1969;Klein, 1982;De la Cruz-Reyna, 1991, 1993. Bebbington and Lai (1996a) applied a Weibull renewal model to describe the patterns of New Zealand volcanoes. Other studies used transition probabilities of Markov chains (Carta et al., 1981;Aspinall et al., 2006;Bebbington, 2007), change-point detection techniques (Mulargia et al., 1987;Burt et al., 1994), Rank-order statistics (Pyle, 1998), Bayesian analysis of volcanic activity (Ho, 1990;Solow, 2001;Newhall and Hoblitt, 2002;Ho et al., 2006;Marzocchi et al., 2008;Sobradelo and Martí, 2010), non-homogeneous models (Ho, 1991;Bebbington and Lai, 1996b), a mixture of Weibull distributions (Turner et al., 2007), geostatistical hazard-estimation methods (Jaquet et al., 2000;, and a mixture of exponential distributions (Mendoza-Rosas and De la Cruz-Reyna, 2009Dzierma and Wehrmann, 2010a,b). Extreme-value methods have been applied to geological and historical eruption time series combined (Mendoza-Rosas and De la Cruz-Reyna, 2008 and historical series of large volcanic magnitudes (Coles and Sparks, 2006).
In this paper we use the historical volcanism to perform hazard assessment for the Canary Islands. Due to the limitations inherent in the available data, including its short sample time and incomplete reporting of small and intermediate magnitudes as well as uncertainties in the age, intensity and magnitude of the eruptions, we will use a method for the best estimate of the volcanic hazard based on a Non-Homogeneous Poisson process with a Generalized Pareto Distribution (GPD) as intensity function (NHGPPP) (Coles, 2001;De la Cruz-Reyna, 2008, 2010). This method has already been applied to other volcanoes for which little or incomplete data exists, like the Citlaltepetl volcano database with only six eruptions, or El Chichón volcano with 12 eruptions (Mendoza-Rosas and De la Cruz-Reyna, 2008. This is the case with our time series of volcanic records for the Canary Islands. The methodology does not require stationarity or completeness for the full eruptive series, since it depends on the number of excesses of eruptions large enough to represent the behavior of the studied volcanoes.
First, we analyze the historical eruptive series to assess independence and homogeneity of the process. Second, we perform a Weibull analysis of the distribution of repose time between successive eruptions. Third, we analyze the nonhomogeneous Poisson process with a generalized Pareto distribution as intensity function.

Geological setting
The Canary Islands are a roughly linear 500 km long chain grown on the passive margin of the African Plate within the eastern Central Atlantic Ocean (Fig. 1). The Canarian archipelago is the result of a long volcanic and tectonic activity that started at around 60 Ma ago (Robertson and Stillman, 1979;Le Bas et al., 1986;Marinoni and Pasquaré, 1994).
Several contrasting models have been proposed to explain the origin of the Canary Islands. These include a hotspot origin (Schmincke, 1982;Hoernle and Schmincke, 1993;Carracedo et al., 1998), a propagating fracture from the Atlas (Le Pichon and Fox, 1971;Anguita and Hernán, 1975), and mantle decompression melting associated with uplift of tectonic blocks (Araña and Ortiz, 1991). However, each and every one of the latter hypotheses presents some inconsistencies with the local and regional geology. A unifying model has been proposed by Anguita and Hernán (2000) who consider the existence of a residual of a fossil plume under North Africa, the Canary Islands, and western and central Europe defined through seismic tomography (Hoernle et al., 1995). Thus, volcanism is assumed to occur where an efficient fracture system allows the magma to ascend (Anguita and Hernán, 2000), i.e. the central European rift system, the volcanic provinces of the westernmost Mediterranean (Balearic and Alborán basins), Iberia, the Canary Islands and Cape Verdes (Hoernle et al., 1995). Although all islands except for La Gomera show Holocene volcanic activity, historical volcanism has been restricted to the La Palma, Lanzarote and Tenerife islands (see Fig. 1). In all cases, historical eruptive activity has been related to basic magmas ranging in intensity from strombolian to violent strombolian, originating scoria cones and lavas. In most cases, historical eruptions occurred on the active rift zones along eruptive fissures occasionally generating alignments of cones. The duration of the eruptions ranged from a few weeks to a few months, except in the case of the Timanfaya eruption in 1730 that lasted for six years. The total volumes of extruded magma range from 0.01 to >1.5 km 3 , in the case of Timanfaya. The eruption sequences that may be deduced from the successions of deposits differ from one eruption to another and reveal that eruptions did not follow a common pattern. In all cases, the resulting volcanic cones were constructed during single eruptive episodes (i.e. they must be referred to as monogenetic), commonly including several distinctive phases that do not show significant temporal separations between them. Table 2 shows the data used for this study. It includes fifteen clearly different volcanic eruptions historically documented between 1430 and 1971 for which the eruption magnitude has been computed using existing information on lava and pyroclasts volume. These data have been compiled from the original information on Table 1 and complemented with data on surface extent and volume of erupted products calculated from the geological maps at 1:25000 of IGME (Spanish Geological Survey, www.igme.es) and from a field revision of the historical eruptions that we have undertaken in this study.

Historical records of volcanic eruptions in the Canary Islands
In compiling the historical dataset of volcanic eruptions for the Canary Islands, only those eruptions with welldocumented references and clearly described eruptive processes have been considered. There are documents that make reference to other possible eruptions, but they will not be included until further documentation sources are confirmed.
The original dataset includes fifteen volcanic eruptions historically documented between 1430 and 1971 on three different islands (Lanzarote, Tenerife, La Palma). We have also considered the eruption of Montaña Cangrejo in Tenerife (Erupción de Colón in Table 1), that is supposed to have been observed by Columbus on his way to America and that has now been confirmed by Carracedo et al. (2007Carracedo et al. ( , 2011. We have considered the eruptions of Siete Fuentes (from 31 Decemebr 1704 to 5 January 1705) and Fasnia (from 5 January 1705 to 16 January 1705) in Tenerife as one unique event. The Arafo eruption (from 2 February 1705 to 27 March 1705) happened later in time and the materials have a slightly different composition than Fasnia and Siete Fuentes, suggesting that this could be a different eruption. The eruptions of Tao, Nuevo del Fuego and Tiguantón in Lanzarote are considered as one unique event. They have been listed in Romero (1991) as different episodes but they are clearly related in terms of timing, petrology and location of vents on the same eruptive fissure, so we may assume the three eruptive events were related to the same batch of magma. We consider them as part of the same eruption event if the location of the vent is not the same, but the eruptions  (Romero, 1991). Results: Dickey-Fuller = −2.5224, p-value = 0.3734 p-value > 5%, not enough statistical evidence to reject the hypothesis that the series is not stationary have the same age and the composition of the extruded products indicate they come from the same batch of magma. This is the case for Tiguantón but not the case for Siete fuentes, Fasnia and Arafo, which have a different magma composition. These criteria have been applied consistently throughout the database. In order to classify the eruptions and apply the NHGPPP, we have calculated for each case the total volume of extruded magma (DRE) based on the volume of exposed materials (lavas and pyroclastic deposits), so our total volumes are minimum estimates (Table 2). Although rapid erosion of tephra and uncertain lava flow thicknesses may cause problems in making accurate calculations, order of magnitude determinations still provide a useful comparison between eruptions. We have computed the volume estimates calculating separately the volumes of tephra and lavas for each eruption. Tephra and lava volumes have been calculated with the help of a DEM at a resolution of 5 m, the digital geological maps of IGME at 1:25000, and checking extension and thicknesses variations of the deposits and lavas in the field.
Volcanic eruptions are natural phenomena where the frequency of the events decreases as their size or magnitude increases. The fact that the magnitude distribution is irregular is not necessarily an indication of incompleteness in the catalog. When we are dealing with a historical catalog, it is very difficult for a high magnitude eruption to go unnoticed. Compared with other volcanoes, a catalog of 15 eruptions in 600 yr seems consistent. There are no records of any more high magnitude eruptions in historical times other than the high magnitude eruption of Timanfaya. For this reason, we assume that the catalog for high magnitude eruptions in historical times is complete. On the other hand, the historical records for the oldest low magnitude eruptions are not as clear and accurate as for the most recent low magnitude eruptions.
To deal with the difficulties derived from the possible lack of catalog completeness for the Canary Islands, we measure the size of each eruption with its magnitude based on the logarithmic scale for magnitude (Pyle, 2000): We have calculated the total magnitude of each eruption (in Kg) assuming a density for the basaltic magma of 2850 kg m −3 (Table 2). To account for the possible missing data catalog, for which historical records are inaccurate or unavailable, observed occurrence rates for magnitudes 4 and 6 were used to extrapolate unobserved records using the best fit to the class magnitude values of eruptions estimated with the power law described by Newhall and Self (1982), which says that the eruption occurrence rate λ M (number of eruptions per unit time) of each class magnitude M is related to the eruption size M as: where a and b are constants that describe the power law decay of occurrences with increasing magnitude over a given time interval.

The method: extreme value theory (EVT)
Volcanic eruption datasets are usually small and the eruptive recurrence is usually very long, and as it happens with other natural phenomena like earthquakes, tsunamis etc., the higher the magnitude the longer the time interval in between events. To face this problem of working with small datasets, and to be able to obtain a mathematical quantification of the volcanic hazard as accurate as possible, we look for methods that allow us to work with databases which are small and sometimes incomplete (Coles, 2001;Davison and Smith, 1990;Beguería, 2005). These methods are part of a branch of statistics called extreme value theory, in which as the name implies, extreme values are atypical and rare events located at the tail of the distribution.
Just as normal distribution proves to be the important limiting distribution for sample sums or averages, as is made explicit in the central limit theorem, another family of distributions proves important in the study of the limiting behaviour of sample extrema. This is the family of extreme value distributions. Extreme value theory and the central limit theory are derived in a similar manner. Both consider the limiting distributions of independent identically distributed (iid) random variables under an affine transformation. In the absence of empirical or physical evidence for assigning an extreme level to a process, an asymptotic argument is used to generate extreme value models. But extreme values are scarce, making it necessary to estimate levels that are much higher than those that have already been observed. In fact, the goal of an extreme value analysis is to quantify the statistical behavior of processes at unusually high levels. In particular, extreme value analysis requires an estimation of the probability of events that are more extreme than any that have ever been observed. This implies an extrapolation from observed levels to unobserved levels. Extreme value theory provides a family of models to make such extrapolation. In fact there are no more serious competitor models than those provided by extreme value theory (Coles, 2001).
There are different extreme value theory methods (Coles, 2001). Depending on how we define our extreme values, we select the method. In our case, it is more convenient to define our values as peaks or exceedances over a threshold, and so we use the Exceedances over a threshold (EOT) method to sample the original data, i.e. X i > u for some value of i. This method is based on a limiting function called GPD, as opposed to the Annual Maximum method with is based on the Gumbel distribution as the limiting function.
The family of GPDs describes the behavior of individual extreme events. It considers observations from a collection of iid random variables where we keep those that exceed a fixed threshold u. As we increase the threshold, the two-parameter GPD family represents the limiting behavior of this new collection of random variables. This makes the family of GPDs a suitable choice for modeling extreme events.
The EOT method includes all the values of the variable that exceed an a-priori established threshold, u, fixed according to the model needs, providing a physically based definition of what must be considered an extreme event. The choice of the threshold value has a strong subjective component. This random variable is defined by the transformed random variable where Y is the excess over the threshold u.
The parameter that will be used as random variable to estimate the probability of occurrence of a future eruption, and thus the volcanic hazard, will be the time interval T between eruptions, also called repose period, together with the magnitude M.
The generalized Pareto distribution can be fitted to data on excesses of high thresholds by a variety of methods including the maximum likelihood method (ML) and the method of probability weighted moments (PWM). We use the Davison and Smith (1990) The NHGPPP is also appropriate for linking historical and geological records, should they become available in the future. So this method sets the base for future analyses and updates should geological records were found.
As a first step before model fitting is undertaken, a number of exploratory graphical methods provide useful preliminary information about the data and in particular, their tail. We explain these methods in the next section. To apply the NHGPPP for volcanic hazard assessment, we first need to examine the data to assess independence between successive events and homogeneity of the process. We will use a serial correlation scatterplot to test for independence and to test for homogeneity, we first assess the stationarity by using the autocorrelation function (ACF) and the Dickey-Fuller unit root test. These tests should be done on a portion of the time series in which no significant eruption data are missing, which in most cases is the historical eruption dataset of intermediateto-high eruption magnitudes.
After independence and homogeneity have been assessed, we do a Weibull analysis of the repose periods between eruptions to quantitatively describe the stationarity of the series through the distribution shape parameter. The further from 1 the shape parameter is, the more evidence that the process is not stationary.
After the data has been analyzed we apply the NHGPPP to estimate the volcanic hazard. The method is applied to an independent, non-overlapping series of events occurring in a space B with an intensity density λ(x i ), where x i are the Bdomain variables where the process develops. In our case x i are the coordinates T (time) and eruption magnitude M of a two-dimensional space.

Statistical analysis of the Canary Islands historic volcanic data using EVT
Assuming that past history of a volcano should reflect at least some relevant features of its expected future behavior, we look at the time series of historical volcanic eruptions on the Canary Islands. The time series dataset has fifteen volcanic events historically documented since 1430 (Romero, 1991).

Independence and stationarity assessment
We want to test the eruptive time series for independence between successive events. We look at the independence of repose intervals between eruptions by means of a serial correlation scatterplot, where the duration of each interval T i +1 between two successive eruptions is plotted against the previous repose interval durations T i . For the particular case of the Canary Islands, the repose interval between eruptions is defined as the elapsed time between the start of one eruption and the start of the next. The diagram in Fig. 2 shows a large dispersion of points and the correlation coefficient between consecutive repose times is 0.3062, indicating a very low serial correlation. We do not have enough evidence to say that consecutive repose intervals are time-dependent. If new geological data arrive in the future we do not rule out the possibility of a new outcome for the time-dependence analysis, but for the time being, based on the available data, we assume independence of repose times based on the above mentioned tests.
Next, we look at the stationarity of the process. A time series is stationary if its underlying statistical structure does not evolve with time. The correlogram is a simple diagram which can help diagnose non-stationarity. If a series is nonstationary then the theoretical autocorrelations will be nearly 1 for all lags k. Thus, if the estimated correlogram fails to die down (or dies down very slowly), the series is non-stationary. The theoretical correlogram is a plot of the theoretical autocorrelations between consecutive repose periods of lag k, corr(x t ,x t−k ), against k. Figure 3 shows the autocorrelation function (ACF) of the Canary Islands' time series.
The argumentation of the non-stationarity based in the shape of the ACF is arguable since the ACF is sensible to seasonal variations, which at the same time could correspond to a stationary process. For this reason, to assess stationarity, we complement the visual ACF analysis in Fig. 3 with the Dickey-Fuller unit root test.
The Dickey-Fuller unit root test was proposed by Dickey and Fuller (1979). In its most basic form, the test compares the null hypothesis H 0 : x t = x t−1 + t , i.e. that the series is a random walk without a drift, against the alternative hypothesis H 1 : x t = c + αx t−1 + t where c and α are constants with | α |< 1. According to H 1 , the process is a stationary AR(1) with mean µ = c/(1 − α). To see this, note that, under H 1 we can write x t = µ(1 − α) + αx t−1 + t , so that Table 3 shows the SAS output for the Dickey-Fuller test. The test statistic has a value of −2.5224, and is associated with a p-value of 0.3734, indicating that there is not enough statistical evidence to reject the null hypothesis that the series is not stationary. This result is consistent with the visual analysis of the ACF, where the series fails to die down. In this preliminary analysis of the time series no significant correlation was found, thus we can assume independence of consecutive repose periods. Additionally, we found no evidence to assume that the series is stationary, so based on the ACF and the Dickey-Fuller test for stationarity, we can say that the Canary Islands' volcanic eruptions time series is not stationary.

Distribution of the repose periods: Weibull versus Exponential
We look at the Weibull distribution to analyze the characteristics of consecutive repose periods and quantitatively describe the stationarity characteristic of the time series through the distribution shape parameter. The Weibull distribution has been widely applied in statistical quality control, earthquake hazard assessment, and many other applications. It has also been used to model volcanic eruption sequences (Bebbington and Lai, 1996b). The 2-parameter cumulative Weibull distribution and survival functions are given by respectively, where α is a scale parameter, and k is a shape parameter.
The shape parameter reflects the stationary or nonstationary character of the time series. In the present paper, we obtain the distribution parameters using a fairly simple graphical method (Bebbington and Lai, 1996a). The probability of having a repose period of duration greater than t has been obtained from the survival function 1 − F (t).
The resulting Weibull distribution parameters are 1.63 for the shape parameter and 4.37 for the scale parameter. Figure 4 shows the comparison between Exponential and Weibull distributions. We see that the Weibull survival function provides a better fit to the repose periods than the exponential function, because the shape parameter accounts for the non-stationarity of the time series. Additionally, the shape parameter value being far from 1 confirms once more the non-stationarity of the process. A shape parameter of 1 would correspond to an exponential, that models very well stationary data, which is not the case here.

Volcanic hazard assessment for the Canary Islands
We then estimate the volcanic hazard for the Canary Islands using the NHGPPP. A Poisson process is a collection  When the rate parameter, or intensity, of the process is not constant, the Poisson process is said to be non-homogeneous, and the generalized rate function is given by λ(t). As seen in a preliminary analysis of the data, the Canary Islands time series is non-stationary, and we will model the volcanic hazard with a non-homogeneous Poisson process (NHPP). Since we use the EOT method to sample the original data, we can use the GPD to model the intensity of the NHPP. Hence, we will be using a NHGPPP to estimate the volcanic hazard for the Canary Islands. (See Mendoza-Rosas and De la Cruz-Reyna (2008) for further details on this methodology).
To calculate the probabilities of occurrence of an eruption in the different magnitude classes, we use the number of excesses inferred from the eruption occurrence rate of each class magnitude, this is, the events above a certain threshold u (X i − u, for some i).
For the particular case of volcanic eruptions, the magnitude of the eruptions and the time of their occurrence are viewed as points in a two-dimensional space, which formally is the realization of a point process (Cox and Isham, 1980). The intensity measure (B) of this two-dimensional Poisson process on the space B = [t 1 ,t 2 ] * (u,∞) with [t 1 ,t 2 ] ⊂ [0,1] is given by where β, and θ are the parameters of the GPD. The GPD is described by a shape parameter β, a scale parameter θ, and a location parameter u (threshold), and has the following cumulative distribution function: Table 4. Probability of at least one event of Magnitude M > x in the next t years in the Canary Islands estimated with a NHGPPP. (Pr(X = 0) and Pr(X ≥ 1) are the probability of having no eruption and the probability of having at least one eruption, respectively, of a certain magnitude in a particular time interval; NHM and HM are the probability estimated with the NHGPPP and the Homogeneous Poisson, respectively, for magnitude M;λ is the estimated parameter rate for the NHGPPP andσ is the estimated standard deviation for the Pr(X ≥ 1) computed with the NHGPPP, based on the Delta method.) Magnitude > 1 NHM>1 HM>1

Years
Pr(X = 0)λ Pr(X ≥ 1)σ Pr(X ≥ 1) Another related property of the GPD refers to the mean excess: if Y = X −u is a GP-distributed variable, then the mean excess over threshold u is: for β > −1, u > 0 and θ − uβ > 0 The sample mean excess plot is given by: where N u is the number of excess x i over a threshold u. Davison and Smith (1990) introduced a diagnostic plot to decide how well the model fits the data. The mathematical basis for this method is Eq. (3), where the key feature is that if Y is GPD, then the mean excess over a threshold u, for any u > 0, is a linear function of u (Coles, 2001;Lin, 2003;Beguería, 2005). In Fig. 5 we plot the mean of the excesses, obtained with Eq. (4), vs their thresholds. The x-axis is the threshold and the y-axis is the sample mean of all excesses over that threshold. As we can see, the mean excess follows a nearly straight line, with a R 2 of 0.8415, suggesting a good fit. A regression line of mean of exceedances over the threshold has been added to confirm that the series follows the GPD. Since we are working with effusive eruptions only, we assumed an upper bound of six for the estimation of the magnitudes, and mapped the probabilities to the [1,6] magnitude interval. Hence, according to Davison and Smith (1990), the preceding results indicate that the NHGPPP is satisfactory and appropriate to model the Canary Islands eruptive time series.
The Pareto generalized parameters for the process, derived from the regression parameters on Fig. 5 and Eq. (3) are 0.104 for shape and 1.711 for the scale. Using Eq. (2) we estimate the intensityλ of the NHGPPP and obtain the probability estimations of at least one eruption of a certain magnitude in a given time interval. Since the approach of exceedances assumes that the scale measuring the phenomena is open, we mapped the probabilities of eruptions to the magnitude interval 1 to 6. Table 4 and Fig. 6 show the probability of having at least one eruption Pr(X ≥ 1) computed as one minus the probability of having no eruption 1 − Pr(X = 0) of a certain magnitude M in a given time window, estimated using the NHGPPP with intensity rateλ. To measure the volatility of the estimated probabilities Prob(X ≥ 1) we compute the standard deviationσ of the estimator. We use the Delta method to determine its asymptotic distribution, and we getσ = 1 n e −2λλ . Additionally, in Fig. 6, we compare the results obtained from the NHGPPP with the volcanic hazard estimates obtained from direct application of the homogeneous Poisson distributions for the same eruption series. We see that the probabilities of occurrence of eruptions in the lower magnitude scales, calculated with the non homogeneous method proposed here, differ very little from the standard Poisson method. The probabilities of occurrence of eruptions exceeding moderate magnitudes are significantly higher with the NHGPPP. This difference is due to the additional information that the GPD adds to the NHGPPP when Based on the existing historical data, the probability of an event in the Canary Islands increases more rapidly in the first 20 yr, with a 99.84 % chance of an event of a magnitude greater than one in the next 20 yr and leveling out after that at 99.99 %. There is a probability of 27.58 % of an event of a magnitude between 1 and 6 in the next 12 months and 3.71 % of an event of a magnitude between 4 and 6 for the same period. There is ongoing work to assess the data in the volcanic eruptions catalog for the Canary Islands more accurately. In this respect, these probability results may vary should new geological records become available.

Discussion and conclusions
The Canary Islands are an active volcanic region densely populated and visited by several million tourists every year. Nearly twenty eruptions have been reported by written chronicles in the last 600 yr. This gives an average of an eruption every 25-30 yr, which suggests that the probability of having a new eruption in the near future is not so low. Under these circumstances and in order to reduce the potential volcanic risk of this region, it is highly recommendable to undertake hazard assessment, and determine the eruption recurrence for the area.
Recent volcanism in the Canary Islands is not well known and is poorly constrained in terms of age of the eruptions. For the island of Tenerife alone, Carracedo et al. (2007Carracedo et al. ( , 2011 have conducted a systematic geochronological study for Teide and the volcanism associated with the rift zones, but this study is still far from being complete. Therefore, the data catalog to be used for statistical and probabilistic assessment of the Canary Islands to establish the eruption recurrence is formed uniquely by historical records. The model can be easily updated in the future should new volcanic records be dated.
As in any data analysis, we should be aware of various layers of uncertainty, perhaps magnified in an extreme value analysis. On one level, there is the parameter uncertainty: even if we had abundant, good quality data to work with and a good model, our parameter estimates are still subject to a standard error. Model uncertainty is also present -we may have good data but a poor model. Using extreme value methods we are at least working with a good class of models, but they are applicable over high thresholds and we must decide where to set the threshold. If we set the threshold too high we have few data and we introduce more parameter uncertainty. If we set the threshold too low we lose our theoretical justification for the model. But even more serious than parameter and model uncertainty is the data uncertainty. It is never possible to have enough data in an extreme value analysis. In Table 4 the expected value of the random variableλ provides a good estimation of the aleatoric uncertainty due to the complexity of the process, and the standard deviationσ provides a good estimation of the epistemic uncertainty, due to our limited knowledge of the process.
Extreme value methods do not predict the future with certainty, but they do offer good models for explaining the extreme events we have seen in the past (McNeil, 1997). Even with a good tail estimate, we cannot be sure that the future does not hold some unexpected catastrophic volcanic eruption. The extreme value method used in this paper to assess the volcanic hazard for the Canary Islands does not predict the future with certainty, but it is a model based on rigorous mathematical theory concerning the behaviour of extrema. Based on past experience (Mendoza-Rosas and De la Cruz-Reyna, 2008), the GPD is a good approximation in the tail for our volcanic data, and the probability results yielded by the extreme value method used here to assess the volcanic hazard for the Canary Islands should not be neglected. It may well be that, by trial and error, some other distribution can be found which fits the data even better in the tail, but such a distribution would be an arbitrary choice, and we would have less confidence in extrapolating it beyond the data.
The probability results obtained are very high. This is partly due to the fact that the area of study is the quasi linear 500 km long chain grown on the passive margin of the African Plate containing the actual Canarian archipelago, signifying that there are several possible vent locations for an eruption. Also, we must consider the fact that we are measuring magnitudes (total erupted volumes) and not VEI (Volcanic Explosivity Index) Newhall and Self (1982). The VEI is a combination of volume of products, eruption cloud height and qualitative observations. It is mainly applied to explosive eruptions and is not appropriate for eruptions which are mainly effusive. This is the reason why we used magnitude instead of VEI and limited the study to levels of magnitude up to six. However, the eruption magnitude, measured as the total erupted volume, only takes into account one of the three measures of the VEI. Hence, the probability estimates for volume alone are expected to be higher than the estimates for volume of products, eruption cloud height and qualitative observations combined. With this in mind, given the current data, it is not surprising to observe a probability of 27.58 % of having a volcanic event of magnitude greater than 1 in the next year in the Canary Islands, most likely in any of the three for which historical data exist (Lanzarote, Tenerife, La Palma) but without excluding the other four islands. Even if there are no historical records documented for all the islands, we cannot rule out the probability of an event forming there since they are part of the same archipelago and there are traces of previous volcanic eruptions. We do not have enough data to do an individual hazard analysis for each island alone.
It is important to highlight the fact that the Canary Islands have a probability greater than zero of undergoing a new volcanic event in the upcoming years. Hence, these results should be taken into account in the assessment of volcanic risk and in the design of prevention and response measures, particularly of major eruptions to which larger areas may be one hundred per cent vulnerable.
The results obtained only apply to the probabilities of having a basaltic eruption in the near future, as historical volcanism has been always represented by this kind of eruption. However, the existence of several eruptions of phonolitic magmas from Teide in Holocene times on the island of Tenerife, the last one having occurred only 1000 yr ago (Carracedo et al., 2007), reminds us that hazard assessment should also consider phonolitic eruptions. Despite their being concentrated during the Holocene exclusively on Tenerife, these eruptions may generate hazards that could have a much greater impact than basaltic eruptions; thus, their potential effects should not be neglected in a more complete volcanic hazard assessment for the Canary region.