Sensitivity to volcanic field boundary

Abstract

Volcanic hazard analyses are desirable where there is potential for future volcanic activity to affect a proximal population. This is frequently the case for volcanic fields (regions of distributed volcanism) where low eruption rates, fertile soil, and attractive landscapes draw populations to live close by. Forecasting future activity in volcanic fields almost invariably uses spatial or spatio-temporal point processes with model selection and development based on exploratory analyses of previous eruption data. For identifiability reasons, spatio-temporal processes, and practically also spatial processes, the definition of a spatial region is required to which volcanism is confined. However, due to the complex and predominantly unknown sub-surface processes driving volcanic eruptions, definition of a region based solely on geological information is currently impossible. Thus, the current approach is to fit a shape to the known previous eruption sites. The class of boundary shape is an unavoidable subjective decision taken by the forecaster that is often overlooked during subsequent analysis of results. This study shows the substantial effect that this choice may have on even the simplest exploratory methods for hazard forecasting, illustrated using four commonly used exploratory statistical methods and two very different regions: the Auckland Volcanic Field, New Zealand, and Harrat Rahat, Kingdom of Saudi Arabia. For Harrat Rahat, sensitivity of results to boundary definition is substantial. For the Auckland Volcanic Field, the range of options resulted in similar shapes, nevertheless, some of the statistical tests still showed substantial variation in results. This work highlights the fact that when carrying out any hazard analysis on volcanic fields, it is vital to specify how the volcanic field boundary has been defined, assess the sensitivity of boundary choice, and to carry these assumptions and related uncertainties through to estimates of future activity and hazard analyses.

Background

Regions of distributed volcanism are often referred to as volcanic fields. As each eruption tends to occur in a different location, they can be further classified as monogenetic volcanic fields and are frequently modelled as spatial or spatio-temporal point processes (Connor and Hill 1995; Marzocchi and Bebbington 2012; Bebbington 2013). Volcanic activity within these regions tends to be predominantly basaltic, with low flux rates. Consequently, volcanic fields tend to grow laterally rather than build up a single edifice in one location (Wood and Shoan 1981), and eruptions typically form cinder cones, tuff cones and rings, maars, lava domes, shield volcanoes, and lava flows (Connor and Conway 2000).

Determination of a distributed volcanic field’s spatio-temporal behaviour is key to understanding its underpinning magmatic processes, evolution, and hazard potential. As these subsurface processes are predominantly unknown, and highly complex, this determination is accomplished via approximation of the known spatio-temporal evolution of the field, and then extrapolated into the future to obtain hazard forecasts. This almost invariably uses spatial or spatio-temporal point processes with model selection and development based on analysis of previous eruption data. The location of future eruptions substantially influences the areas at risk from eruption related hazards. The specific location may also affect the eruption style, for example, the presence of ground-water may lead to an explosive phreatomagmatic eruption (Ross et al. 2011). Thus, the probable location is one of the most important pieces of information for volcanic hazard estimation.

To obtain probable future locations, the first step is to assess the existence of patterns in previous eruption locations and timings. This is accomplished using various exploratory statistics including lineament identification (Wadge and Cross 1988; Lutz and Gutmann 1995) and cluster analyses (Clark and Evans 1954; Hopkins and Skellam 1954; Ripley 1979). For these, a surface region representing the volcanic field must be defined, both for area calculations and Monte Carlo point placement steps. A region must also be defined for the evaluation of spatio-temporal models (Connor and Hill 1995; Cronin et al. 2001), and for assessment of all spatial and spatio-temporal models during verification and validation steps.

The spatial boundaries of a volcanic field are conventionally defined as the extent of the existing lava flows, scoria cones, or other volcanic edifices. These are usually well constrained via geological mapping. This provides a present-day view of the location of past volcanic material from a resource-mapping perspective. However, for assessment of future volcanic activity, this definition is not sufficient as it does not account for the locations of new eruption vents and may confound the locations of previous volcanic eruptions with factors controlling lava distribution. Depending on rheology and underlying topography, lava may extend tens of kilometres from a source vent and accumulate in depressions giving false weight to these regions far from potential eruption locations.

Volcanic fields form above mantle regions where partial melting produces magma, which may be related to low rates of crustal extension (Takada 1994). In these sites, magma supply is insufficient to maintain a single conduit, such that each new magma batch has to find its own path to the surface (Fedotov 1981). Due to the variety of geological conditions under which mantle materials melt to form magmas (e.g., within subduction zones, in areas of weak crustal thinning, or hot-spot upwelling), boundary parameters remain unknown, thus a well-informed geological approach to boundary definition is currently impossible. For forecasting purposes however, the volcanic field boundary has a significant influence on estimates of event recurrence rate and spatial intensity. Where a hard boundary is required (e.g., for clustering statistics, or spatio-temporal models), it is assumed that there is zero probability of an eruption occurring outside of this region. Often, even when a soft boundary is feasible, a hard boundary is imposed for convenience in spatial models (Lindsay et al. 2010; Sandri et al. 2012).

Estimation of the region of potential volcanic vent breakout in a field represents a typical bias-variance trade-off. Underestimating the area may result in future eruptions occurring outside the area of analysis and an overestimate of average vent density, which may in turn bias statistical diagnostics (Clark and Evans 1954; Hopkins and Skellam 1954; Ripley 1979). Overestimation ensures a conservative ‘safe’ solution, but potentially underestimates average vent-density and consequently spatial intensity. It also increases the coverage area over which supplementary data is required for a hazard analysis, thus either increasing costs, or lowering the data resolution.

In this work, various approaches for the definition of volcanic field boundaries are evaluated for spatial intensity estimates at two quite different volcanic fields: Harrat Rahat, Kingdom of Saudi Arabia, and Auckland Volcanic Field (AVF), New Zealand. This paper begins with an introduction to these two volcanic fields selected as case studies. Five boundary options are then outlined which comprise of three commonly applied fitting methods (an ellipse, a rectangle, and a convex hill), alongside two vent-density based suggestions (isotropic and anisotropic kernel smoothing methods). Five geological considerations are then discussed for the two volcanic fields that then form the underlying assumptions and related uncertainties for each fitted boundary. These considerations combined with the five boundary options result in ten potential boundary options for the volcanic field case studies. These provide the basis for the sensitivity analyses for each field and the affect of boundary definition is assessed using exploratory methods (summary statistics) and spatial intensity estimates.

Harrat Rahat

Harrat Rahat is one of at least nineteen volcanic fields within the western Arabian Plate (Brown 1970). The harrat (‘volcanic-field’) has a total erupted volume of ~2000 km³ (Blank and Sadek 1983) from at least 968 identified vents (Runge et al. 2014, Fig. 1a) and is predominantly basaltic with some more evolved compositions. The cause of volcanism across the Arabian Shield is still very much under discussion (Moufti et al. 2013; Rolandone et al. 2013), and although there are other volcanic fields, and a major zone of crustal rifting (the Red Sea rift) proximal to Harrat Rahat, little is known about their respective influences on magma generation within the harrat.

Fig. 1

Volcanic field examples. a Harrat Rahat vent locations, shaded area represents lava flow field; b Harrat Rahat vent density distribution; c Auckland Volcanic Field vent locations; d Auckland Volcanic Field vent density distribution. Vents are shown as black marks, North is up. Vent density contours at [0.05, 0.5, 1, 5, 10, 20, 50, 80] % of maximum density

Erupted products include extensive lava flows (Murcia et al. 2014), as well as scoria and spatter cones, shield volcanoes, tuff rings, and lava domes (Camp and Roobol 1989). Eruptions in Harrat Rahat have been dated between 10 Ma and 1256 AD, with age data broadly suggesting a northward focussing of volcanism with time (Camp and Roobol 1989). A caveat on this, however, is that definitive ages are only available for 3 % of vents, and recent ⁴⁰Ar/³⁹Ar ages suggest there are also older flows in the north, and younger flows in the centre of the field (Moufti et al. 2013).

Previous boundary definitions for Harrat Rahat have included the mapped limits of lava flows (Moufti et al. 2010), a convex hull of the previous known eruption locations (Runge et al. 2014), and a rectangular area of interest within the northern region known as Harrat Al-Madinah (Moufti et al. 2010; El Difrawy et al. 2013).

Auckland volcanic field

The Auckland Volcanic Field (AVF) is a basaltic field located on the North Island of New Zealand. The field has a total erupted volume of at least 1.7 km³ (Kereszturi et al. 2013) from approximately 50 volcanoes (Fig. 1c), Hayward et al. 2011). Volcanic features include lava flows, scoria cones, maars, and tuff rings (Allen and Smith 1994). Eruption order and locations are relatively well constrained (Bebbington and Cronin 2011) with events dated between ~250 and 0.6 ka (Lindsay et al. 2011).

The AVF boundary has been defined as a minimum area ellipse fitted to the outer vent locations (Spӧrli and Eastwood 1997), and was used for evacuation planning (Tomsen et al. 2014). An elliptical boundary was also used by Sandri et al. (2012). This approach was expanded by Bebbington (2013) to include an outer ellipse with an exponentially decaying intensity. Rectangular bounds were employed by Lindsay et al. (2010), and a convex hull with a ~1 km buffer zone was employed for exploratory analyses by Le Corvec et al. (2013a).

Boundary options

While the area is an integral part of the description of a volcanic field, very few studies define how the boundary is represented, thus comparison of exploratory statistics (e.g., average vent densities) between fields is often meaningless, and repeatability of the study impossible. Here, three of the most commonly applied volcanic field boundary definitions will be investigated: (1) convex hull, (2) minimum area ellipse, and (3) minimum area rectangle. Two additional boundary options based on isotropic and anisotropic kernels are also explored. A description of how each of these shapes were fitted to the data is also provided.

Convex hull

For a volcanic field, a convex hull is the smallest convex set that contains all of the vents (any straight line between any two vents must be entirely contained within the hull). Various algorithms can be used to obtain this shape, such as gift-wrapping (Chand and Kapur 1970), or divide-and-conquer (Preparata and Hong 1977), and the majority of statistical packages (e.g., MATLAB, or R) have inbuilt functions for this.

For this reason, and because the fitting of a convex hull is quick, and always renders the same answer, a convex hull is often used as the default shape when applying initial exploratory methods on volcanic fields. One of the first uses of a convex hull in describing volcanic systems was by Zhang and Lutz (1989), for the assessment of the spatial extent in kimberlites. Since then it has been frequently used in volcanic hazard evaluations, e.g., Springerville volcanic field, Arizona (Condit and Connor 1996), Mount Gambier, Newer Volcanic Province, Australia (Bishop 2007) and Auckland Volcanic Field (Le Corvec et al. 2013c). In this study, the convex hull shape was determined via the inbuilt MATLAB function (‘convhull’).

Minimum area ellipse

The main rationale for an elliptical, or curved boundary to a volcanic field is based primarily on thermodynamic arguments (Fedotov 1981). A mantle source zone below a field is likely to be susceptible to partial melting from a relatively stable (in terms of location) heat source (Condit and Connor 1996). These tend to be semi-circular in plan-view, because heat radiates equally in all directions as long as the lithology is broadly homogeneous.

A minimum area ellipse can be fitted to a set of points (in this case, volcanic vents) via a variety of methods such as non-linear least squares (Angell and Barber 1977), maximum likelihood (Kanatani 2005), or bootstrapping (Cabrera and Meer 1996). Ellipse fitting algorithms are slightly more complex to run than a convex hull algorithm because there are five parameters to vary within the chosen optimisation process: the centre of the ellipse (x_o, y_o), the major and minor radii (r_a, r_b), and the orientation angle (θ). In this work, the minimum area ellipse method of Khachiyan (1996) was adopted which greatly reduces the computational time by finding the bounding convex hull first, reducing the number of points required to fit an ellipse.

Minimum area rectangle

A minimum bounding rectangle can be fitted to the points of a convex hull from a data set, similarly to an ellipse. Several relatively straightforward methods are available such as rotating callipers (Toussaint 1983), or minimum perimeter (Freeman and Shapira 1975). Due to their simplicity both in terms of proportional area calculations, and presentation of results, rectangular shapes are often used when analysis of a sub-region of a field is required for example, around an area being assessed as a nuclear waste facility (Connor and Hill 1995).

Although rectangular volcanic field bounds are difficult to justify geologically, similar shapes have been used in areas where there is thought to be some geological structural control. Kear (1964) fitted oblong structures to the Bay of Islands Volcanic Zone, NZ, and a rift-subparallel structure known as the Mekkah-Madinah line was inferred in the Kingdom of Saudi Arabia (Camp and Roobol 1989). This shape also has several advantages over an ellipse: it is easier to fit a minimum area rectangle, and it tends to truncate the field at one or more ends which minimises empty space. In this work, the minimum area rectangle was calculated from the convex hull of the vents, and a rectangle fitted based on the orientations of each edge, the best-fit is then the rectangle corresponding to the minimum area.

Isotropic and anisotropic kernels (Vent Density Contours)

The final two boundary methods included in this work are based on Gaussian isotropic (Eq. 1), and anisotropic (Eq. 2) kernels:

$$ {\lambda}_s=\frac{1}{N}{\displaystyle {\sum}_{i=1}^N\frac{1}{2\pi h} \exp}\left[\frac{{\left\Vert {x}_i-x\right\Vert}^2}{2{h}^2}\right] $$

(1)

$$ {\lambda}_s(x)=\frac{1}{2N\pi \sqrt{\left|H\right|}}{\displaystyle {\sum}_{i=1}^N \exp}\left[-0.5{\left({x}_i-x\right)}^T{H}^{-1}\left({x}_i-x\right)\right] $$

(2)

Where N: total number of eruptions, h: 1D kernel bandwidth, H: 2D kernel bandwidth matrix, |H|: is the determinant of H, and ‖x _i  − x‖²: distance between ith vent at x_i and location x.

These methods have been frequently applied for spatial intensity estimates (Connor and Hill 1995; Martin et al. 2004; Connor and Connor 2009; Kiyosugi et al. 2010; Bebbington and Cronin 2011; Germa et al. 2013; Bebbington 2013) or to represent vent density variation (Connor 1990; Lutz and Gutmann 1995). Despite this potential reflection of spatial variation in magma fertility or persistent melt columns (as represented by high vent density regions, Germa et al. 2013), kernel methods have been suggested (Germa et al. 2013) but are not yet widely used to define a volcanic field boundary. In this study, bandwidth estimation (h, H) is accomplished via SAMSE (‘Squared-asymptotic-mean squared error’, Duong and Hazelton 2003) using the R package ‘ks’ (Duong 2007).

The use of kernel density estimation is consistent with the probabilistic approach taken for intensity estimates by assuming that future volcanic activity is strongly linked to past activity, and thus more fertile regions of the fields (as represented by a greater number of vents) are more highly weighted, pushing the boundary out farther from them. As the exponential decay function (with distance) has infinite support, a user-defined cut-off value is required to define a boundary. Here, this truncation was assumed as the maximum isocontour that includes all vents within the same shape, without containing any empty regions (simply connected).

Geological considerations

Regions of distributed volcanism are the surface representation of highly complex subsurface systems which vary from system to system. Thus, each volcanic field must be considered independently with respect to any potential constraints or known system properties (such as major crustal lithological or fault boundaries). There are several assumptions that need to be made when imposing a boundary and since some may be subjective, it is important that their accompanying uncertainties are propagated through to intensity estimates and hazard forecasts. Assumptions need to address the following:

• (i) Temporal invariance

• (ii) Unobserved (hidden) vents

• (iii) Anomalous vents

• (iv) Geological constraints

• (v) Field maturity

These are considered here one by one.

Temporal invariance

In the case of volcanic field boundaries, temporal invariance suggests a fixed region within which all eruptions have, could have, or will occur. Questions that need to be asked include whether the eruption locations have shifted in a systematic way with time, or if the locus (or loci) of activity has drifted towards, or away from, a particular area.

Shifts in the locus of activity within a volcanic field have been noted by Tanaka et al. (1986) from west to east in the San Francisco volcanic field over the last 2.5 My, and Connor and Hill (1995) found substantial movement in volcanic activity (east to west) within the Yucca Mountain region, Nevada. Kiyosugi et al. (2010) noted a constant location of volcanism within the Abu Monogenetic Volcano Group, Japan over the last 0.46 My, and Condit and Connor (1996) observed continuous waxing and waning of certain areas across Springerville volcanic field, Arizona.

If it is thought that the locus of activity is moving in a specific direction, then a volcanic field boundary may need to be extended, or widened, in that direction to accommodate future eruptions. Or, if no young eruptions have been identified in a region of the volcanic field, this area could be removed from the current or future boundary region by shortening, or narrowing, the boundary in these regions. Recognising that tectonic stress regimes, upward mantle flow (driving partial melting), and plate motions all vary with time, a time-invariant volcanic field boundary is unlikely. However, for young fields, or if the influence of these on local volcanism is minor, then assuming temporal invariance may have little effect on subsequent analysis.

For Harrat Rahat, the two most recent eruptions occurred in the northern most part of the field. However, due to the high levels of uncertainty regarding the ages of the older eruptions, northward migration of the source region cannot be conclusively demonstrated. Thus, temporal invariance will be assumed. While this may be a conservative estimate for the boundary definition, it may lead to underestimation of the hazard, especially in the northern areas of the harrat.

Spӧrli and Eastwood (1997) found no shift in activity with time for the AVF, neither did more recent work by Bebbington and Cronin (2011) nor Bebbington (2013). Given that this field is very young (<0.25 Ma), temporal invariance is a justifiable assumption.

Unobserved (Hidden) vents

For the majority of volcanic fields, the eruption record is typically incomplete as eruptive centres can be both hidden (e.g., due to surface water, vegetation, or severe erosion), or buried (e.g., by sedimentation or subsequent eruptions). Geological, geophysical, seismic, and remote sensing studies can be used to reduce the number of unidentified vents/eruptions (Perry et al. 2005), but this is not always feasible or affordable. The main problem with an incomplete eruption record is that of underestimating temporal recurrence rates and spatial intensities (Runge et al. 2014). Further, hidden vents may lie significantly outside boundaries derived from observed vents, leading to underestimation of the field area. This argument can be used to justify the inclusion of a lava flow field’s extent within a volcanic field boundary, although surface topography and magma rheology are likely uncorrelated with the spatial location of source areas and magmas.

For Harrat Rahat, Runge et al. (2014) used the ages and extents of lava flows with an erosion function to estimate that two thirds of eruptive vents are missing from the record. This estimate correlated well with estimates of lava flow volumes (Camp and Roobol 1989) and volcanic pile thickness (Blank and Sadek 1983). The maximum thickness of the eruptive sequence across the harrat commonly reaches > 200 m and along the centre of the field can be up to 400 m (Blank and Sadek 1983). More than 50 % of observed cones and related vent structures are < 50 m in height, hence they are easily buried by lava in central regions of the field. Since total lava thickness dramatically reduces with distance from the central spine of the field (Fig. 1a), Blank and Sadek 1983), the likelihood of vent burial at its edges is low. Thus it is assumed that no unobserved vents are located outside of the region enclosed by those observed.

While the number of known eruption centres within the Auckland Volcanic Field has recently increased (e.g., Hayward et al. 2011), these were re-discovered and were present on an existing map of the field (Von Hochstetter 1864). The development of Auckland city almost entirely on top of this field suggests that there may be a few hidden eruption centres towards the centre of the field but the decrease in urbanisation towards the external vents means that it is assumed here that no unobserved vents are located outside of the region enclosed by those observed.

Anomalous vents

This section considers whether atypical volcanic vents within or near a volcanic field should be included or ignored. Anomalous eruptions might have distinct geochemistry, erupted volume, or age. There are also regions where polygenetic volcanoes coexist within, or overlap, the region of a volcanic field, for example, the Cascade Range (Guffanti and Weaver 1988), and the Eastern Volcanic Zone, Iceland (Jakobsson 1979).

Condit and Connor (1996) discounted three outlying vents from the Springerville volcanic field on the basis of previous cluster analysis work (Connor et al. 1992). Lutz and Gutmann (1995) discounted any distal vents deemed to be located at the edges of the Pinacate volcanic field when performing lineament analysis using a circular region. For island-based distributed volcanism (Jeju, Brenna et al. 2012; Ambae Island, Vanuatu, Németh and Cronin 2009), offshore eruptions are often discounted when forming a volcanic field boundary for spatial intensity estimation (Tenerife, Martí et al. 2009; El Hierro, Canaries, Melían et al. 2014), although they were included in work by Becerril et al. (2013). This may considerably underestimate the hazard, especially where magma-water interaction may cause the most explosive eruptions from such fields (Németh and Cronin 2009).

The most widely dispersed volcanic vents of Harrat Rahat occur at the most north-western, and south-eastern extents of the field (Fig. 1a). The vents in the north-west are among the most recent eruptions (641 AD), and those in the south-east among the oldest (though to be 8.5–2.5 Ma, Camp and Roobol 1989; Runge et al. 2014). Since there is no evidence other than their location to suggest these vents are anomalous, they must be included in the field boundary estimations.

Evidence of more complex eruption behaviour is also present in Harrat Rahat in the form of more evolved (trachytic) eruptions within a small region around 24°5′N. These eruptions have erupted through older more mafic volcanic structures. While anomalous in eruption style, these vents and their magma compositions are related to the basaltic lavas in this field, and are common to other volcanic fields on the Arabian Shield (Khaybar, Camp et al. 1991), as well as other similar volcanic fields in other parts of the world (Jeju, Brenna et al. 2012). The co-genetic relationship of these eruptions of differing composition leads to the decision to include all observed vents within the volcanic field boundary. However, more sophisticated modelling techniques that consider the behaviour as two separate processes could lead to multiple (likely overlapping) boundary definitions for the same field.

In the Auckland Volcanic Field, the most recent eruption and volcanic structure, Rangitoto, is anomalous in terms of size (Kereszturi et al. 2013), potential polygenetic eruption style, and composition (Needham et al. 2011). However, the most recent eruption may be the best indicator of future eruptions, and it is not uncommon for volcanic fields to contain large, shield-like volcanic structures (Mount Gambier, Bishop 2007; Jeju, Brenna et al. 2012), which may represent slightly different manifestations of volcanism from the same source area. Thus, in the AVF, all vents are included in the field boundary estimations.

Geological considerations

Knowledge of any field-limiting factors should be incorporated into the volcanic field boundary definition, such as crustal/fault structure, geological boundaries, or specific impenetrable lithologies. For example, the Uinkaret volcanic field, Colorado is bounded to the east and west by the Toroweap and Hurricane faults respectively (Karlstrom et al. 2007), and in the Garrotxa volcanic field, Spain, structural control on magma at depth was inferred and used to provide field-limiting constraints (Barde-Cabusson et al. 2014).

Another consideration may be the presence of ‘holes’ within a volcanic field representing localised regions of depleted mantle or crustal conditions that hinder magma reaching the surface. Stalled, or failed eruptions, where magma is intruded into the shallow crust may also provide additional boundary information. At Harrat Rahat, a 1999 AD earthquake swarm, was interpreted as a site of potential dyke intrusion (Lindsay and Moufti 2014), similar to that of the 2009 AD intrusion in nearby Harrat Lunayyir (Pallister et al. 2010). The location of these inferred intrusions could either be treated as sites where eruptions are now less likely, or conversely, where magma may be stalled, and instead erupt as a more evolved, and therefore potentially more explosive magma in the future. The incorporation of ‘holes’ in volcanic field boundary definitions would require a substantial change to how boundaries are construed beyond a simply connected polygon. However, with evidence for localised magma depletion, an initial attempt could begin with holes located at previous eruption sites.

Observed structures across and proximal to Harrat Rahat include suspected shear zones, terrane boundaries, and Precambrian sutures and faults (Camp and Roobol 1989). There are also volcanic fields nearby, which have been identified as separate systems based on the variation in vent clustering within and between regions, and relative geochemistry (Coleman et al. 1983). Although identified structures may control the orientations of vent lineaments (Connor 1990; Magill et al. 2005), little is known about their influence on magma pathways for Harrat Rahat, so no definitive field-limiting factors can be assumed.

The eastern boundary of the Auckland Volcanic Field was hypothesised at a contact with up-faulted greywacke basement rocks (Allen and Smith 1994). Consequently, Lindsay et al. (2010) extended their AVF boundary to the North where no outcropping of basement is observed. However, this geological bound was not considered here during boundary construction in the interest of a conservative approach, but this assumption should be revisited if more substantial support for this idea comes to light.

Field maturity

Fitting a boundary to the observed vents based on a minimum constraint, e.g., a convex hull, or minimum area polygon, is used to obtain a unique solution to a mathematical construct. However, this process means that those outermost vents are likely to lie on, or very close to, the volcanic field edge. This assumes that all future eruptions will occur in the currently confined area, an assumption that may not be valid for a very young field, or one with evidence of persistent vent migration, or where the most recent activity is focussed on a sub-area.

For the Auckland Volcanic Field, Spӧrli and Eastwood (1997) argued that vents at the outer limit of the current field extent varied substantially with age, potentially defining a constant magma source area. In many other settings however, this argument is hard to justify. Specifics of magma transport through the crust below volcanic fields are debated (Lister and Kerr 1991; Buck 2006; Tait and Taisne 2013), but generally it is thought to move via dyke propagation (Valentine and Hirano 2010). The geometry and orientation of dykes may be influenced by existing tectonic or geological boundary structures within the crust (Le Corvec et al. 2013b). The degree of lateral forcing of such structures is little known, however, if dyke propagation is forced too close to horizontal (e.g., via reaching a neutral buoyancy within the upper crust), a sill will form and magma is unlikely to reach the surface (Kavanagh et al. 2006).

As a conservative compromise to cover source uncertainty, a 5 km buffer zone was added to each fitted boundary. This value comes from Runge et al. (2014) and represents the > 95 % upper limit of dyke length from posterior distributions adapted to Harrat Rahat, with initial prior distributions taken from an eroded analogue volcanic field (San Rafael, Utah, Delaney and Gartner 1997), and expert elicitation. Although adapted specifically for Harrat Rahat, the same distance was applied to the Auckland Volcanic Field, while Le Corvec et al. (2013a) used ~1 km. In this work, Harrat Rahat was used to illustrate the sensitivity of spatial intensity results to the inclusion of a buffer. Thus, the buffer distance is applied to the AVF is purely to show the variation in scale (Fig. 2) and is not used in any sensitivity studies.

Fig. 2

Volcanic field boundary options, solid lines represent buffered boundaries, dashed lines represent minimum fitted shapes (see text for shape definitions), North is up. Harrat Rahat: convex hull (a), ellipse (b), rectangle (c), isocontour (d), anisocontour (e), red shaded area represents the extent of Al-Madinah city (Fig. 6a). Auckland Volcanic field: convex hull (f), ellipse (g), rectangle (h), isocontour (i), anisocontour (j)

For the minimum area ellipse and rectangle boundaries, a buffer zone is straightforward to implement, with 10 km added to both the major and minor axes of the shape. For a convex hull, each of the hull points is moved 5 km further away from the polygon’s centroid. For the contour boundaries, which already include a form of buffer zone, the boundary cut-off criteria was set to 10 % of the optimally determined value (‘no holes’ criterion) of the spatial density, thus increasing the bounds in all directions, but not necessarily by the same distance (Fig. 2).

Results – summary statistics

Five volcanic field boundaries for the two volcanic fields were constructed based on the above assumptions (Fig. 2). This section assesses the sensitivity of four commonly applied exploratory methods (analysing vent clustering and preferred orientations) to the choice of boundary definition. First-order statistics are based on a single measurable variable (e.g., the mean of a set of values). The methods of Clark and Evans (1954) and Hopkins and Skellam (1954) have been used in previous volcanic field work (Connor and Hill 1995; Martin et al. 2004) to determine the existence of clustering. Second-order statistics are based on variation, or distribution, of a measured value (e.g., the standard deviation). Here, the K-function (Ripley 1979) is applied to identify variations in clustering behaviour at different scales (Fig. 3) which has been used previously for the Auckland Volcanic Field (Magill et al. 2005).

Fig. 3

Visualisation of clustering statistics. Red triangles represent vents. Small blue circles represent randomly placed points. a Clark-Evans, b Hop-F, c K-function. See text for mathematical formula and explanation of individual statistics

The final exploratory method considered here is designed for lineament identification, the two-point azimuth (TPA) method of Lutz and Gutmann (1995) who applied it to the Pinacate volcanic field, Mexico. Connor (1990) also used TPA within and across inferred clusters of the TransMexican volcanic belt. All the exploratory methods were carried out before a buffer zone was applied. A buffer zone would cause an increase in area, a corresponding decrease in average vent density, and stronger clustering results.

Clark-Evans: 1st order randomness statistic

The Clark-Evans test compares the mean distance between neighbouring samples (${\overline{r}}_a$) to the typical distance of an equivalent random (homogeneous Poisson) distribution (${\overline{r}}_e$, Fig. 3a). The random distribution (${\overline{r}}_e$) requires an average vent density and therefore relies heavily on the estimated volcanic field area (A):

$$ {\overline{r}}_a = \frac{{\displaystyle {\sum}_{i=1}^N}{r}_i}{N},\kern0.5em {\overline{r}}_e=\frac{1}{2\sqrt{\rho }},\ where\kern0.5em \rho =\frac{N}{A} $$

where N is the number of vents, ρ is average vent density, r_i is the distance between the ith vent and nearest neighbour vent, and A is the area. ${\overline{r}}_e$ can also be estimated by placing the same number of points within the defined bounds via a Monte Carlo (MC) approach.

$$ R=\frac{{\overline{r}}_a}{{\overline{r}}_e} $$

(3)

A ratio (R) of the two of R = 1 suggests a random distribution, while R < 1 indicates clustering, and R > 1, regular spacing or dispersion (Clark and Evans 1954). The statistical significance can be assessed with the standard normal variate (c)

$$ c = \frac{{\overline{r}}_a-{\overline{r}}_e}{\sigma_{{\overline{r}}_e}},\kern0.5em where\kern0.5em {\sigma}_{{\overline{r}}_e} = \frac{0.26136}{\sqrt{N\rho }} $$

(4)

where ${\sigma}_{{\overline{r}}_e}$ is the standard deviation for the equivalent random distribution (${\overline{r}}_e$). The constant 0.26136 was derived from the probability of the number of points found within an area of specified size given a mean population density using a Poisson exponential function (derivation found in the appendix of Clark and Evans 1954). Where c lies outside of the range ± 2.58, the result can be assumed more than 99 % significant (Clark and Evans 1954), although Hamilton et al. (2010) suggest this critical value of c should be increased for analysis of low sample sizes.

All Clark-Evans results (R, Table 1) suggest significant clustering (R < 1) in Harrat Rahat. The elliptical volcanic field boundary suggests much stronger clustering than the other three due to the substantially larger area encompassed, much of which is empty. The Auckland Volcanic Field results (Table 2) are highly dependent on the choice of volcanic field boundary. The contour based boundaries suggest clustering, the ellipse and rectangular boundaries suggest random dispersion, and the convex hull suggests a more regularly spaced vent distribution. However, none of these results appear to be statistically significant (−2.58 < c < 2.58).

Table 1 Harrat Rahat boundary summary statistics. Clark-Evans (R) values are reported alongside the standard variate (c), Hop-F (A_g) values are reported with the standard error on the mean from 1000 simulations

Table 2 Auckland Volcanic Field boundary summary statistics. Clark-Evans (R) values are reported alongside the standard variate (c), Hop-F (A_g) values are reported with the standard error on the mean from 1000 simulations

Hop-F: 1st order randomness statistic

The Hop-F statistic compares the distance between vent I and its nearest neighbour vent (r_i), and the distance between a random point and the nearest vent (d_i) (Fig. 3b).

$$ {A}_g = \frac{{\displaystyle {\sum}_{i=1}^N}{d}_i^2}{{\displaystyle {\sum}_{i=1}^N}{r}_i^2} $$

(5)

where N is the number of vents. The resulting coefficient of aggregation (A_g) indicates randomness if A_g = 1, clustering if A_g > 1, and regular spacing or dispersion if A_g < 1 (Hopkins and Skellam 1954). The statistic A_g is most often reported with the standard error on the mean using the results of at least 100 simulations (Cressie 1991). It is this simulation step that requires definition of the specific field boundary.

The Hop-F results from 1000 simulations (A_g, Table 1) for the ellipse, rectangle, and convex hull suggest random dispersal of vents across Harrat Rahat with the standard error bounds all including unity (random). The isocontour boundary suggests slightly dispersed vent locations (A_g < 1). For the Auckland Volcanic Field, results vary dramatically as with the Clark-Evans results, however, similar to the Harrat Rahat, almost all results include the null hypothesis of randomness. The exception is again that of the isocontour, which suggests dispersed vent locations, and the ellipse boundary (although the confidence bound includes one) that more strongly suggests clustering than any of the other bounds.

While the Clark-Evans statistic tests the degree of departure from a Poisson distribution (Fig. 3a), the Hop-F statistic suggests clustering (Fig. 3b). These results suggest that the Clark-Evans statistic is more sensitive to departures from randomness other than clustering, such as lineaments and anisotropy. The presence of multiple-vent events (Runge et al. 2014) may have more influence on the Hop-F statistic when using a more closely fitting boundary within which random points can fall. Since contour boundaries tailor the boundary to the vent locations, they effectively discard empty space furthest from vents and the false-positive clustering results are avoided. However, without knowledge of the actual volcanic field boundary, it is difficult to determine if this makes the contour results more reliable for awkward volcanic field shapes, such as that of Harrat Rahat, or instead, overly biased towards observed eruption locations.

K-Function: 2nd order randomness statistic

This second order test compares the number of vents within distance d of another vent, with the number expected for a random distribution with the same density (Fig. 3c)

$$ \widehat{K}(d) = \frac{A}{N^2}{\displaystyle {\sum}_{j=1}^N{\displaystyle {\sum}_{i\ne j}^N\frac{H\left(d-{d}_{i,j}\right)}{w_{i,j}}}} $$

(6)

where N is the number of vents, A is the area of the field (requires a boundary), d_i,j is the distance between ith and jth vents, H is the Heaviside function, and w_i,j is an edge correction (Martinez and Martinez 2001). This statistic is most easily interpreted as a graph of the L-function

$$ \widehat{L}(d) = \sqrt{\frac{\widehat{K}(d)}{\pi }}-d $$

(7)

where $\widehat{L}(d)\kern0.5em =\kern0.5em 0$ suggests randomness, $\widehat{L}(d)\kern0.5em >\kern0.5em 0$implies clustering, and $\widehat{L}(d)\kern0.5em <\kern0.5em 0$, regular spacing (Ripley 1979). Significance levels are determined via random Poisson (MC) simulations run over the defined area.

K-function (or L-function) results are often diluted or smoothed if the method is applied over larger areas containing highly varying vent densities. The strong spine-type structure of Harrat Rahat (Fig. 1a) thus produces unusual results when this method is applied across the whole field and, in practice, much smaller regions should be assessed to look for specific clustering distances. All the volcanic field boundaries fitted to Harrat Rahat show clustering ($\widehat{L}(d)\kern0.5em >\kern0.5em 0$) at all distances (Fig. 4a), with similar shapes following similar trends (e.g., both contour boundaries follow similar lines, as do the convex hull and rectangular bounds).

Fig. 4

L-function ($\widehat{L}(d)$) statistic for Harrat Rahat (HR) and Auckland Volcanic field (AVF) vent locations within various field boundaries. Thin lines represent upper and lower 95 % significance levels from 1000 random (Poisson) simulations. Thick black lines represent vent location based values. a HR All boundaries, b AVF convex hull, c AVF ellipse, d AVF rectangle, e AVF isocontour, f AVF anisocontour. See Fig. 2 for boundary illustrations, and text for L-function definition

For the Auckland Volcanic Field, results are more informative for actual clustering distances due to the more evenly dispersed eruption locations, and lower vent numbers (Fig. 4). All results show the minimum vent distance at ~ 300 m (first trough), and suggest clustering at distances > 11 km (vent results lie above the upper significance bound), which likely reflects the maximum/regular vent-spacing, and consequently, field dimensions. The convex hull, ellipse, and rectangular bounds suggest clustering between ~ 0.5 and 2.5 km (above upper bound), and then maximum/regular spacing between 3 and 6 km (below lower bound, Fig. 4c and d) (3 and 8 km for the convex hull, Fig. 4b), and then clustering at larger distances. The two contour boundaries however, suggest clustering at all distances above ~ 0.5 km. Previous results from application of the K-function (Ripley 1979) to 49 vents of the Auckland Volcanic Field by Magill et al. (2005) found clustering between 0.9 and 1.6 km, and then maximum spacing at 4.6 km, which closely resembles the convex hull results presented here. The significant variation in results between boundaries is important. Although the boundary shapes of the Auckland Volcanic Field are broadly similar, the extra space incorporated around the external vents in the contour methods dilutes the second-order statistics by effectively forcing a cluster in the centre of the region (Fig. 2 and Fig. 4e-f).

TPA: lineament identification method

Lutz and Gutmann (1995)’s method for lineament identification (Two Point Azimuth: TPA) is based on the azimuth between each sample and every other sample and provides a distribution of the frequency of alignment angles, usually presented within groups of 10° intervals (Connor 1990; Wadge and Cross 1988). To accommodate any non-circular field shape, the results are then compared to Monte Carlo simulations with the same number of samples placed within the defined area, thus the results are dependent on the specific field boundary definition.

Harrat Rahat TPA results for the five volcanic field boundaries show approximately the same results, but with, in some cases, substantially varying proportions (Fig. 5). For example, the convex hull results show that vent alignments are only slightly more likely to be between 160° and 170° than randomly placed points within the defined field area.

Fig. 5

TPA frequency histograms for Harrat Rahat (HR) and Auckland Volcanic Field (AVF) vent locations within various field boundaries. Grey lines represent the 95 % significance levels from 1000 random (Poisson) simulations. Black lines represent vent location based values. a HR convex hull, b HR ellipse, c HR rectangle, d HR isocontour, e HR anisocontour, f AVF convex hull, g AVF ellipse, h AVF rectangle, i AVF isocontour, j AVF anisocontour

In contrast, the elliptical boundary results show ~ 50 % more alignments than would be expected for random points within the elliptical field (Fig. 5b, the difference in vent-vent alignment frequency for the thick black line representing vents, and the grey lines representing significance bounds). Had only a convex hull boundary been used, these apparent NNW/SSE vent alignments may have been discarded as an artefact of overall volcanic field shape. However, the fact that there are consistently more vents than would be expected aligned between 160° and 170° may suggest similar structural controls on at both large (the region of volcanism) and small (vent lineaments/structures) scales. The only variation is that of the anisocontour boundary which does not pick up the peak that the other four bounds do (160° and 170°). This is a result of the base kernel (as defined by H, Eqn. 2) being pre-aligned (via SAMSE) to the overall dominant alignment direction, i.e., this alignment direction of vents is not anomalous within this boundary definition, as all the simulations are based on, and therefore reproduce, it.

Across the Auckland Volcanic Field, all boundaries (except the isocontour) suggest a slight tendency for alignments between 110° and 120° (Fig. 5). This matches the work of Von Veh and Németh (2009) who used the Hough Transform (Wadge and Cross 1988) and an elliptical boundary. The rectangular bounds suggest a negative preference for alignments between 40° and 50°. Although several of the other results get close, none breach the 95 % significance levels. Unlike the previous clustering results (Table 2, Fig. 4), the TPA results are insensitive to the volcanic field boundary of the Auckland Volcanic Field because this analysis is based on field shape (relative proportions) rather than absolute dimensions.

Results – spatial intensity estimates

For exploratory methods, a buffer is superfluous as it is the existing locations and spacings of eruptions that are assessed. However, for forecasts, field maturity and potential eruption migration need to be considered, with a buffer zone included in the definition of a boundary as a conservative measure. However, the addition of a buffer zone is non-trivial, especially where specific spatial densities are important as excessively increasing the volcanic field area will result in an underestimate of average vent-density and consequently spatial intensity. This is illustrated here using the city of Al-Madinah, located immediately to the north-west of Harrat Rahat (Fig. 1a) which is spread over an area of ~30 km² (Fig. 6a), with a population of 1.5 million.

Fig. 6

Harrat Rahat boundary variations overlain by the city limits of Al-Madinah (red). Thin black lines represent fitted shapes, thick black lines represent boundaries with buffer, North is up. a 1256 AD lava flow extent shown in red, Al-Madinah city boundary shown in blue. b convex hull, c ellipse, d rectangle, e isocontour, f anisocontour

A large number of methods can be employed to estimate spatial density within volcanic fields (Marzocchi and Bebbington 2012; Bebbington 2013). Here, a uniform spatial intensity is assumed (the simplest of the available methods), whereby all regions within the volcanic field boundary are assumed equally likely to host future eruptions. This allows the sensitivity of the spatial intensity results to boundary choice to be assessed without assumptions inherent to more complex methods confounding the comparison. For a uniform spatial density model, given an eruption occurs (probability across the whole region must sum to one), the spatial density at a point (γ(x)) is equal to 1/A where A is the total field area. The probability of an eruption occurring within the area of interest then becomes the ratio between the proportion of the area of interest (in this case, the city of Al-Madinah) within the volcanic field boundary and the total area of the field.

Four buffer variations were assessed for each boundary type: (1) a fitted boundary (no buffer), (2) an equal intensity buffered boundary, and (3),(4) two inner-outer boundary approaches. For these latter two, the uniform spatial density estimate is imposed across the region, but with the values in the outer boundary set as a proportion (p ≤ 1) (50 and 10 % here) of that of the inner boundary.

$$ \upgamma \left(\mathrm{x}\right)=\Big\{\begin{array}{c}\frac{1}{\eta_o{A}_1}x\kern0.5em \in \kern0.75em {A}_1(inner)\hfill \\ {}\kern1.5em \frac{p}{\eta_o{A}_1}x\kern0.5em \in \kern0.5em {A}_2\&x\kern0.5em \notin \kern0.5em {A}_1(outer)\kern1em \end{array},\kern1em {\eta}_o={\left[p\left(\frac{A_2}{A_1}-1\right)+1\right]}^{-1} $$

(8)

where A₁ is the area within the inner boundary, A₂ is the area within the outer boundary (i.e., A₁ ⊆ A₂, and A₂ ≥ A₁). This is a generalised form of that suggested in Bebbington (2013) in which the inner area was assumed to have twice the probability of the outer.

As the city of Al-Madinah borders the current flow field extents, the effect of volcanic field boundary definition on the area susceptible to future eruptions is substantial (Fig. 6). For the fitted boundaries (without a buffer) the anisocontour boundary already includes the whole of Al-Madinah, and both the isocontour and elliptical bounds incorporate almost the entire extent of Al-Madinah within the area deemed susceptible to a future eruption. Thus, the addition of the buffer zones results in a decrease in eruption probability within the city limits as the spatial density is spread over a larger region (Table 3). Conversely, the convex hull and rectangular boundaries enclose less than half of Al-Madinah without a buffer zone, and thus when the buffer is added, the proportion of the city within the boundary increases sufficiently to also increase the eruption probability, even though the spatial intensity at any specific point is reduced (Table 3).

Table 3 Eruption probability estimates for an eruption within the area of interest, the city of Al-Madinah. Estimates are based on a uniform distribution of eruption probability, given that an eruption occurs within the region defined for Harrat Rahat

The estimate of the intensity within Al-Madinah therefore is a balance between the extent of the city enclosed by the bound, and the spread of the intensity across the whole field. Comparing across boundaries, the isocontour boundary with no buffer resulted in the highest probability of eruption within the area of Al-Madinah city, with values more than twice that of the corresponding unbuffered convex hull and rectangular boundaries (Table 3). The probability of an eruption within the city limits (given that an eruption occurs) was assessed as 0.0206, or a 2.1 % chance. OF the five boundary shapes tested, the convex hull results were most sensitive to the addition of a buffer zone, varying up to ~160 % (Table 3), predominantly due to the change in proportion of the area of interest falling within a susceptible region (Fig. 6b). This is closely followed by the rectangular volcanic field boundary with a variation of ~140 % for the same reason (Fig. 6d).

Discussion

The results above show the significant effect that the definition of a volcanic field boundary has on even the simplest of exploratory methods (Tables 1 and 2). The effect on spatial density estimation is also substantial, with the estimate of eruption probability within Al-Madinah varying significantly even between the few boundaries assessed here (Table 3). The use of a constant (fixed) boundary is also questionable, especially in those cases (such as Harrat Rahat), where a possible shift in eruption locations with time is noted. In these cases, a temporally varying boundary may be more appropriate by allowing the boundary definition to expand with time, or even move with time for those fields where regions of older activity are considered no longer active – a limited-memory boundary could be imposed which is fitted only to the eruptions from the most recent x years.

For forecasting purposes, a buffer zone is recommended as a conservative approach. It allows some accommodation of spatio-temporal migration or focussing (e.g., Harrat Rahat), and on-going field development (e.g., Auckland Volcanic Field). However, the size of the buffer zone was not investigated here and while the 5 km buffer zone imposed for Harrat Rahat is a relatively small increase in total area, this value has a substantial effect on the area designated to the Auckland Volcanic Field (Fig. 2). A buffer zone distance based on the dimensions of a volcanic field may be more appropriate (e.g., the ~1 km buffer zone imposed by Le Corvec et al. 2013a).

The 5 km buffer zone for Harrat Rahat was taken from previous estimates of dyke length (Runge et al. 2014), however it is likely that the distribution of dyke lengths, and subsequently the distance along which magma may erupt from, may vary with field thus, at the very least, the tectonic regime should be considered (Pinel and Jaupart 2004). The presence of fissure eruptions (as noted in Harrat Rahat) may also justify an increase in buffer zone distance, and, if multiple fissures show similar alignments, anisotropic buffer zones may be more appropriate. Evidence of waning may also suggest a decrease (or even exclusion) of a buffer zone as the future eruptions may be smaller, and from shorter subsurface structures (Valentine and Perry 2006).

Smaller scale structures should also be considered during volcanic field definition. For example, in Harrat Rahat, where there is strong anisotropic structure (Runge et al. 2014), or dominant alignment directions, the anisotropic kernel based boundary added a large amount of empty space in the directions determined by the SAMSE optimised base kernel, which in turn is strongly influenced by regional vent alignments. The isotropic kernel based contour resulted in a much more irregular volcanic field boundary than that of the anisotropic kernel (Fig. 1), mainly due to the smaller optimum bandwidth (also reflected in the TPA results for Harrat Rahat, Fig. 5). This is consistent with previous work that suggests vent-specific locations are influenced by smaller-scale local, rather than regional, structures (Connor et al. 1992).

For Harrat Rahat, none of the boundary methods would have included the 641 AD eruption (the most north-westerly vent area) should this analysis have been carried out in 640 AD. This suggests an evolving field boundary and a buffer zone should be imposed for Harrat Rahat, and other volcanic fields exhibiting strong spatio-temporal migration.

For the Auckland Volcanic Field, with the exception of the rectangular boundary, the range of fitting processes resulted in similar shapes. However, had the bounds been fitted prior to the most recent eruption (Rangitoto), only the elliptical and rectangular bounds would have included it, reflecting the unusual number of vents within a short distance of the elliptical boundary (Spӧrli and Eastwood 1997; Bebbington 2015). While the field is still relatively young (Lindsay et al. 2010), previous work has suggested that the spatial distribution of eruptions within the AVF has been relatively stable since the formation of the 15th eruption (Bebbington 2013; Le Corvec et al. 2013a). This may suggest the addition of a buffer zone is unnecessary. However, as the most recent eruption was anomalous both in terms of polygenetic behaviour (Needham et al. 2011) and large erupted volume (Kereszturi et al. 2013), a buffer zone is advised to conservatively accommodate, at least in part, these potential changes in eruption characteristics.

Conclusions

Forecasts of long-term volcanic activity should be accompanied by an estimate of their precision and uncertainties imposed by the assumptions made. However, hazard analyses are required for incompletely understood volcanic regions threatened by future eruptions. This requires parameterisation of the unknown. This study showed the substantial effect that the definition of a volcanic field boundary may have on even the simplest exploratory methods for hazard forecasting.

As research on volcanic fields and geophysical methods progress, boundary fitting may become better informed by more extensive knowledge of specific field constraints, however, a region definition will always be vital for probabilistic approaches. When carrying out any analysis on volcanic fields, it is vital to specify how the volcanic field boundary has been defined, assess the sensitivity of boundary choice, and to carry these assumptions and related uncertainties through to estimates of future activity and related hazard analyses.

For Harrat Rahat, sensitivity of results to boundary definition is substantial and a buffer zone is highly recommended due to the spatio-temporal behaviour of the field. For the Auckland Volcanic Field, the range of options resulted in similar shapes, nevertheless, some of the statistical tests still showed substantial variation in results.

This work highlights the fact that when carrying out any hazard analysis on volcanic fields, it is vital to specify how the volcanic field boundary has been defined, assess the sensitivity of boundary choice, and to carry these assumptions and related uncertainties through to estimates of future activity and hazard analyses.