Brief communication: A magma depletion alternative for vent distribution in volcanic fields

The location of a volcanic vent controls an eruption's hazards, intensities, and impact. Current kernel density estimation methods of future vent locations in volcanic fields assume that locations with more past-vents are more likely to produce future-vents. We examine an alternative hypothesis that an eruption depletes the magma source, causing holes or dips in the spatial density estimate for future vent locations. This is illustrated with the Auckland Volcanic Field, Aotearoa-New Zealand, where both magmatic and phreatomagmatic eruptions have occurred, according to the vent location, with the latter resulting in more explosive eruptions and hence hazard.

1 Introduction

Volcanic fields are regions of distributed volcanism where each vent usually erupts only once, and new eruptions typically occur through new vents in new locations. Approximately 75 volcanic fields have been active in the last 10 000 years with the number of vents in an individual field ranging from tens to thousands (Global Volcanism Program, 2024). Where the next eruption will occur is vital information as the hazards, intensities, and subsequent impact depends on near-surface hydrology, geology and topography, as well as magma composition and eruptive rate (Kereszturi et al., 2014). Furthermore, their typically low average eruption rates (10⁻⁴–10⁻⁶ years; Valentine and Connor, 2015), and fertile soils, mean that settlements are often located proximal to (or on) these fields. The accuracy of existing approaches to eruption location forecasting, and indeed the fundamental underlying assumptions, remain unknown and unvalidated without any results from prospective forecasts. This technical note outlines a new alternative, magma depletion, modification of existing kernel-based spatial density estimates.

Long-term probabilistic eruption forecasting for volcanic fields is essentially spatial smoothing (Connor and Hill, 1995; Marzocchi and Bebbington, 2012; Connor et al., 2015) whereby 2D probability density surfaces are built from the location of known eruption vents. This inherently assumes that locations with more past-vents are more likely to produce future-vents. Kernels are placed coincident with past eruptive centres, and their bandwidths estimated via optimisation algorithms (Connor et al., 2018), e.g., the isotropic Gaussian kernel:

$$\begin{array}{}\text{(1)}& {\mathit{\kappa }}_{\mathrm{0}}\left(x\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{h}^{\mathrm{2}}}}\exp \left[-\mathrm{0.5}{\left({\displaystyle \frac{\left|x_i-x\right|}{h}}\right)}^{\mathrm{2}}\right].\end{array}$$

Spatial density at a point is then estimated with

$$\begin{array}{}\text{(2)}& \mathit{\lambda }\left(x\right)={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{\mathit{\kappa }}_{\mathrm{0}}\left(x\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{h}^{\mathrm{2}}N}}\sum _{i=\mathrm{1}}^N\exp \left[-\mathrm{0.5}{\left({\displaystyle \frac{\left|x_i-x\right|}{h}}\right)}^{\mathrm{2}}\right].\end{array}$$

Symbology is consistent throughout this note as follows: N: Total number of vents, h: isotropic kernel bandwidth, |x_i−x|: Distance between ith vent at x_i and location x (in 2D space), λ(x): Spatial density estimate at location x, V_i: Volume erupted at ith vent, r_i, l_i: radius and height of a volume-equivalent cylinder at ith vent.

2 Magma depletion alternative for vent distribution

We examine here an alternative (but not necessarily better) hypothesis of magma depletion, i.e., that after an eruption, the magma source at depth is depleted in this area, causing holes or dips to appear in the spatial density estimate. These mechanics are aligned to current tectono-magmatic frameworks for volcanic fields (Valentine and Perry, 2006). We assume that the available source region is spatially heterogenous, with past eruption locations representing regions of higher magma fertility (e.g., McGee et al., 2015), but as eruptions occur, they deplete their immediate locality, with larger eruptions depleting the available source region more than smaller ones. This depletion then reduces the likelihood of a further eruption from that area, although it may still be more likely than some other areas.

Three (of many potential) alternate models are provided here (Fig. 1) that consider different manifestations of how this hypothesis may be pragmatically implemented. These approaches require some changes to the standard mathematical formulae, specifically the calculation of normalisation factors. We also restrict ourselves for this note to variations on isotropic Gaussian kernels (variations of Eq. 1, implemented through Eq. 3) to aid visualisation. The kernels below are used to calculate a spatial density, using a leave-one-out cross validation approach (Bebbington, 2015) that maximises the Kullback-Leibler score

$$\begin{array}{}\text{(3)}& S=\sum _k\log {\widehat{\mathit{\lambda }}}_{}\left(x_k\right)-\int \widehat{\mathit{\lambda }}\left(x\right)\mathrm{d}x,\end{array}$$

in order to estimate bandwidths and any scaling parameters, where $\widehat{\mathit{\lambda }}\left(x_k\right)={\sum }_{i\ne k}^N{\mathit{\kappa }}_i/Z_i$ (i.e. is computed using all locations except x_k), and Z_i is a normalisation factor to ensure the density integrates to 1. Note that the formulae below are independent of bandwidth estimator (see e.g., Connor et al., 2018 for alternatives to Eq. 3). For the following, κ₀ is the base kernel (Eq. 1), α is the scaling parameter [0, ∞] (so the simpler case with α=0 is nested in the model), and to collapse the erupted (3D Dense Rock Equivalent) volume (v_i) to 2D we set

$$\begin{array}{}\text{(4)}& r_i=\sqrt{{\displaystyle \frac{v_i}{\mathrm{2}\mathit{\pi }l}}},\end{array}$$

representing a disk of radius r_i, where disk height l acts as a scaling parameter [0, ∞] and is optimised as a constant for the volcanic field.

1. Disk-based depletion

   $$\begin{array}{}\text{(5)}& \mathit{\kappa }={\mathit{\kappa }}_{\mathrm{0}}\exp \left(-\mathit{\alpha }I\left(\left|x_i-x\right|<r_i\right)\right),\end{array}$$

   where depletion is applied at all distances |x_i − x| less than r_i (Eq. 4), and with no effect at distances greater than r_i, and spatial density at a point is

   $$\begin{array}{}\text{(6)}& \begin{array}{rl}{\mathit{\lambda }}_{}\left(x\right)& =\sum _{i=\mathrm{1}}^N{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{h}^{\mathrm{2}}N{Z}_i}}\exp \left[-\mathrm{0.5}{\left({\displaystyle \frac{{\left|x_i-x\right|}_{}}{h}}\right)}^{\mathrm{2}}\right]\\ & \times \exp \left(-\mathit{\alpha }I\left(\left|x_i-x\right|<r_i\right)\right),\end{array}\end{array}$$

   where $Z_i=\exp \left(-\mathrm{0.5}\frac{r_i^{\mathrm{2}}}{h^{\mathrm{2}}}\right)+$ exp (−α)$\left(\mathrm{1}-\exp \left(-\mathrm{0.5}\frac{r_i^{\mathrm{2}}}{h^{\mathrm{2}}}\right)\right)$. This partial depletion becomes complete depletion when α=∞.

2. Inverse-volume weighted

   $$\begin{array}{}\text{(7)}& \mathit{\kappa }={\mathit{\kappa }}_{\mathrm{0}}{v}_i^{-\mathit{\alpha }},\end{array}$$

   where kernel contributions are down-weighted by volume (larger volume past eruptions contribute less than smaller ones), and spatial density at a point is

   $$\begin{array}{}\text{(8)}& {\mathit{\lambda }}_{}\left(x\right)=\sum _{i=\mathrm{1}}^N{\displaystyle \frac{v_i^{-\mathit{\alpha }}}{\mathrm{2}\mathit{\pi }{h}^{\mathrm{2}}{\sum }_{k=\mathrm{1}}^N{v}_k^{-\mathit{\alpha }}}}\exp \left[-\mathrm{0.5}{\left({\displaystyle \frac{{\left|x_i-x\right|}_{}}{h}}\right)}^{\mathrm{2}}\right].\end{array}$$

3. Distance-weighted kernel

   $$\begin{array}{}\text{(9)}& \mathit{\kappa }={\mathit{\kappa }}_{\mathrm{0}}\left|x_i-x\right|,\end{array}$$

   where kernel contributions are down-weighted by volume and multiplied by a distance term (|x_i−x|), with spatial density at a point calculated by

   $$\begin{array}{}\text{(10)}& {\mathit{\lambda }}_{}\left(x\right)=\sum _{i=\mathrm{1}}^N{\displaystyle \frac{\left|x_i-x\right|}{N\sqrt{\mathrm{2}}{\mathit{\pi }}^{\mathrm{3}/\mathrm{2}}{h}^{\mathrm{3}}{v}_i^{\mathrm{3}\mathit{\alpha }}}}\exp \left[-\mathrm{0.5}{\left({\displaystyle \frac{{\left|x_i-x\right|}_{}}{h{v}_i^{\mathit{\alpha }}}}\right)}^{\mathrm{2}}\right].\end{array}$$

Figure 1 Cross sections through the centre of the isotropic kernels. (a) Isotropic Gaussian (Eq. 1), and Inverse-volume weighted (Eq. 7). (b) Disk-based depletion (Eq. 5). (c) Distance-multiplied (Eq. 9). The dashed curves represent a bandwidth roughly twice that of the solid curves, noting that for (a), the shape of the kernel is identical for Eqs. (1) and (7) so only two variations are presented.

3 Example application: Auckland Volcanic Field, Aotearoa-New Zealand

The Auckland Volcanic Field (AVF) has c. 51 eruptive centres of varying ages (∼ 193 ka to ∼ 500 years BP; Hopkins et al., 2020), and volumes (0.0001 to 0.7 km³, Bebbington, 2015) and lies beneath the City of Auckland (population ∼ 1.7 M). Previous eruptions span the magmatic-phreatomagmatic spectrum with vent locations both on- and off-shore. The most recent eruptive centre was also the largest (Rangitoto, ∼ 500 years BP, 0.7 km³), and the most North-Easterly vent. Spatial density estimates for the AVF have been published with both isotropic (Bebbington, 2013) and anisotropic kernels (Bebbington and Cronin, 2011). For illustrative simplicity, we assume a fixed region of volcanic activity, noting that the delineation of a volcanic field boundary is problematic at best (Bebbington, 2015; Runge et al., 2015). Loglikelihoods for each of the fitted models (Fig. 2) show the fit for all models is similar, with the disk-based depletion model (Eqs. 5 and 6) providing the best-fit to the Auckland Volcanic Field data (Table 1). In terms of the Akaike Information Criterion ($\mathrm{AIC}=\mathrm{2}k-\mathrm{2}\times \mathrm{Loglikelihood}$, where k is the number of fitted parameters), the additional parameters provide either equal (inverse-volume weighted model) or significant improvement over the Isotropic Gaussian kernel (Table 1: AIC). The distance multiplied kernel is a poor fit in both relative and absolute terms.

Table 1 Model loglikelihoods, Akaike Information Criterion (AIC), and fitted parameters.

4 Discussion

This brief communication introduces a new paradigm, including example formulae and code to produce alternate spatial density estimates for volcanic fields under a magma depletion hypothesis. This shows that the baseline isotropic Gaussian kernel approach can be improved simply by depleting the likelihood in the immediate vicinity of a vent, at the very least, for the AVF. Some of these models can produce estimates with questionable degrees of smoothness (e.g., the ridges in Fig. 2b), future work could introduce additional smoothing parameters – the outstanding problem for this is how to robustly estimate the appropriate level of smoothing. The formulae trialled here can likely be improved with additional data, for example, the disk-based method could incorporate depletion geometries calculated from the modelling of geochemical data (e.g., major and trace elements).

In the AVF, the widely accepted elliptical boundary (Spörli and Eastwood, 1997; Bebbington, 2015), provides the justification for finite support, i.e., a fixed subsurface region that reflects the magma source region, from which depletion can then be modelled. This is unusual for volcanic fields (Runge et al., 2014) and is why the AVF was selected as our depletion exemplar. Additionally, vent locations in the AVF appear to be randomly located (homogeneous Poisson distribution – Le Corvec et al., 2013a), rather than clustered or uniformly spaced (Runge et al., 2015), which supports the use of isotropic (symmetric) kernels. This may suggest a uniform, homogeneous source region that is then depleted over time, unlike most volcanic fields that tend to have clustered vent locations (Le Corvec et al., 2013b).

Figure 2 Fitted spatial densities to the Auckland Volcanic Field, Aotearoa-New Zealand. (a) Isotropic Gaussian (Eq. 2), (b) Disk-based depletion (Eq. 6), (c) Inverse-volume weighted (Eq. 8), (d) Distance-multiplied (Eq. 10). Vent location symbols (triangles) scale with square root of volume (except in a as volume is not considered in the model).

The main assumption in our method is that the process exhibits spatial stationarity. To test this, we refit the four models without Rangitoto (the most recent eruption in the AVF), and got similar ranges and ordering of loglikelihoods, with expected minor changes in fitted parameters (see Supplement). Beyond the statistical properties of the model, any underlying hypothesis of spatial density needs to account for stalled eruptions (where magma has left the source region but does not reach the surface). White et al. (2006) suggest that only every 5th–10th dyke reaches the surface (intrusive:extrusive ratios). Thus, for every observable vent in the depletion model, there may be many others that have depleted the source region but that we cannot see where, or by how much. Hence, we need to assume that stalled eruptions have the same spatial distribution as actual eruptions. In contrast, a clustering model such as kernel density estimation can arise from a homogenous source by assuming that stalled eruptions occur in a different spatial pattern as actual eruptions. The inclusion of eruptive volume as a covariate may improve vent distribution, however, eruptive volumes are often accompanied by large uncertainties (e.g., Kereszturi et al., 2013), especially in older volcanic fields with more burial of the deposits and erosion. However, as long as the relative volumes between vents do not change significantly, then these uncertainties will be absorbed within model fitting parameters, noting that uncertainties in erupted volume estimates will propagate through the model, and consequently may affect uncertainties in the fitted parameters.

These proposed ideas herein have used perturbed kernel-based models (i.e., centred on vent locations), which is obviously a major restriction. The scope for models outside of this class is undelimited, but our results suggest that further reducing the likelihood clustering around previous vent locations might prove fruitful. Future work should also consider a broader family of kernels, including anisotropic, and spatio-temporal options, and the grouping of multi-vent eruptions into single events before spatial density models are deployed (Runge et al., 2014). The effect of such groupings will be most pronounced on those estimation methods that rely heavily on minimum inter-vent distances (e.g., the distance-multiplied model), where likelihoods will be over-penalised by using vents, rather than events. These alternatives may provide insights for those volcanic fields that are purely monogenetic. The appropriate use for caldera systems remains unknown and may prove an interesting avenue for future research, although account would need to be taken of the effect of the caldera structure on vent locations (e.g., Kósik et al., 2020). Volume-based temporal forecasts for polygenetic eruptive centres can be found in Bebbington (2008).

These model outputs are two-dimensional representations of a highly multi-dimensional reality. As such, the inclusion of any available covariates that inform magma source region (e.g., geochemistry – Rowe et al., 2020), system dynamics (e.g., eruptive ages or order – Bebbington and Cronin, 2011), or multi-vent events (e.g., Magill et al., 2005) may be a valuable next step.

5 Conclusions

Here, we provide an alternative depletion-based hypothesis for spatial density estimates of vent locations at volcanic fields, alongside the mathematical framework(s), derivations (see Supplement) and MATLAB codes required to implement them. Until we can directly observe, monitor and forecast magma source regions in the subsurface, we must investigate as many alternative hypotheses as are plausible, to produce transparent forecasts which together convey accurate descriptions of variability and uncertainty.