 Research
 Open Access
Sensitivity to volcanic field boundary
 Melody G. Runge^{1}Email author,
 Mark S. Bebbington^{2},
 Shane J. Cronin^{1},
 Jan M. Lindsay^{1} and
 Mohammed Rashad Moufti^{3}
https://doi.org/10.1186/s136170150040z
© Runge et al. 2015
 Received: 20 April 2015
 Accepted: 7 December 2015
 Published: 15 December 2015
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 spatiotemporal point processes with model selection and development based on exploratory analyses of previous eruption data. For identifiability reasons, spatiotemporal 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 subsurface 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.
Keywords
 Probabilistic
 Hazard forecasting
 Harrat Rahat
 Auckland volcanic field
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 spatiotemporal 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 spatiotemporal 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 spatiotemporal evolution of the field, and then extrapolated into the future to obtain hazard forecasts. This almost invariably uses spatial or spatiotemporal 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 groundwater 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 spatiotemporal models (Connor and Hill 1995; Cronin et al. 2001), and for assessment of all spatial and spatiotemporal 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 presentday view of the location of past volcanic material from a resourcemapping 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 hotspot upwelling), boundary parameters remain unknown, thus a wellinformed 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 spatiotemporal 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 biasvariance tradeoff. 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 ventdensity 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 ventdensity 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
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 ^{40}Ar/^{39}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 AlMadinah (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^{3} (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 giftwrapping (Chand and Kapur 1970), or divideandconquer (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 semicircular in planview, 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 nonlinear 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 subregion 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 riftsubparallel structure known as the MekkahMadinah 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 bestfit is then the rectangle corresponding to the minimum area.
Isotropic and anisotropic kernels (Vent Density Contours)
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‖^{2}: 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 (‘Squaredasymptoticmean 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 userdefined cutoff 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

(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 timeinvariant 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 rediscovered 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 islandbased 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 magmawater 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 northwestern, and southeastern extents of the field (Fig. 1a). The vents in the northwest are among the most recent eruptions (641 AD), and those in the southeast 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 cogenetic 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, shieldlike 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 fieldlimiting 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 fieldlimiting constraints (BardeCabusson 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 fieldlimiting factors can be assumed.
The eastern boundary of the Auckland Volcanic Field was hypothesised at a contact with upfaulted 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 subarea.
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).
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 cutoff 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
The final exploratory method considered here is designed for lineament identification, the twopoint 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.
ClarkEvans: 1st order randomness statistic
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.
Harrat Rahat boundary summary statistics. ClarkEvans (R) values are reported alongside the standard variate (c), HopF (A_{g}) values are reported with the standard error on the mean from 1000 simulations
Boundary  Area (km^{2})  Vent density (km^{−2})  R (ClarkEvans)  A_{g} (HopF) 

Convex hull  12,535  0.077  0.652, c = −20.77  1.014 ± 0.047 
Ellipse  22,430  0.043  0.487, c = −30.58  1.006 ± 0.043 
Rectangle  14,717  0.066  0.601, c = −23.77  0.978 ± 0.048 
Isocontour  13,496  0.072  0.628, c = −22.18  0.914 ± 0.038 
Anisocontour  15,292  0.063  0.5899, c = −24.45  0.956 ± 0.047 
Auckland Volcanic Field boundary summary statistics. ClarkEvans (R) values are reported alongside the standard variate (c), HopF (A_{g}) values are reported with the standard error on the mean from 1000 simulations
Boundary  Area (km^{2})  Vent density (km^{−2})  R (ClarkEvans)  A_{g} (HopF) 

Convex hull  324.5  0.1527  1.134, c = 1.829  0.887 ± 0.160 
Ellipse  385.7  0.129  1.04, c = 0.547  1.191 ± 0.242 
Rectangle  400.2  0.127  1.021, c = 0.287  1.0582 ± 0.226 
Isocontour  585.7  0.087  0.844, c = −2.132  0.793 ± 0.141 
Anisocontour  579.4  0.088  0.849, c = −2.069  0.896 ± 0.173 
HopF: 1st order randomness statistic
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 HopF 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 ClarkEvans 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 ClarkEvans statistic tests the degree of departure from a Poisson distribution (Fig. 3a), the HopF statistic suggests clustering (Fig. 3b). These results suggest that the ClarkEvans statistic is more sensitive to departures from randomness other than clustering, such as lineaments and anisotropy. The presence of multiplevent events (Runge et al. 2014) may have more influence on the HopF 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 falsepositive 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.
KFunction: 2nd order randomness statistic
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.
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 ventspacing, 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 Kfunction (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 secondorder statistics by effectively forcing a cluster in the centre of the region (Fig. 2 and Fig. 4ef).
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 noncircular 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.
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 ventvent 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 prealigned (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
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 AlMadinah) within the volcanic field boundary and the total area of the field.
where A_{1} is the area within the inner boundary, A_{2} is the area within the outer boundary (i.e., A_{1} ⊆ A_{2}, and A_{2} ≥ A_{1}). 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.
Eruption probability estimates for an eruption within the area of interest, the city of AlMadinah. Estimates are based on a uniform distribution of eruption probability, given that an eruption occurs within the region defined for Harrat Rahat
Boundary  No buffer  Equal intensity buffer  50 % intensity buffer  10 % intensity buffer 

Convex hull  0.0092  0.0147  0.01227  0.0099 
Ellipse  0.0132  0.0119  0.0125  0.0131 
Rectangle  0.0100  0.0146  0.0125  0.0106 
Isocontour  0.0206  0.0190  0.0198  0.0204 
Anisocontour  0.0200  0.0184  0.0192  0.0199 
The estimate of the intensity within AlMadinah 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 AlMadinah 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 AlMadinah 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 limitedmemory 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 spatiotemporal migration or focussing (e.g., Harrat Rahat), and ongoing 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 ventspecific locations are influenced by smallerscale 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 northwesterly 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 spatiotemporal 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 longterm 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 spatiotemporal 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.
Declarations
Acknowledgements
This work is part of the VORiSA (Volcanic Risk in Saudi Arabia) collaborative project between King Abdulaziz University and The University of Auckland. Valuable comments were also provided by Laura Sandri and an anonymous reviewer.
Open AccessThis 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.
Authors’ Affiliations
References
 Allen S, Smith I. Eruption styles and volcanic hazard in the Auckland Volcanic Field, New Zealand. Geosci Rep Shizuoka Univ. 1994;20:5–14.Google Scholar
 Angell I, Barber J. An algorithm for fitting circles and ellipses to megalithic stone rings. Sci Archaeol. 1977;20:11–6.Google Scholar
 BardeCabusson S, Gottsmann J, Martí J, Bolós X, Camacho A, Geyer A, et al. Structural control of monogenetic volcanism in the Garrotxa volcanic field (Northeastern Spain) from gravity and selfpotential measurements. Bull Volcanol. 2014;76(1):1–13.View ArticleGoogle Scholar
 Bebbington M. Assessing spatiotemporal eruption forecasts in a monogenetic volcanic field. J Volcanol Geotherm Res. 2013;252:14–28.View ArticleGoogle Scholar
 Bebbington MS. Spatiovolumetric hazard estimation in the Auckland Volcanic Field. Bull Volcanol. 2015;77:39.View ArticleGoogle Scholar
 Bebbington M, Cronin S. Spatiotemporal hazard estimation in the Auckland Volcanic Field, New Zealand, with a new eventorder model. Bull Volcanol. 2011;73:55–72.View ArticleGoogle Scholar
 Becerril L, Cappello A, Galindo I, Neri M, Del Negro C. Spatial probability distribution of future volcanic eruptions at El Hierro Island (Canary Islands, Spain). J Volcanol Geotherm Res. 2013;257:21–30.View ArticleGoogle Scholar
 Bishop M. Point pattern analysis of eruption points for the Mount Gambier volcanic subprovince: a quantitative geographical approach to the understanding of volcano distribution. Area. 2007;39(2):230–41.View ArticleGoogle Scholar
 Blank H, Sadek H. Spectral analysis of the 1976 aeromagnetic survey of Harrat Rahat, Kingdom of Saudi Arabia. 1983. Tech Report USGSOF0367, USGS.Google Scholar
 Brenna M, Cronin S, Smith I, Sohn Y, Maas R. Spatiotemporal evolution of a dispersed magmatic system and its implications for volcano growth, Jeju Island Volcanic Field, Korea. Lithos. 2012;148:337–52.View ArticleGoogle Scholar
 Brown G. Eastern margin of the Red Sea and the coastal structures in Saudi Arabia. Phil Trans R Soc Lond Ser A Math Phys Sci. 1970;267:75–87. doi:10.1098/rsta.1970.0024.View ArticleGoogle Scholar
 Buck W. The role of magma in the development of the AfroArabian Rift System. Geol Soc Lond Spec Pub. 2006;259(1):43–54.View ArticleGoogle Scholar
 Cabrera J, Meer P. Unbiased estimation of ellipses by bootstrapping, pattern analysis and machine intelligence. IEEE Trans. 1996;18(7):752–6.Google Scholar
 Camp VE, Roobol MJ. The Arabian continental alkali basalt province: Part I. Evolution of Harrat Rahat, Kingdom of Saudi Arabia. Geol Soc Am Bull. 1989;101(1):71–95.View ArticleGoogle Scholar
 Camp VE, Roobol MJ, Hooper PR. The Arabian continental alkali basalt province: Part II. Evolution of Harrats Khaybar, Ithnayn, and Kura, Kingdom of Saudi Arabia. Geol Soc Am Bull. 1991;103(3):363–91.View ArticleGoogle Scholar
 Chand D, Kapur S. An algorithm for convex polytypes. J ACM. 1970;17(1):78–86.View ArticleGoogle Scholar
 Clark P, Evans F. Distance to nearest neighbour as a measure of spatial relationships in populations. Ecology. 1954;35:s–53.View ArticleGoogle Scholar
 Coleman RG, Gregory RT, Brown GF. Cenozoic volcanic rocks of Saudi Arabia. 1983. Tech Rep Open File Report USGSOF93, Saudi Arabian Deputy Minister of Mineral Resources.Google Scholar
 Condit C, Connor C. Recurrence rates of volcanism in basaltic volcanic fields, an example from the Springerville volcanic field, Arizona. Geol Soc Am Bull. 1996;108:1225–41.View ArticleGoogle Scholar
 Connor C. Cinder cone clustering in the TransMexican volcanic belt: implications for structural and petrologic models. J Geophys Res Solud Earth (1978–2012). 1990;95(B12):19,395–405.View ArticleGoogle Scholar
 Connor C, Connor L. Estimating spatial density with kernel methods. Volcanic and tectonic hazard assessment for nuclear facilities. Cambridge: Cambridge University Press; 2009. p. 346–68.View ArticleGoogle Scholar
 Connor C, Conway F. Basaltic volcanic fields. In: Encyclopedia of volcanoes. New York: Academic; 2000. p. 331–43.Google Scholar
 Connor C, Hill B. Three nonhomogeneous Poisson models for the probability of basaltic volcanism: application to the Yucca Mountain region, Nevada. J Geophys Res Solid Earth (1978–2012). 1995;100(B6):10,107–25.View ArticleGoogle Scholar
 Connor C, Condit C, Crumpler L, Aubele J. Evidence of regional structural controls on vent distribution: Springerville volcanic field, Arizona. J Geophys Res Solid Earth (1978–2012). 1992;97(B9):12,349–59.View ArticleGoogle Scholar
 Cressie N. Statistics for spatial data, Wiley Series in Probability and Mathematical Statistics: Applied Probability and Statistics. New York: Wiley; 1991.Google Scholar
 Cronin S, Bebbington M, Lai C. A probabilistic assessment of eruption recurrence on Taveuni volcano, Fiji. Bull Volcanol. 2001;63(4):274–88.View ArticleGoogle Scholar
 Delaney P, Gartner A. Physical processes of shallow mafic dike emplacement near the San Rafael Swell, Utah. Geol Soc Am Bull. 1997;109(9):1177–92.View ArticleGoogle Scholar
 Duong T. ks: Kernel density estimation and kernel discriminant analysis for multivariate data in R. J Stat Softw. 2007;21(7):1–16.View ArticleGoogle Scholar
 Duong T, Hazelton ML. Plugin bandwidth selectors for bivariate kernel density estimation. J Nonparametric Stat. 2003;15:17–30.View ArticleGoogle Scholar
 El Difrawy M, Runge M, Moufti MR, Cronin S, Bebbington M. A first hazard analysis of the Quaternary Harrat AlMadinah volcanic field, Saudi Arabia. J Volcanol Geotherm Res. 2013;267:39–46.View ArticleGoogle Scholar
 Fedotov S. Magma rates in feeding conduits of different volcanic centres. J Volcanol Geotherm Res. 1981;9(4):379–94.View ArticleGoogle Scholar
 Freeman H, Shapira R. Determining the minimumarea encasing rectangle for an arbitrary closed curve. Commun ACM. 1975;18(7):409–13.View ArticleGoogle Scholar
 Germa A, Connor L, CanonTapia E, Le Corvec N. Tectonic and magmatic controls on the location of postsubduction monogenetic volcanoes in Baja California, Mexico, revealed through spatial analysis of eruptive vents. Bull Volcanol. 2013;75(12):1–14.View ArticleGoogle Scholar
 Guffanti M, Weaver C. Distribution of late Cenozoic volcanic vents in the Cascade Range: volcanic arc segmentation and regional tectonic considerations. J Geophys Res Solid Earth (1978–2012). 1988;93(B6):6,513–29.View ArticleGoogle Scholar
 Hamilton CW, Fagents SA, Thordarson T. Explosive lavawater interactions ii: selforganization processes among volcanic rootless eruption sites in the 1783–1784 laki lava flow, Iceland. Bull Volcanol. 2010;72(4):469–85.View ArticleGoogle Scholar
 Hayward BW, Murdoch G, Maitland G. Volcanoes of Auckland: the essential guide. N Z Sci Rev. 2011;68:123.Google Scholar
 Hopkins B, Skellam J. A new method for determining the type of distribution of plant individuals. Ann Bot. 1954;18(2):213–27.Google Scholar
 Jakobsson S. Petrology of recent basalts of the Eastern Volcanic Zone, Iceland. Acta Nat Islandica. 1979;26:1–103.Google Scholar
 Kanatani K. Statistical optimization for geometric computation: theory and practice. 2005. Courier Dover Publications.Google Scholar
 Karlstrom K, Crow R, Peters L, McIntosh W, Raucci J, Crossey L, et al. 40Ar/39Ar and field studies of Quaternary basalts in Grand Canyon and model for carving Grand Canyon: Quantifying the interaction of river incision and normal faulting across the western edge of the Colorado Plateau. Geol Soc Am Bull. 2007;119(11–12):1,283–1312.View ArticleGoogle Scholar
 Kavanagh J, Menand T, Sparks R. An experimental of sill formation and propogation in layered elastic media. Earth Planet Sci Lett. 2006;245(3):799–813.View ArticleGoogle Scholar
 Kear D. Volcanic alignments north and west of New Zealand’s central volcanic region. N Z J Geol Geophys. 1964;7(1):24–44.View ArticleGoogle Scholar
 Kereszturi G, Németh K, Cronin S, AgustínFlores J, Smith I, Lindsay J. A model for calculating eruptive volumes for monogenetic volcanoesimplication for the Quaternary Auckland Volcanic Field, New Zealand. J Volcanol Geotherm Res. 2013;266:16–33.View ArticleGoogle Scholar
 Khachiyan L. Rounding of polytopes in the real number model of computation. Math Oper Res. 1996;21(2):307–20.View ArticleGoogle Scholar
 Kiyosugi K, Connor C, Zhao D, Connor L, Tanaka K. Relationships between volcano distribution, crustal structure, and Pwave tomography: an example from the Abu Monogenetic Volcano Group, SW Japan. Bull Volcanol. 2010;72(3):331–40.View ArticleGoogle Scholar
 Le Corvec N, Bebbington M, Lindsay J, McGee L. Age, distance, and geochemical evolution within a monogenetic volcanic field: analyzing patterns in the Auckland Volcanic Field eruption sequence. Geochem Geophys Geosyst. 2013a;14(9):3,648–65.View ArticleGoogle Scholar
 Le Corvec N, Menand T, Lindsay J. Interaction of ascending magma with preexisting crustal fractures in monogenetic basaltic volcanism: an experimental approach. J Geophys Res Solid Earth. 2013b;118(3):968–84.View ArticleGoogle Scholar
 Le Corvec N, Spӧrli K, Rowland J, Lindsay J. Spatial distribution and alignments of volcanic centers: clues to the formation of mono genetic volcanic fields. Earth Sci Rev. 2013c;124:96–114.View ArticleGoogle Scholar
 Lindsay J, Moufti MR. Assessing volcanic risk in Saudi Arabia. Eos Trans Am Geophys Union. 2014;95(31):277–8.View ArticleGoogle Scholar
 Lindsay J, Marzocchi W, Jolly G, Constantinescu R, Selva J, Sandri L. Towards realtime eruption forecasting in the Auckland Volcanic Field: application of BETEF during the New Zealand national disaster exercise ‘Ruaumoko’. Bull Volcanol. 2010;72(2):185–204.View ArticleGoogle Scholar
 Lindsay J, Leonard G, Smid E, Hayward B. Age of the Auckland Volcanic Field: a review of existing data. N Z J Geol Geophys. 2011;54(4):379–401.View ArticleGoogle Scholar
 Lister J, Kerr R. Fluidmechanical models of crack propagation and their application to magma transport in dykes. J Geophys Res Solid Earth (1978–2012). 1991;96(B6):10,049–77.View ArticleGoogle Scholar
 Lutz T, Gutmann J. An improved method for determining and characterizing alignments of pointlike features and its implications for the Pinacate volcanic field, Sonora, Mexico. J Geophys Res Solid Earth (1978–2012). 1995;100(B9):17,659–70.View ArticleGoogle Scholar
 Magill C, McAneney K, Smith I. Probabilistic assessment of vent locations for the next Auckland Volcanic Field event. Math Geol. 2005;37(3):227–42.View ArticleGoogle Scholar
 Martí J, Ortiz R, Gottsmann J, García A, De La CruzReyna S. Characterising unrest during the reawakening of the central volcanic complex on Tenerife, Canary Islands, 2004–2005, and implications for assessing hazards and risk mitigation. J Volcanol Geotherm Res. 2009;182(1):23–33.View ArticleGoogle Scholar
 Martin A, Umeda K, Connor C, Weller J, Zhao D, Takahashi M. Modeling longterm volcanic hazards through Bayesian inference: An example from the Tohoku volcanic arc, Japan. J Geophys Res Solid Earth (1978–2012). 2004;109:B10.Google Scholar
 Martinez W, Martinez A. Computational statistics handbook with MATLAB. Boca Raton: CRC press; 2001.View ArticleGoogle Scholar
 Marzocchi W, Bebbington MS. Probabilistic eruption forecasting at short and long time scales. Bull Volcanol. 2012;74(8):1,777–805.View ArticleGoogle Scholar
 Melían G, Hernández P, Padrón E, Pérez NM, Barrancos J, Padilla G, et al. Spatial and temporal variations of diffuse CO2 degassing at El Hierro volcanic system: Relation to the 2011–2012 submarine eruption. J Geophys Res Solid Earth. 2014;119(9):6,976–91.View ArticleGoogle Scholar
 Moufti MR, Matsah MI, Soliman MA, Moghazi A. Arabian plume dynamics beneath AlMadinah AlMunawwarah region and its related geohazards, final report. 2010. Technical Report ARP2679, King Abdulaziz City for Science and Technology.Google Scholar
 Moufti MR, Moghazi A, Ali K. 40Ar/39Ar geochronology of the NeogeneQuaternary Harrat AlMadinah intercontinental volcanic field, Saudi Arabia: implications for duration and migration of volcanic activity. J Asian Earth Sci. 2013;62:253–68.View ArticleGoogle Scholar
 Murcia H, Németh K, Moufti MR, Lindsay J, ElMasry N, Cronin S, et al. Late Holocene lava flow morphotypes of northern Harrat Rahat, Kingdom of Saudi Arabia: Implications for the description of continental lava fields. J Asian Earth Sci. 2014;84:131–45.View ArticleGoogle Scholar
 Needham A, Lindsay J, Smith I, Augustinus P, Shane P. Sequential eruption of alkaline and subalkaline magmas from a small monogenetic volcano in the Auckland Volcanic Field, New Zealand. J Volcanol Geotherm Res. 2011;201(1):126–42.View ArticleGoogle Scholar
 Németh K, Cronin S. Phreatomagmatic volcanic hazards where riftsystems meet the sea, a study from Ambae Island, Vanuatu. J Volcanol Geotherm Res. 2009;180(2):246–58.View ArticleGoogle Scholar
 Pallister J, McCausland W, Jónsson S, Lu Z, Zahran H, El Hadidy S, et al. Broad accommodation of riftrelated extension recorded by dyke intrusion in Saudi Arabia. Nat Geosci. 2010;3(10):705–12.View ArticleGoogle Scholar
 Perry F, Cogbill A, Kelley R. Uncovering buried volcanoes at Yucca Mountain. Eos Trans Am Geophys Union. 2005;86(47):485–8.View ArticleGoogle Scholar
 Pinel V, Jaupart C. Magma storage and horizontal dyke injection beneath a volcanic edifice. Earth Planet Sci Lett. 2004;221(1):245–62.View ArticleGoogle Scholar
 Preparata F, Hong SJ. Convex hulls of finite sets of points in two and three dimensions. Commun ACM. 1977;20(2):87–93.View ArticleGoogle Scholar
 Ripley B. Tests of randomness for spatial point patterns. J R Stat Soc Ser B. 1979;41:368–74.Google Scholar
 Rolandone F, Lucazeau F, Leroy S, Mareschal JC, Jorand R, Goutorbe B, et al. New heat flow measurements in Oman and the thermal state of the Arabian Shield and Platform. Tectonophysics. 2013;589:77–89.View ArticleGoogle Scholar
 Ross PS, Delpit S, Haller M, Németh K, Corbella H. Influence of the substrate on maar–diatreme volcanoes—an example of a mixed setting from the Pali Aike volcanic field, Argentina. J Volcanol Geotherm Res. 2011;201(1):253–71.View ArticleGoogle Scholar
 Runge M, Bebbington M, Cronin S, Lindsay J, Kenedi C, Moufti MR. Vents to events: determining an eruption event record from volcanic vent structures for the Harrat Rahat, Saudi Arabia. Bull Volcanol. 2014;76(3):1–16.View ArticleGoogle Scholar
 Sandri L, Jolly G, Lindsay J, Howe T, Marzocchi W. Combining longand shortterm probabilistic volcanic hazard assessment with costbenefit analysis to support decision making in a volcanic crisis from the Auckland Volcanic Field, New Zealand. Bull Volcanol. 2012;74(3):705–23.View ArticleGoogle Scholar
 Spӧrli K, Eastwood V. Elliptical boundary of an intraplate volcanic field, Auckland, New Zealand. J Volcanol Geotherm Res. 1997;79(3):169–79.Google Scholar
 Tait S, Taisne B. The dynamics of dike propagation. In: Fagents S, Gregg T, Lopes R, editors. Modeling volcanic processes: the physics and mathematics of volcanism. 2013. p. 33–54.Google Scholar
 Takada A. The influence of regional stress and magmatic input on styles of monogenetic and polygenetic volcanism. J Geophys Res Solid Earth (1978–2012). 1994;99(B7):13,563–73.View ArticleGoogle Scholar
 Tanaka K, Shoemaker E, Ulrich G, Wolfe E. Migration of volcanism in the San Francisco volcanic field, Arizona. Geol Soc Am Bull. 1986;97(2):129–41.View ArticleGoogle Scholar
 Tomsen E, Lindsay J, Gahegan M, Wilson T, Blake D. Evacuation planning in the Auckland Volcanic Field, New Zealand: a spatio temporal approach for emergency management and transportation network decisions. J Appl Volcanol. 2014;3(1):6.View ArticleGoogle Scholar
 Toussaint G. Solving geometric problems with the rotating calipers. Proc IEEE Melecon. 1983;83:A10.Google Scholar
 Valentine G, Hirano N. Mechanisms of lowflux intraplate volcanic fields—Basin and Range (North America) and northwest Pacific Ocean. Geology. 2010;38(1):55–8.View ArticleGoogle Scholar
 Valentine G, Perry F. Decreasing magmatic footprints of individual volcanoes in a waning basaltic field. Geophys Res Lett. 2006;33:14.View ArticleGoogle Scholar
 von Hochstetter F. Geology of New Zealand. In: Delattre T, editor. Explanation of the geographical & topographical atlas of New Zealand. 1864.Google Scholar
 Von Veh MW, Németh K. An assessment of the alignments of vents based on geostatistical analysis in the Auckland Volcanic Field, New Zealand. Geomorphology. 2009;3/2009:175–86.Google Scholar
 Wadge G, Cross A. Quantitative methods for detecting aligned points: an application to the volcanic vents of the MichoacanGuanajuato volcanic field, Mexico. Geology. 1988;16(9):815–8.View ArticleGoogle Scholar
 Wood C, Shoan W. Growth patterns of monogenetic volcano fields. Eos Trans Am Geophys Union. 1981;62:1061.Google Scholar
 Zhang D, Lutz T. Structural control of igneous complexes and kimberlites: a new statistical method. Tectonophysics. 1989;159(1):137–48.View ArticleGoogle Scholar