 Research
 Open Access
 Published:
Phreatic eruptions at crater lakes: occurrence statistics and probabilistic hazard forecast
Journal of Applied Volcanologyvolume 6, Article number: 4 (2017)
Abstract
Phreatic eruptions, although posing a serious threat to people in crater proximity, are often underestimated and have been comparatively understudied. The detailed eruption catalogue for Ruapehu Volcano (New Zealand) provides an exceptional opportunity to study the statistics of recurring phreatic explosions at a crater lake volcano. We performed a statistical analysis on this phreatic eruption database, which suggests that phreatic events at Ruapehu do not follow a Poisson process. Instead they tend to cluster, which is possibly linked to an increased heat flow during periods of a more shallowseated magma column. Larger explosions are more likely to follow shortly after smaller events, as opposed to longer periods of quiescence. The absolute probability for a phreatic explosion to occur at Ruapehu within the next month is about 10%, when averaging over the last 70 years of recording. However, the frequency of phreatic explosions is significantly higher than the background level in years prior to magmatic episodes.
Combining clast ejection simulations with a Bayesian event tree tool (PyBetVH) we perform a probabilistic assessment of the hazard due to ballistic ejecta in the summit area of Ruapehu, which is frequently visited by hikers. Resulting hazard maps show that the absolute probability for the summit to be affected by ballistics within the next month is up to 6%. The hazard is especially high on the northern lakeshore, where there is a mountain refuge.
Our results contribute to the local hazard assessment as well as the general perception of hazards due to steamdriven explosions.
Introduction and background
Not all explosive activity at volcanoes or hydrothermal areas is directly linked to the ascent and eruption of fresh magma. There are also events caused by the explosive expansion of water, steam and/or other gases; these are described by many terms, including phreatic eruption, hydrothermal explosion, steamblast or nonjuvenile eruption, and the definitions are not always consistent. In Barberi et al. (1992), phreatic explosions are defined as explosions of a confined pocket of steam and gas with no direct involvement of magma. Other authors include eruptions where minor amounts of juvenile material are ejected, as long as gas is the main driver of the eruption (e.g. Kilgour et al. 2010). Browne and Lawless (2001) distinguish between hydrothermal, phreatic and magmatichydrothermal eruptions depending on the main energy source, while Mastin (1995) uses the underlying process to define three types of nonjuvenile eruptions: gas eruptions, mixing eruptions and boilingpoint eruptions. The latter are the most common and are caused by the depressurisation of water near its boiling point. In this study, we use the term "phreatic eruption" for all explosive activity caused by the explosive expansion of water, steam and/or gas  independent of the fluid source.
These explosions occur at many different volcanoes; Soufrière de Guadeloupe (Le Guern et al. 1980), Meakandake (Ogiso and Yomogida 2012), Usu (Miura et al. 2012), Bandai (Glicken and Nakamura 1988), Karkar (Barberi et al. 1992), Nyos (Neri et al. 1999), Kilauea (Dvorak 1992), Mayon, Pinatubo, Mount St Helens (Newhall et al. 2001), Ruapehu (Kilgour et al. 2010) and Poás (Fischer et al. 2015) are just a few examples. However, as their magnitude is generally smaller than that of magmatic events, phreatic eruptions are comparatively understudied. Due to minor volumes of erupted material, the deposits are typically not well preserved in the geological record and phreatic eruptions are also likely to be overlooked in observational reports, especially when they are precursors to a magmatic eruption. But even on their own, phreatic eruptions pose some significant hazards: ballistic blocks, toxic gases, base surges, lahars, directional blasts and debris avalanches have been reported (Barberi et al. 1992). The particular danger lies in the difficulty to recognise the commonly very shortlived and highly localised precursors; in some cases precursors are not identifiable (e.g. Barberi et al. 1992; Kato et al. 2015; Maeda et al. 2015). As a consequence, the forecasting of phreatic eruptions is extremely challenging.
Due to their comparatively small impact radius, direct hazards from phreatic explosions are concentrated close to the crater. However, recent events at Mayon in 2013 with 5 casualties (Venzke 2013) and Ontake in 2014 with at least 57 casualties (Kato et al. 2015) demonstrate that the impact of localised phreatic eruptions is not to be underestimated by the authorities and that the risk at frequently visited, apparently quiescent, volcanoes might be higher than anticipated.
There are very few detailed, quantitative studies of hazards related to phreatic activity. An exception is the generic investigation by Neri et al. (1999) that adresses the hazards posed by phreaticrelated surges, blasts and toxic gases using simple numerical simulations. A second example is the case study presented by Fitzgerald et al. (2014), who combined detailed field data analysis from a phreatic eruption at the Upper Te Maari Crater in 2012 with the simulation of 3D ballistic trajectories to assess ballistic hazards at this volcano.
This study focuses on phreatic eruptions at Ruapehu, New Zealand (Fig. 1), an active andesitic stratovolcano with frequent phreatic and phreatomagmatic eruptions, as well as two larger historic magmatic episodes (Kilgour et al. 2014). The crater hosts a hot acidic crater lake covering an active hydrothermal system. Phreatic explosions are caused by the volatile and heat release from shallowseated magma into this system, and/or a pressure build up in the hydrothermal reservoir by sealing mineralisation and the rupture of the resulting cap (e.g. Christenson and Wood 1993; Christenson et al. 2010). The latter occurred for example in September 2007, producing a surge, lahars and ballistics and seriously injured a hiker that was present on the summit at the time (Kilgour et al. 2010). Seismic anomalies have been reported prior to previous phreatic eruptions at Ruapehu (e.g. Dibble 1969; Nairn et al. 1979; Jolly et al. 2010; Mordret et al. 2010; Carniel et al. 2013), but have only retrospectively been identified as precursors and therefore could not be utilised as timely warning systems. For the 2007 eruption, using hindsight analysis Jolly et al. (2010) found seismic precursors starting 10 min before the eruption  therefore, even if a monitoring system had identified these signals, little time remained to evacuate potential visitors from the summit.
The very detailed catalogue of Ruapehu’s activity (Scott 2013) offers an exceptional opportunity to study the statistics of recurring phreatic events and thereby support the assessment of hazards that might affect people in the summit area, which is frequently visited by hikers. We present an analysis of the recurrence statistics of past events at Ruapehu, followed by a probabilistic hazard assessment. The methodology and results are presented individually, and then discussed together.
Recurrence statistics
Methods
Scott (2013) compiled the most detailed record of the activity at Ruapehu since 1830 based on scientific and public sources. The activity is classified into 6 eruption scales from 0 to 5 (Table 1) based on their impact range. Scott (2013) identified several periods of higher than normal activity. Two magmatic episodes, September 1944 to July 1946 and June 1995 to November 1997, make up about 65% of the observations and Scott (2013) suggests removing these for the assessment of the volcano’s background behaviour. The completeness of the data set was assessed by Scott (2013) through comparing the cumulative frequencies of eruptive scales, which suggests that eruptions of Scale 3 and smaller are underreported. Due to improved monitoring, the data set is considered more complete from 1940 onwards (Scott 2013), although events confined to the lake (i.e. Scale 01) can still be missed due to the lack of continuous visual observation.
For our purposes we extracted all reports of phreatic eruptions from the dataset from 1940 onwards, hence neglecting all magmatic events. However, due to the varying definitions and terms used in literature and in the dataset, it is often unclear whether an event was purely phreatic or phreatomagmatic. Table 2 lists all terms that we interpreted as phreatic eruptions, as well as some ambiguous terms that could have been phreatic but have been interpreted as magmatic and are therefore not included in our analysis.
After identifying all phreatic events from the original catalogue, only Scale 1 to 4 remain in our database. The original database has a daily resolution; for our purposes we use a monthly resolution, i.e. a month is counted as a success, if there is at least one day with a reported phreatic explosion. We use monthbased bins because there is often a large uncertainty in the database concerning how many events occurred in severaldayslong periods. Hence, all subsequently calculated probabilities should be understood within this timeframe, i.e. they represent the probability for at least one event occurring within the next month. A basic statistical analysis of our database, like tests for a Poissondistributed behaviour or a scaletime relationship, was realised with simple MATLAB®; routines.
Results part one: recurrence statistics for phreatic eruptions at Ruapehu
The cumulative number of phreatic eruptions with time (Fig. 2) shows periods of increased activity and periods of relative quiescence. Some periods of phreatic quiescence coincide with magmatic eruptions, but there are also real periods of quiescence with neither phreatic nor magmatic activity (e.g., August 1959 to May 1964). The known magmatic phases are preceded by an increase in phreatic activity. However, the opposite does not hold true, i.e. not all increases in phreatic activity are followed by a magmatic phase.
In total, 90 months with eruptions were counted. If we average this over the whole time span from 1940 until January 2016, the overall relative occurrence of months with eruptions is 10%. However, there are periods with a higher eruption frequency than others. Periods 2, 3 and 4 indicated in Fig. 2 precede magmatic phases and show a significantly higher relative occurrence of months with eruptions of 31, 32 and 17%, respectively. For comparison, Period 1  between July 1946 and March 1966  has an eruption frequency of only 6%. In the time since the last magmatic phase in 1997, the average eruption occurrence has been even lower: 3%.
The average of the intervals between two phreatic eruptions, i.e. the mean of intereruptive time (“IET”) μ _{IET}, is 9 months, with a standard deviation σ of 14 months. Therefore, the coefficient of variation (Cox and Lewis 1966) is \(\text {CV}=\frac {\sigma }{\mu _{\textrm {IET}}}=1.6\). This value is an indicator of the statistical nature of processes: A Poisson process must have a CV value close to one, while for regularly occurring events CV is <<1 and clustering processes typically show a value >1. Therefore, it seems that phreatic eruptions at Ruapehu cluster. A χ ^{2}test confirms that the phreatic events do not follow a Poisson distribution: the Poissonhypothesis is rejected at the 0.05 significance level (pvalue of 0.0135).
Figure 3 a shows the cumulative number of phreatic eruptions individually for the different eruption scales and illustrates the expected much higher frequency of smallerscaled events. The intervals between two eruptions of Scale 3 or 4, respectively, can be much longer than those between smaller scaled eruptions (Fig 3 b). However, no statistical significant relationship between IET and eruption size can be proven (see Appendix).
Probabilistic hazard assessment using PyBetVH
Methods
Following the statistical analysis, we performed a probabilistic assessment of the hazard due to ballistic ejecta in the summit area of Ruapehu using PyBetVH (Tonini et al. 2015). Probabilistic volcanic hazard assessment (PVHA) aims to include the intrinsic uncertainty in forecasting volcanic activity. PyBetVH is the software implementation, coded in Python programming language, for the BET_VH model (Marzocchi et al. 2010; Tonini et al. 2016), a Bayesian event tree model for longterm PVHA. Prior knowledge, e.g. from model results or a priori beliefs, are combined with frequencies from the catalogue and the procedure accounts for both the aleatoric and epistemic uncertainty. The final products are conditional and absolute probability maps for volcanic hazardous phenomena. In the following, we describe the creation of the input files for PyBetVH, following Tonini et al. (2015).
Nodes 13
Usually, these nodes represent the probability that a volcano enters unrest (node 1), that the unrest is of magmatic origin (node 2) and that the unrest leads to an eruption (node 3), respectively. Since phreatic eruptions often occur independently of any magmatic unrest (Rouwet et al. 2014), we combine these nodes such that node 1 simply represents the probability of a phreatic eruption and ignore nodes 2 and 3 by assigning them a probability of 1. The assigned prior probability θ for an eruption is 0.5, with λ=1. λ is the “equivalent number of data”, which expresses the confidence in the prior probability value: a low value indicates that only very few data would suffice to change our guess of the prior probability. Our chosen combination of θ and λ for nodes 1  3 expresses the maximum ignorance probability distribution, i.e. we do not know anything about the real probability distribution. This prior probability is combined with past data: our phreatic eruption catalogue spans a time period of 789 months (counting only those that did not begin already in a state of phreatic eruption), and we count a total of 90 months with phreatic explosions.
Node 4
This node represents the probability for an eruption originating from a certain vent position. PyBetVH has two categories for the geometry of possible vent locations: a “volcanic field” corresponds to a rectangular grid of possible vents, while a “volcanic cone” is characterised by a central vent position over the summit crater area and 4 lateral sectors. Ruapehu falls in the latter category, with the crater lake as the central vent. Since it seems highly unlikely that the eruption does not originate from the crater lake, we assign a prior probability of 0.99 to the central and 0.0025 to each of the four radial sector positions, respectively. These prior beliefs are combined with past data: all of the recorded eruptions originated from the crater lake.
Node 5
Node 5 represents the probability for an eruption to be of a certain scale. To obtain prior probabilities, we make use of the general assumption that the scalefrequency relation of eruptions follows a power law. From the global magmatic eruption catalogue, Marzocchi et al. (2004) determined the relation log(N)=5.811−0.7794∗VEI, with N being the number of eruptions of the respective VEI in 1000 years. As a best guess, we assume that the scalefrequency relation of phreatic eruptions at Ruapehu also follows this law, leading to a probability of 0.834 for Scale 1, 0.139 for Scale 2, 0.023 for Scale 3 and 0.004 for Scale 4. It should be noted that the relation from Marzocchi et al. (2004) is based on a database of eruptions with VEI ≥2 and therefore might not be valid for phreatic eruptions at Ruapehu. Consequently, we assigned this prior probability a very low confidence value, i.e. λ=1. In combination with the large number of recorded eruptions, this means that the posterior probability is mainly determined by the data from the catalogue. As input for past data, we extracted 42 Scale 1, 28 Scale 2, 13 Scale 3 and 7 Scale 4 eruptions from the catalogue  corresponding to relative frequencies of 0.467, 0.311, 0.144 and 0.078, respectively.
Node 6
This node reflects the probability for the occurrence of specific hazards. Here, we consider the hazard of ballistic ejecta, as these pose a serious threat to people on the summit. By definition, all eruptions of Scale 2 and larger deposit material outside of the lake, while Scale 1 eruption products are confined to the lake. Hence we assign a probability of 1 as the prior probability for Scales 2 to 4, and include all recorded eruptions of these scales as past data for the occurrence of ballistics. The prior probability is set to 0 accordingly for Scale 1, and all past eruptions of Scale 1 are included as a recorded lack of ballistics. λ is set to 1.
PyBetVH calculates the probability of a hazard to overcome intensity thresholds, which are defined in this node. As intensity measure we choose the density of ballistic impacts F (following Fitzgerald et al. (2014)), and set the thresholds 0, 0.001, 0.005 and 0.01 clasts/m^{2}. We assume an average hazardous area of 10 m^{2} per clast (following estimations in Fitzgerald et al. (2014)) that is multiplied by the impact density to calculate the average number of clasts that affect a person on the summit.
The definition of areas on the map, where the hazard is to be assessed, is also included in this node by the definition of grid points; we use a radial grid (see Fig. 4). First, we define radial areas (rings) A _{ n } surrounding the crater lake by their minimum radial distance r _{ n } from the lake centre (250, 400, 550, 700, 850, 1000, 1150, 1300, 1450, 1600, 1750, 1900 and 2050m distance, the upper bound is 3950 m). For each of these rings, 24 equally spaced points are then defined on its central radius as the grid points for BET (see Fig. 4).
Nodes 7&8
In this combined node, the prior probability for a hazard to overcome an intensity threshold in a certain area of interest is defined. To obtain the probabilities for ballistic impacts on Ruapehu’s summit to overcome impact density thresholds, we ran 60,000 simulations of ballistic trajectories (see section “Assessment of probabilities of ballistic impacts on Ruapehu’s summit”) with random input parameters, providing a statistical distribution for distances covered by ejected ballistic clasts. In combination with the number of ejected clasts that make it out of the crater (m), we can calculate the expected impact density F in different areas on the summit (see section “Assessment of probabilities of ballistic impacts on Ruapehu’s summit”) for a given eruption scale. Since m is poorly constrained, we use 10,000 randomly sampled m to obtain an empirical probability density function for F and use this to define the probability that the impact density in an area of interest overcomes a certain threshold (see section “Assessment of probabilities of ballistic impacts on Ruapehu’s summit”). Again, λ is 1. This prior knowledge is combined with data from the 2007 phreatic eruption (Kilgour et al. 2010), where ballistics were mapped: they occurred up to 2 km away from the crater, with an impact density of F=0.014 clasts/m^{2} and are constrained to the northern side of the lake within a pitch circle with an opening angle of about 40°.
Assessment of probabilities of ballistic impacts on Ruapehu’s summit
To assess the hazard due to ballistic ejecta in the summit area that can be produced by phreatic eruptions of Scale 2 to 4, we developed a MATLAB®; procedure, which performs the following steps:

1.
Generate random input parameters for ballistic flightpaths (based on eruption scale)

2.
Calculate ballistic trajectories and infer the impact distances by comparing the trajectory with Ruapehu’s topography

3.
Derive the probability for one ejected clast to reach a certain area on the summit

4.
Derive a probability density function for expected impact densities, given an eruption of a certain scale, for the different areas.
The steps are explained in more detail below.
Step 1: Random input parameters
Ballistic trajectories are mainly determined by eruption parameters, i.e. the ejection velocity and angle, and physical properties of the projectiles such as their size, shape and density. Ballistics from a phreatic eruption at Ruapehu are described for the 2007 eruption: they were predominantly angular to subangular fragments of hyalocrystalline to hypocrystalline pyroxene andesite, with a lesser amount of andesite breccia cemented by a medium grey mud (Christenson et al. 2010; Kilgour et al. 2010). Projectiles of a variety of volcanolacustrine lithologies made up a subordinate group. The clasts showed varying degrees of hydrothermal alteration and vesicularity; their densities ranged from 1700 kg/m^{3} to 2700 kg/m^{3} (Kilgour et al. 2010). In our simulations, we assume cubeshaped projectiles and fix clast density at the mean value determined for ballistics during the 2007 eruption (Kilgour et al. 2010). The relative uncertainty of clast density is small in comparison to that regarding other input parameters, and has a comparatively small influence on resulting ballistic trajectories.
We randomly varied clast diameter, initial velocity and ejection angle, depending on the eruption scale. We only consider ballistic projectiles larger than 20 cm in diameter, because smaller clasts are likely to be affected by the eruptive column and their dispersal can therefore not be accurately simulated by our ballistic ejection model (AlatorreIbargüengoitia et al. 2016). The largest ballistic clasts of the 2007 eruption were 2 m in diameter, and eruption fragmentation is assumed to follow a power law (Girault et al. 2014). Hence we invoke three clast diameters with relative proportions based on a power law distribution: 60% of all simulated clasts are 20 cm, 30% are 1 m and 10% of clasts are 2 m in diameter.
Eruptions at Ruapehu are scaled according to the distance covered by their deposits. For each scale of eruption (Table 1) we used this maximum distance reached by volcanic ejecta to infer a maximum eruption velocity of the ballistic clasts. The maximum velocity for ballistics produced by eruptions of Scale 4 is also constrained by Kilgour et al. (2010), who inferred a maximum velocity of 135 m/s for the 2007 phreatic eruption. The minimum velocity for Scale 2 eruptions is based on the minimum distance necessary for ejecta to make it out of the lake. The minimum ejection velocities for Scale 3 and 4 are set to the maximum velocity of Scale 2 and 3, respectively. Within these intervals (see Table 3), we sample a uniform distribution of velocities. In doing so, we assume that exit velocity and grain size are independent parameters, which of course is a simplification. More sophisticated fluid dynamic models of the conduit acceleration phase take into account the effect of the carrier flow field on the eruptive mixture more accurately (e.g. de’Michieli Vitturi et al. 2010), but this is beyond the scope of this study.
The ejection angle depends on the erupting vent: Ruapehu Crater Lake hosts two eruptive vents, a northern and a central one (Christenson 2000). While the central vent is upright, the northern vent is inclined to the North at 45–60° from the vertical (Fig. 1, Kilgour et al. (2010)). To obtain both an ejection and a direction angle, we sample the von MisesFisher distribution (e.g. Fisher et al. 1993), a probability distribution on a sphere, using a MATLAB®; procedure following instructions in Wenzel (2015). The distribution has two parameters, the mean direction and the concentration parameter κ. For the central vent, the mean direction is upright with κ=5; for the northern vent we use a vector inclined 52.5° (from the vertical) to the North as the mean direction and κ=10, since explosions from this vent tend to be more directed. Figure 5 illustrates the distributions for the two vents. Sampling the von MisesFisher distribution delivers unit vectors on the whole sphere, from which we infer both the azimuth and ejection angle. This gives ejection angles between 90 and 90° from the horizontal. Only clasts that make it out of the crater lake are of interest to the hazard assessment for the summit. Furthermore, the total number of ejected clasts that make it onto the summit is needed to calculate the impact density and data for this number are available, while no data are available for clasts that remain in the lake. We therefore only simulate trajectories of clasts that make it out of the crater lake. This requires the ejection angle to be larger than a minimum value, which depends on the vent location (i.e. North or central) and ejection direction, and lies between 17 and 36°. To truncate the distribution, we therefore test whether the obtained ejection angles are within the correct interval, and if not, then the angle is discarded and we resample the distribution for a new angle until it is in the desired interval.
Step 2: Calculate ballistic trajectories and infer impact distances
Ballistic trajectories accounting for drag effects were calculated by a RungeKutta procedure with the program “Eject!” (Mastin 2001; 2011), based on the clast’s size, mass, ejection angle, initial velocity and other input parameters (Table 3). We translated the Visual Basic code to MATLAB®; in order to call it automatically and easily combine it with our evaluation scripts. The calculated trajectory is then compared with an interpolated topography function, depending on the random ejection direction of the clast (North, East, South or West), to determine its landing point. Clastclast interactions as presented in Tsunematsu et al. (2014) and the effects of tailwind are ignored in these simulations.
Step 3: Derive the probability for one ejected clast to reach a certain area
The topography of the summit and the location and eruption behaviour of the northern vent have a strong directiondependent influence on impact distances and therefore impact probabilities. To take this into account in the evaluation of results, we divide the areas A _{ n } into 4 pitch rings A _{ n−dir }, with dir being the azimuth (North, East, South or West) (see Fig. 4). For each eruption scale and both vents, we ran 10,000 trajectory simulations that were used to infer the probability P _{ n−dir }(c) of one ejected clast to reach the area A _{ n−dir }. We count a success for A _{ n−dir } if a clast’s ejection direction is dir and its impact distance is between radius r _{ n } and r _{ n+1} (see Fig. 4). The number of successes is then divided by 10,000 to obtain probabilities (conditional to an eruption scale and vent location). Subsequently, we average the probabilities obtained for the northern and central vent, resulting in probabilities depending only on eruption scale and respective area. This implies that eruptions occur from the two vents with an equal probability, reflecting our complete ignorance on whether one of them is more active. Unfortunately, there are no informative data available.
Step 4: Derive a probability density function for impact densities
The above described procedure provides the probability for an ejected clast to reach a certain area on the summit. We are, however, interested in the total number of ballistic impacts per area and this depends on the number of clasts that get ejected onto the summit. The more clasts are ejected, the higher are expected impact densities. From the derived probabilities, one can calculate the ballistic impact density F _{ n−dir } for area A _{ n−dir } with
with m being the number of erupted clasts that make it out of the crater. There are few data available for m in phreatic eruptions. The only comprehensive study, presented by Fitzgerald et al. (2014), estimates 15,000 ejected clasts (outside of the crater) during the Upper Te Maari eruptions, Tongariro, in 2012. We sample a uniform distribution between 10,000 and 20,000, to obtain 10,000 random values for m and derive a mean value and standard deviation for impact densities in the different areas. Counting how often the derived impact density overcomes a certain threshold in the simulations and dividing by 10,000 gives the prior probabilities used as input in nodes 7&8 in PyBetVH.
Results part two: probabilistic hazard assessment for Ruapehu
Results from Eject! simulations
From the Eject! simulations, we obtained the probabilities P _{ n−dir }(c) for a clast reaching an area A _{ n−dir }, given the eruption of a certain scale. These probabilities vary with distance from the crater as well as the cardinal direction, which is shown for distances up to 1 km away from the lake centre in Figs. 6, 7 and 8. The variation with ejection direction is due to the varying topography of the summit and the influence of the northern vent that produces directed blasts. From these probabilities, we calculated the impact density in the respective areas (Figs. 6, 7 and 8), which can be multiplied by 10 m^{2} to obtain the average number of clasts affecting a person. The hazard is substantially greater on the northern side of the lake. For example, the expected impact density on the northern shore during an ongoing eruption of Scale 3 is about 0.03 clasts/m^{2} when one is located less than 400 m from the lake centre, and 0.025 clasts/m^{2} when one is located less than 550 m from the lake centre  compared to just about 0.01 and 0.001 clasts/m^{2}, respectively, on the shores in other directions. Comprehensive results from the simulations, including all calculated mean impact densities, their standard deviation and the probabilities for the ballistic impacts in areas of interest to overcome intensity thresholds, can be found in the Additional file 1 to this paper.
PyBetVH results
Using PyBetVH, we combined the information from the phreatic eruption dataset with prior knowledge and the Eject! simulation results. Since we have a high number of recorded eruptions, the probability calculated by PyBetVH for an eruption of a certain scale occurring within the next month is very similar to the frequency directly calculated from the database: The probability for an eruption occurring within the next month is 10.8%, the conditional probabilities for the eruption being of Scale 1,2, 3 or 4 are 48.0, 30.4, 14.0 and 7.6%, respectively.
The final products from PyBetVH are the monthly probability and hazard maps for ballistic impact densities. Figure 9 shows hazard maps for the ballistic impact density at a probability of 1% conditional to an eruption of Scale 2, 3 and 4, respectively. Most significant is the extended hazardous area on the northern shore for eruptions of Scale 4, where the impact density is 0.01 clasts/m^{2} up to 2 km away from the vent. The maps also indicate infrastructure on the summit, which are likely whereabouts of people and can be damaged by ballistics. The location of the dome shelter is within the most hazardous area for Scale 3 and 4 eruptions; ballistics of a Scale 4 eruption can even reach ski lifts on Ruapehu’s northern flank. Figure 10 shows the absolute hazard map for ballistic impact density at a probability of 1%, which is up to 0.01 clasts/m^{2} everywhere up to 400 m away from the lake centre, i.e. about 150 m from the lake shore.
Figure 11 shows the probabilities for areas to experience a ballistic impact density F>0.01 clasts/m^{2} conditional to an eruption of Scale 2, 3 and 4, respectively. This impact density would correspond to one clast per 100 m^{2}, or a person being hit by 0.1 clasts on average, so the chance to be hit is 10%. For Scale 2, this probability is 1 on the northern and southern lake shore up to 400 m from the lake centre, and between 0.1 and 0.6 on the eastern and western shore, decreasing with distance. For Scale 3 eruptions, the probability is 1 up to 600 m away from the centre on the northern shore, while the other shores show probabilities between 0.125 and 0.75, again decreasing with distance from the vent. Probabilities directly at the lake shore for Scale 4 eruptions are generally slightly lower, however, the area with probabilities larger than 0.1 extends significantly farther to the north, up to 2 km away from the vent. Figure 12 shows the absolute probability for areas to be affected by a ballistic impact density F>0.001 clasts/m^{2} within the next month, which is larger than 1% up to 900 m away from the lake centre on the northern shore, and up to 500 m away on the other shores. It increases to about 6% closer to the lake.
Discussion
Limitations of the presented methodology
The presented results for Ruapehu are subject to several uncertainties. First, the dataset is very likely not complete, especially for eruptions of Scale 1 or 2. Additional errors originate from the lack of robust examination and determination of magmatic vs phreatic explosions  not all have been studied in sufficient detail with regards to amounts and origin of juvenile material. The possible misinterpretation of terms used in the catalogue can lead to a false inclusion or exclusion of an event in the phreatic database.
The most substantial errors in the statistical assessment of ballistic hazards stem from our ignorance of any wind effects when calculating the ballistic trajectories. Strong winds can occur on Ruapehu and significantly influence ballistic flight paths, especially of smaller clasts. We focused here on clasts larger than 20 cm to ensure that the effects of wind and the eruption column are of second order. Smaller clasts are very likely affected even by weak winds and the eruptive column, but they can still cause severe damage when colliding with a person and should be included in future studies.
We furthermore invoke several simplifications, such as an arbitrary radius of reduced drag and an ejection velocity independent of the grain size. More sophisticated numerical models, e.g. as presented in de’Michieli Vitturi et al. (2010), take processes during the acceleration phase in the conduit into account and more accurately simulate the effect of the carrier flow field on ballistic trajectories. Future studies should make use of these more advanced models.
We also assume that the grain size distribution of ballistic projectiles can be represented by 3 distinct clast diameters weighted according to a power law, which is likely an oversimplification. However, to date, no comprehensive total grain size distribution of ballistic deposits could be determined. An important error source in our calculations of expected impact densities is the uncertainty related to the number of ejected clasts m. Only one comprehensive study on a hydrothermal eruption provides an estimate on this number (Fitzgerald et al. 2014). Additionally, since m inevitably depends on eruption dynamics and cratertopography, using a different volcano as an analogue is problematic. However, no better data are available at the moment and we therefore have to view our results under this caveat.
The created hazard maps are of course somewhat dependent on the chosen areal grid resolution, since the evaluation of simulation results necessarily involved averaging the discontinuous ballistic impacts over areas of interest on the summit. The overall spatial resolution of the expected impact density on the summit is more accurate the finer the grid. To test this influence, we decreased the areas A _{ n−dir } by adding another radius halfway between all of the previously defined r _{ n } and recalculated impact densities from the ballistic simulations. This affected results in the most proximal ring (from the lake shore up to 150 m radial distance), where the mean impact densities in the halved ring deviated significantly from those calculated in the original, larger ring. However, the resulting impact density still correctly mirrors the averaging of ballistic impacts and adequately gives the number of expected clasts in the defined areas. For hazard assessment purposes, it is sufficient to know the average number of expected ballistics in the defined areas.
Comparison of simulation results with the 2007 phreatic eruption
To assess validity of our ballistic simulations, we can compare simulated impact densities for Scale 4 eruptions with the observed values of the 2007 phreatic eruption, reported in Kilgour et al. (2010). This northward directed blast produced a ballistic apron on the northern lake shore up to a distance of 2 km, with an almost homogeneous impact density of 0.014 clasts/m^{2}. The calculated impact density from simulations using random ejection angles corresponding to a northern vent eruption is negligible on the other shores, and lies between 0.009 and 0.02 clasts/m^{2} on the northern shore up to a distance of 850 m from the lake, which is a reasonably good agreement with the observations. However, at larger distances, the simulations predict significantly lower impact densities (<0.005 clasts/m^{2}) than observed. According to the isopleth maps in Kilgour et al. (2010), the distal ballistic field is almost entirely small clasts. While we have not modelled the trajectories of clasts <20 cm, the reported impact density takes the smaller clasts into account. Southerly winds were strong during the 2007 eruption and could therefore have transported the smaller clasts, which are most affected by wind, to the distal field. Hence, the lack of wind as well as omitting clasts smaller than 20 cm in our simulations are probably the main causes of the discrepancy in the distal impact densities. Different ejection velocities could be another contributing factor. If smaller ballistic clasts are ejected at higher velocities, they are transported to larger distances, which is not accounted for in our simulations. However, while they cannot reproduce the finergrained, distal ballistic field of an eruption in windy conditions, our simulations likely produce valid results for the coarsergrained part of the ballistic field.
Hazard implications
Previous studies on ballistic hazards have deterministically calculated the maximum distance reached by ballistics in different eruption scenarios, in order to define safety zones around a volcano (e.g. AlatorreIbargüengoitia and DelgadoGranados 2006; AlatorreIbargüengoitia et al. 2016; Konstantinou 2015). Probabilistic approaches for ballistic hazard assessment have been proposed only recently (Fitzgerald et al. 2014; Biass et al. 2016). They provide a more accurate and detailed assessment of the hazard due to ballistic projectiles and are therefore of increased use for longterm riskreduction strategies on volcanoes. When producing hazard maps, the chosen measure of hazard intensity is an important factor. While we followed Fitzgerald et al. (2014) in using impact density as a proxy for hazard intensity, Biass et al. (2016) focused on the kinetic energy on impact, which is a commonly used intensity measure for ballistics (e.g. Tsunematsu et al. 2016). Hail storms are a natural hazard that is somewhat comparable to volcanic ballistics, in the sense that the threat consists of discontinuous impacts over an area. While hail storms are classified based mostly on grain size and impact energy, it has long been acknowledged that other factors such as the number of hailstones, which translates into an impact density, also affect the severity of hail hazards (Webb et al. 2009). Impact energy is the more relevant hazard proxy for studies regarding the vulnerability of built infrastructure, because only projectiles of a sufficient energy are capable of damaging buildings and penetrating roofs. However, the only building in the proximity of Ruapehu’s crater lake is dome shelter, a mountain refuge, and only the most distal ballistics pose a threat to some ski lifts. Similar to the study of Fitzgerald et al. (2014), our main interest is the threat to unprotected people on the summit, which is frequented by hikers and skiers. Any encounter with a ballistic projectile is a health hazard and potentially lethal for a person. Therefore, impact density is the more relevant hazard proxy in this case as it directly translates into the probability to be hit and injured by a projectile via multiplication with the hazardous area per clast, as outlined in Fitzgerald et al. (2014).
We can derive a few important implications for eruption forecasting and hazard assessment at Ruapehu volcano in particular, and possibly similar systems prone to phreatic eruptions. Based on our catalog, magmatic phases at Ruapehu are preceded by an increase in phreatic activity, with up to a 10fold increase in explosion frequency. However, the duration of phreatic activity phases and the rate of phreatic events are not clearly correlated with the timing of magmatic events. No statistically significant relationship between phreatic eruption scale and duration of relative quiescence between two phreatic phases can be inferred from our data. Interestingly, Phillipson et al. (2013) also found no statistically significant relationship between the length of quiet periods and the duration of following unrest periods at volcanoes. However, the phreatic events tend to cluster, so one can expect more phreatic explosions to follow once the system is reactivated. Since shorter IET are generally more likely than longer periods of quiescence, this also suggests that a bigger phreatic event is more likely within a few years after a phreatic explosion than during a prolonged quiet phase of the volcano. It remains for future research to determine whether phreatic eruption catalogues at other systems  with and without crater lakes  show similar statistics and to determine the physical reason behind the similarities.
The overall probability for a phreatic eruption to occur at Ruapehu within the next month is about 10%  although in the current state of the volcano, the probability seems to be lower (about 3%). The absolute probability for areas within 1 km from the lake centre to be affected by a ballistic impact density larger than 0.001 clasts/m^{2} within the next month is up to 6%. Conditional to a phreatic eruption, some areas on the summit have a probability of 100% to experience an impact density larger than 0.01 clasts/m^{2}  these areas are larger, the stronger the eruption. Following estimations in Fitzgerald et al. (2014), this impact density implies that a person standing in these summit areas roughly has a 10% chance of being hit by a ballistic large enough to cause injury. To put these values into perspective, impact densities during the lethal 2014 Mount Ontake eruption exceeded 0.4 clasts/m^{2} up to 750 m distance from the erupting vent, and were between 0.1 and 0.4 clasts/m^{2} in the zone up to 850 m distance (Kaneko et al. 2016). These values translate to a 100% chance to be hit by a ballistic, explaining the large number of casualties during this event. Based on the results of our simulations and data from the 2007 eruption, the threat for injury by ballistic projectiles during a phreatic eruption on Ruapehu seems to be an order of magnitude less than on Ontake. However, it is still a nonnegligible danger for people on the summit.
Generally, the hazard is higher on the northern shore of the lake  which is important, because the dome shelter on the summit is built on that side of the lake and the presence of people is therefore more likely in that area. While ski lifts are only threatened by ballistics from a Scale 4 eruption, with an expected impact density of less than 0.006 clasts/m^{2} during an eruption, the dome shelter lies within the most hazardous zone for Scale 3 and 4 eruptions, where the expected impact density is larger than 0.01 clasts/m^{2} (Fig. 9). There is an absolute probability of more than 2% that the shelter is affected by ballistics within the next month (Fig. 12). However, our hazard and probability maps show that the hazard due to ballistic ejecta on the eastern and western lake shore is mostly constrained to a ring of a few hundred metres around the crater lake. If a precursory signal for a phreatic explosion is identified, even a short warning time could therefore be sufficient to relocate to safer areas.
In this study, we provide a probabilistic hazard assessment for ballistic impacts. Expanding our work for a risk assessment for Ruapehu’s summit would require an additional vulnerability and exposure analysis that needs to be combined with the hazard assessment. For example, the average number of people on the summit per month is an important parameter in determining the probability of casualties on Ruapehu. It is likely that there is a seasonality in this risk, because more hikers visit the summit during the summer months. The time spent near the lake on the other hand is dominant in determining the risk for an individual to be affected by a ballistic clast. Our study focused on a hazard analysis with proxies relevant for human injuries. In order to study the risk for built infrastructure, such as the mountain refuge or the ski lifts, it would be necessary to (i) use a different hazard proxy, i.e. impact energy and (ii) determine the energy thresholds necessary to damage this infrastructure, as has been done for example in Jenkins et al. (2014) and Pomonis et al. (1999). Our hazard and probability maps provide the first and important step towards comprehensive and quantitative risk analyses.
Possible reasons for clustering of phreatic explosions and potential use for eruption forecasting
Ruapehu Crater Lake shows thermal cycling, with temperatures varying irregularly between 10 and 60 °C, and increased phreatic as well as phreatomagmatic activity typically coincides with periods of higher lake temperatures (e.g. Christenson and Wood 1993). The elevated heat flow through the hydrothermal system indicates a more shallowseated magma source that could enhance eruptive activity. Since 1960, when regular monitoring began, only 3 eruptions have been observed during periods of decreased heat flow. These have been interpreted to result from pressure build up in a sealed hydrothermal system (e.g. Christenson et al. 2010).
There is no obvious correlation between the cumulative eruption number and the observed lake temperatures (Fig. 13). To further test whether there is a link between lake temperature and phreatic eruptions, we define temperature bins in 5 °C intervals and compare the relative frequency of prephreatic eruption temperatures to a baseline dataset (Fig. 14). There is a wide range of preeruptive temperatures, but the mean preeruption temperature is significantly higher than the mean recorded lake temperature, and more than half of eruptions follow lake temperatures higher than 40 °C. This is in line with findings of previous studies, e.g. Hurst (1981). While there might be some sampling bias, this supports the theory that phreatic explosions are more frequent during periods of a more shallowseated magma, which is consistent with an increased frequency prior to (or indeed within) magmatic phases. The differing heat input from the magma column at depth is a possible reason for the clustering of phreatic explosions.
The possible correlation between lake temperature and eruption frequency is important information for volcano monitoring and eruption forecasting purposes. An increase in lake temperatures, which are now continuously monitored at Ruapehu, can be a warning sign of an impending eruption. To provide a more quantitative estimate, Fig. 15 a shows the frequency of eruption occurrence within 2 months after a temperature measurement was performed (this time window was chosen because temperature measurements were only made on a roughly monthly basis in the past). Eruption frequencies (ef) of 0.3 for lake water temperatures below 15 °C are significantly higher than those for the temperature range between 15 to 30 °C (ef<0.2). However, increased frequencies of between 0.4 and 0.8 are also evident for lake temperatures higher than 40 °C. This indicates that the eruption frequencies follow a bimodal distribution that most likely stems from the two mechanisms leading to phreatic eruptions: (i) increased heat flow due to more shallowseated magma, and (ii) pressure buildup in a sealed hydrothermal system. The distributions for different eruption scales are shown in Fig. 15 b. Distributions for eruptions of Scales 1 and 3 only show the increased eruption frequency for hotter lake temperatures, while the distribution for Scale 2 eruptions indicates bimodality. Although its significance is questionable due to few data points, the distribution for Scale 4 eruptions suggests a higher eruption frequency for colder lakes. In summary, these graphs show that lake temperatures both significantly above and below the longterm background level can indicate an increased eruption probability.
Conclusions and outlook
Our results shed light on the underlying statistics of phreatic eruptions at Ruapehu and possibly elsewhere. The most important conclusions for Ruapehu are:

1.
Phreatic events do not follow a Poisson process but tend to cluster, which is possibly due to increased heat input and hence increased activity in periods of a more shallowseated magma column.

2.
Larger events are most likely to follow other events after a short IET, rather than to occur during a longer period of quiescence.

3.
Magmatic episodes are likely preceded by an increase in the number of phreatic events.

4.
The probability of a phreatic eruption occurring at Ruapehu Crater Lake within the next month is about 10% and conditional to a phreatic eruption, large areas are likely to experience an impact density sufficient to threaten people present on the summit.

5.
At any time, there is a nonnegligible probability for the summit area to be affected by ballistics. The hazard is especially high on the northern shore, where infrastructure exists.

6.
If the crater lake is significantly warmer or colder than the background temperature, the probability for an impending eruption is increased.
To improve our understanding of recurrence statistics, future research should test whether the data in the discussed phreatic eruption catalogue can be described with a CoxHawkes process, a probability distribution that describes cluster processes. Future work concerning the ballistic hazard on Ruapehu’s summit needs to consider the effect of wind on ballistic trajectories by incorporating information on the statistical distribution of wind velocities and directions. To further improve the hazard assessment, other hazards from phreatic eruptions should be included in the investigations, most importantly lahars and base surges. Finally, the accuracy of ballistic trajectory simulations and resulting hazard assessments will significantly benefit from a better understanding of processes during ballistic ejection, as well as ground truthing in the field of poorly constrained parameters such as impact densities and grain size distributions.
This study highlights the necessity to better understand steamdriven explosions, which often occur with little or no detectable precursors or advance warning. Statistical approaches developed in this study should support local hazard assessment at Ruapehu and provide a possible framework for similar studies at other volcanoes.
Appendix
Relationship between IET and eruption size
Figure 16 shows IET versus the scale of the eruption preceding or following the quiet interval, respectively. In both cases, no clear relationship is visible and linear regression analysis confirms that no linear fit to the data can be proven: The pvalue for the Ftest for a linear relationship between scale and IET is 0.76 (preceding scale) and 0.84 (following scale), respectively, so the null hypothesis of no relationship cannot be rejected at the 0.05 significance level  meaning that there is no statistically significant linear relationship. From visual inspection, it seems that only eruptions of smallermedium scales precede or follow a longer IET, while all scales occur for shorter IET. However, this is very likely a pure sampling effect: generally, longer IET are less likely, as are larger eruptions scales. In combination, this leads to few large eruptions corresponding to a long IET. Linear regression analysis shows that the hypothesis of no linear relationship between the maximum IET and corresponding eruption scale cannot be rejected at the 0.05 significance level (pvalue of 0.42 and 0.22, respectively). We also tested for linear relationships between eruption scales and log(IET), but again, the linear models are not significant at the 0.05 significance level (pvalues are >0.05).
Abbreviations
 A _{ n } :

Radial areas (rings) surrounding the crater lake, defined by r _{ n }
 A _{ n− dir } :

Pitch ring (quarter of A _{ n }), with dir being the azimuth (N, E, S, W)
 CV:

Coefficient of variation
 ef :

Eruption frequency
 F :

Ballistic impact density (clasts per unit area)
 F _{ n−dir } :

Ballistic impact density for area A _{ n− dir }
 IET:

Intereruptive time (i.e. time between two eruptions)
 m :

Number of ejected clasts that make it out of the crater lake during a phreatic eruption
 N :

Number of eruptions
 PVHA:

Probabilistic volcanic hazard assessment
 P _{ n−dir }(c):

Probability of one ejected clast to reach the area A _{ n− dir }
 r _{ n } :

Minimum radial distance of A _{ n } from the lake centre
 VEI:

Volcanic Explosivity Index
 κ :

Concentration parameter of the von MisesFisher distribution
 λ :

Number of equivalent data
 σ :

Standard deviation
References
AlatorreIbargüengoitia, MA, DelgadoGranados H. Experimental determination of drag coefficient for volcanic materials: calibration and application of a model to Popocatépetl volcano (Mexico) ballistic projectiles. Geophys Res Lett. 2006:;33(11). doi:10.1029/2006GL026195.
AlatorreIbargüengoitia, MA, MoralesIglesias H, RamosHernández SG, JonSelvas J, JiménezAguilar JM. Hazard zoning for volcanic ballistic impacts at El Chichón Volcano (Mexico). Nat Hazards. 2016; 81(3):1733–44.
Barberi, F, Bertagnini A, Landi P, Principe C. A review on phreatic eruptions and their precursors. J Volcanol Geoth Res. 1992; 52(4):231–46.
Biass, S, Falcone JL, Bonadonna C, Di Traglia F, Pistolesi M, Rosi M, Lestuzzi P. Great Balls of Fire: A probabilistic approach to quantify the hazard related to ballistics  A case study at La Fossa volcano, Vulcano Island, Italy. J Volcanol Geoth Res. 2016; 325:1–14.
Browne, PRL, Lawless JV. Characteristics of hydrothermal eruptions, with examples from New Zealand and elsewhere. EarthSci Rev. 2001; 52:299–331.
Carniel, R, Jolly A, Barbui L. Analysis of phreatic events at Ruapehu volcano, New Zealand using a new SOM approach. J Volcanol Geoth Res. 2013; 254:69–79.
Christenson, B. Geochemistry of fluids associated with the 1995–1996 eruption of Mt. Ruapehu, New Zealand: signatures and processes in the magmatichydrothermal system. J Volcanol Geoth Res. 2000; 97:1–30.
Christenson, B, Wood C. Evolution of a venthosted hydrothermal system beneath Ruapehu Crater Lake, New Zealand. B Volcanol. 1993; 55:547–65.
Christenson, B, Reyes A, Young R, Moebis A, Sherburn S, ColeBaker J, Britten K. Cyclic processes and factors leading to phreatic eruption events: Insights from the 25 September 2007 eruption through Ruapehu Crater Lake, New Zealand. J Volcanol Geoth Res. 2010; 191:15–32.
Cox, D, Lewis P. The Statistical Analysis of Series of Events. New York: Methuen; 1966.
Dibble, RR. Seismic Power Recordings during Hydrothermal Eruptions from Ruapehu Crater Lake in April 1968. J Geophys Res. 1969; 74(27):6545–51.
Dvorak, J. Mechanism of explosive eruptions of Kilauea Volcano, Hawaii. B Volcanol. 1992; 54:638–45.
Fischer, TP, Ramírez C, MoraAmador RA, Hilton DR, Barnes JD, Sharp ZD, Le Brun M, de Moor JM, Barry P, Füri E, Shaw AM. Temporal variations in fumarole gas chemistry at Poás volcano, Costa Rica. J Volcanol Geoth Res. 2015; 294:56–70.
Fisher, N, Lewis T, Embleton B. Statistical analysis of spherical data. Cambridge: Cambridge University Press; 1993.
Fitzgerald, RH, Tsunematsu K, Kennedy BM, Breard ECP, Lube G, Wilson TM, Jolly AD, Pawson J, Rosenberg M, Cronin SJ. The application of a calibrated 3D ballistic trajectory model to ballistic hazard assessments at Upper Te Maari, Tongariro. J Volcanol Geoth Res. 2014; 286:248–62.
Girault, F, Carazzo G, Tait S, Ferrucci F, Kaminski É. The effect of total grainsize distribution on the dynamics of turbulent volcanic plumes. Earth Planet Sci Lett. 2014; 394:124–34.
Glicken, H, Nakamura Y. Restudy of the 1888 eruption of Bandai volcano, Japan. In: Proc Kagoshima Int Conf Volcanoes. Japan: Kagoshima International Conference on Volcanoes: 1988. p. 392–95.
Hurst, A. Pattern recognition applied to temperature and chemical measurements of Crater Lake, Ruapehu. Tech. Rep. 181, Wellington, Geophysics Division, Dept. of Scientific and Industrial Research; 1981.
ICAO. Manual of the ICAO Standard Atmosphere: extended to 80 kilometres (262 500 feet), 3rd edn. International Civil Aviation Organization. Quebec: International Civil Aviation Organization; 1993.
Jenkins, S, Spence R, Fonseca J, Solidum R, Wilson T. Volcanic risk assessment: Quantifying physical vulnerability in the built environment. J Volcanol Geoth Res. 2014; 276:105–20.
Jolly, AD, Sherburn S, Jousset P, Kilgour G. Eruption source processes derived from seismic and acoustic observations of the 25 September 2007 Ruapehu eruption  North Island, New Zealand. J Volcanol Geoth Res. 2010; 191:33–45.
Kaneko, T, Maeno F, Nakada S. 2014 Mount Ontake eruption: characteristics of the phreatic eruption as inferred from aerial observations. Earth Planets Space. 2016; 68(1):1–11.
Kato, A, Terakawa T, Yamanaka Y, Maeda Y, Horikawa S, Matsushiro K, Okuda T. Preparatory and precursory processes leading up to the 2014 phreatic eruption of Mount Ontake, Japan. Earth Planets Space. 2015:;67(111). doi:10.1186/s406230150288x.
Kilgour, G, Manville V, Della Pasqua F, Graettinger A, Hodgson K, Jolly G. The 25 September 2007 eruption of Mount Ruapehu, New Zealand: Directed ballistic surtseyan jets, and iceslurry lahars. J Volcanol Geoth Res. 2010; 191:1–14.
Kilgour, G, Saunders K, Blundy J, Cashman K, Scott B, Miller C. Timescales of magmatic processes at Ruapehu volcano from diffusion chronometry and their comparison to monitoring data. J Volcanol Geoth Res. 2014; 288:62–75.
Konstantinou, K. Maximum horizontal range of volcanic ballistic projectiles ejected during explosive eruptions at Santorini caldera. J Volcanol Geoth Res. 2015; 301:107–15.
Le Guern, F, Bernard A, Chevrier RM. Soufriere of Guadeloupe 1976–1977 eruption—Mass and energy transfer and volcanic health hazards. Bull Volcanol. 1980; 43(3):577–93.
Maeda, Y, Kumagai H, Lacson R, Figueroa II M, Yamashina T, Ohkura T, Baloloy A. A phreatic explosion model inferred from a very long period seismic event at Mayon Volcano, Philippines. J Geophys ResSol Ea. 2015; 120(1):226–42.
Marzocchi, W, Sandri L, Gasparini P, Newhall C, Boschi E. Quantifying probabilities of volcanic events: The example of volcanic hazard at Mount Vesuvius. J Geophys Res. 2004; 109:B11, 201. doi:10.1029/2004JB003155.
Marzocchi, W, Sandri L, Selva J. BET_VH: a probabilistic tool for longterm volcanic hazard assessment. B Volcanol. 2010; 72:705–16. doi:10.1007/s0044501003578.
Mastin, L. Thermodynamics of gas and steamblast eruptions. B Volcanol. 1995; 57:85–98.
Mastin, LG. A simple calculator of ballistic trajectories for blocks ejected during volcanic eruptions: U.S. Geological Survey OpenFile Report 0145, 16p; 2001. https://pubs.usgs.gov/of/2001/0045/.
Mastin, L. eject model. 2011. https://vhub.org/resources/455. Accessed 25 Sept 2015.
Miura, K, Ban M, Ohba T, Fujinawa A. Sequence of the 1895 eruption of the Zao volcano, Tohoku Japan. J Volcanol Geoth Res. 2012; 247–248:139–57.
Mordret, A, Jolly A, Duputel Z, Fournier N. Monitoring of phreatic eruptions using Interferometry on Retrieved CrossCorrelation Functions from Ambient Seismic Noise: Results from Mt. Ruapehu, New Zealand. J Volcanol Geoth Res. 2010; 191:46–59.
Nairn, IA, Wood CP, Hewson CAY. Phreatic eruptions of Ruapehu: April 1975. J Geol Geophys. 1979; 22(2):155–73.
Neri, A, Macedonio G, Gidaspow D. Phreatic explosion hazard assessment by numerical simulation. Phys Chem Earth Pt A. 1999; 24(11):989–95.
Newhall, C, Albano S, Matsumoto N, Sandoval T. Roles of groundwater in volcanic unrest. J Geol Soc Philippines. 2001; 56:69–84.
NIWA. Climate Data and Activities. 2014. http://www.niwa.co.nz/educationandtraining/schools/resources/climate. Accessed 22 May 2014.
Ogiso, M, Yomogida K. Migration of tremor locations before the 2008 eruption of Meakandake Volcano, Hokkaido, Japan. J Volcanol Geoth Res. 2012; 217:8–20.
Phillipson, G, Sobradelo R, Gottsmann J. Global volcanic unrest in the 21st century: An analysis of the first decade. J Volcanol Geoth Res. 2013; 264:183–96. doi:http://dx.doi.org/10.1016/j.jvolgeores.2013.08.004.
Pomonis, A, Spence R, Baxter P. Risk assessment of residential buildings for an eruption of Furnas Volcano, Sao Miguel, the Azores. J Volcanol Geoth Res. 1999; 92(1):107–31.
Rouwet, D, Sandri L, Marzocchi W, Gottsmann J, Selva J, Tonini R, Papale P. Recognizing and tracking volcanic hazards related to nonmagmatic unrest: a review. J Appl Volcanol. 2014:;3(17). doi:10.1186/s1361701400173.
Scott, B. A revised catalogue of Ruapehu volcano eruptive activity: 1830–2012. GNS Science Report 2013/45. New Zealand: GNS Science; 2013.
Tonini, R, Sandri L, Thompson M. PyBetVH: A Python tool for probabilistic volcanic hazard assessment and for generation of Bayesian hazard curves and maps. Comput Geosci. 2015; 79:38–46.
Tonini, R, Sandri L, Selva J, Marzocchi W. PyBetVH. 2016. https://vhub.org/resources/betvh. Accessed 28 Sept 2015.
Tsunematsu, K, Chopard B, Falcone JL, Bonadonna C. A numerical model of ballistic transport with collisions in a volcanic setting. Comput Geosci. 2014; 63:62–9.
Tsunematsu, K, Ishimine Y, Kaneko T, Yoshimoto M, Fujii T, Yamaoka K. Estimation of ballistic block landing energy during 2014 Mount Ontake eruption. Earth Planets Space. 2016:;68(1). doi:10.1186/s4062301604638.
Venzke E, (ed).Global Volcanism Program, 2013. Volcanoes of the World, v. 4.4.1.; 2013. Smithsonian Institution. Downloaded 23 Oct 2015. http://dx.doi.org/10.5479/si.GVP.VOTW42013.
de’Michieli Vitturi, M, Neri A, Esposti Ongaro T, Lo Savio S, Boschi E. Lagrangian modeling of large volcanic particles: Application to Vulcanian explosions. J Geophys ResSol Ea. 2010:;115(B8). doi:10.1029/2009JB007111.
Webb, J, Elsom D, Meaden G. Severe hailstorms in Britain and Ireland, a climatological survey and hazard assessment. Atmos Res. 2009; 93(1):587–606.
Wenzel, J. Numerically stable sampling of the von Mises Fisher distribution on S ^{2} (and other tricks). 2015. On, https://www.mitsubarenderer.org/wenzel/files/vmf.pdf. Accessed 6 Oct 2015.
Acknowledgments
The authors would like to thank Larry Mastin from USGS for providing the Eject! source code and his assistance in translating the code to MATLAB®;.
Funding
The research leading to these results has received funding from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/20072013) under the project NEMOH, REA grant agreement number 289976. Additional funding was provided by the MEDSUV project, under the grant agreement number 308665, and the VUELCO project under grant agreement number 282759, both part of the European Union’s Seventh Framework programme.
NEMOH has provided salary (2 years), research and travel costs for KS, including secondments at INGV Bologna to collaborate with LS and RT. MEDSUV and VUELCO provided a studentship (1 year) and travel costs for KS. NEMOH and VUELCO provided computational facilities for KS.
Availability of data and materials
The eruption catalogue for Ruapehu is published by Scott (2013) and available online. PyBetVH and Eject! are freely available on the VHub cyberinfrastructure (vhub.org), where further information on the software can be found. Results of our clast ejection simulations are uploaded as an additional supporting file (excel spreadsheet).
Authors’ contributions
KS planned and conducted the research, with substantial input and assistance from LS, ACR and JHG. RT provided support with PyBetVH. GK provided expertise concerning Ruapehu Volcano, as well as crater lake temperature data and the basis for Fig. 1. KS drafted the initial version of the manuscript, which has been reviewed and approved by all authors.
Competing interests
The authors declare they have no competing interests.
Author information
Additional file
Additional file 1
Results of ballistic flightpath simulations. (XLSX 55 KB)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Received
Accepted
Published
DOI
Keywords
 Phreatic eruptions
 Ruapehu
 Probabilistic hazard assessment
 PyBetVH
 Ballistic ejecta