Deep Magma Transport Control on the Size and Evolution of Explosive Volcanic Eruptions

Explosive eruptions are the surface manifestation of dynamics that involve transfer of magma from the underground regions of magma accumulation. Evidence of the involvement of compositionally different magmas from different reservoirs is continuously increasing to countless cases. Yet, models of eruption dynamics consider only the uppermost portion of the plumbing system, neglecting connections to deeper regions of magma storage. Here we show that the extent and efficiency of the interconnections between different magma storage regions largely control the size of the eruptions, their evolution, the causes of their termination, and ultimately their impact on the surrounding environment. Our numerical simulations first reproduce the magnitude-intensity relationship observed for explosive eruptions on Earth and explain the observed variable evolutions of eruption mass flow rates. Because deep magmatic interconnections are largely inaccessible to present-day imaging capabilities, our results motivate the need to better image and characterize extant magma bodies.

1 Introduction

Volcanic eruptions are the surface manifestation of complex dynamics in underground systems involving one or more magma chambers, their interconnecting dykes, the surrounding rocks and aquifers, and the conduit or dyke connecting the shallow reservoir to the surface. With the exception of the volcanic crater region, the flow conditions are subsonic everywhere in the underground system, implying that any process or disturbance anywhere in the fluid system propagates both downflow and upflow, influencing the dynamics in the entire system. As a consequence, the existence and activation of multiple magmatic reservoirs during volcanic eruptions cannot be neglected without losing substantial knowledge and understanding. The occurrence of an explosive eruption as a consequence of new injections of magma into a shallow chamber has been proposed at several volcanoes worldwide (e.g., Askja, Sparks et al. (1977); Piñatubo, Pallister et al. (1992); Krakatau, Mandeville et al. (1996); Mount Ranier, Venezky and Rutherford (1997); SoufriereHills, Murphy et al. (2000); El Chichon, Tepley et al. (2000); Nevado del Ruiz, Stix et al. (2003); St. Helens, Waite and Moran (2009); Tungurahua, Samaniego et al. (2011); Eyjafjallajkull, Keiding and Sigmarsson (2012); Etna, Bonaccorso et al. (2014); Sakurajima, Araya et al. (2019)). At Mount St. Helens, a P-wave velocity model based on local earthquake data images a low-velocity zone located beneath the edifice, between about 2 and 3.5 km below sea level, that is interpreted as a shallow magma storage zone; low velocities between about 5.5 and 8 km are, instead, modeled as the location of a large magma storage volume (Waite and Moran, 2009). At Nevado del Ruiz, Londoño and Kumagai (2018) infer from 4D seismic tomography a magma chamber at 2–3 km depth, refilled by a larger magma source at 7–8 km depth. At Piñatubo, Pallister et al. (1992) report evidence for the mixing of basaltic and dacitic magmas shortly before the Plinian eruption on June 15, 1991, suggesting that the injection of basaltic magma into the base of a magma chamber might have triggered the eruption. Similar pre-sin-eruptive injections by deep magma into a shallow reservoir, and consequent magma mixing, are described for many calderic and post-calderic eruptions at Campi Flegrei (Di Renzo et al., 2011). The 39 ka caldera-forming Campanian Ignimbrite eruption, Campi Flegrei, is described as originating from a 6 km deep large magma chamber connected to shallow, 2–3 km deep, partially degassed “magma apophysis” from which the eruption initially took place (Marianelli et al., 2006). The November 23, 2013 violent explosive lava fountain eruption at Mount Etna volcano confirmed the existence of a relatively small source region located at sea level (about 3 km depth), where a volume of magma, allowed to accumulate before erupting, was invaded by more primitive, fast rising magma resulting in pressurization and subsequent lava fountaining (Bonaccorso et al., 2014).

In spite of widespread evidence suggesting the activation of multiple reservoirs during volcanic eruptions, physical models often attempt to explain the eruption dynamics by referring to just one magma chamber hosting a given magma type, progressively emptied during the eruption (e.g., Bower and Woods (1998); Folch and Martí (1998); Martí et al. (2000); Macedonio et al. (2005); Folch and Martí (2009); Colucci et al. (2014)). In all such cases, pressure decrease upon magma extraction is invariably computed, terminating the eruption when the magmatic pressure drops below the local lithostatic pressure by an amount dictated by rock tensile strength, at which conditions rock failure is assumed to occur. A volcanic eruption in such a set up is relatively simple: the eruption intensity (or mass flow-rate) invariably decreases with time; the larger the magma chamber, the longer the eruption, and the larger the mass of discharged magma. With such an approach, Bower and Woods (1998), for example, find a two orders of magnitude increase in the erupted magma volume when the chamber volume is also increased by two orders of magnitude. More recently, Melnik and Costa (2014) explored the long-term (years) dynamics of a system made of two magma chambers connected each other by a dyke and connected to the Earth surface through a conduit: the study shows a cyclic behavior dependent on influx rate, volume, shape and degree of connectivity of the two reservoirs. Reverso et al. (2014) develop an analytic model for the long-term dynamics of two connected magma chambers surrounded by an elastic medium and fed by a constant basal magma inflow, modeling the magma during the ascent as an incompressible fluid.

In this work we investigate, through a global sensitivity analysis, the short-term (from hours to weeks) eruptive dynamics for a plumbing system characterized by a deep large reservoir hosting andesitic magma, connected to a shallow dacitic magma chamber. The available data on real eruption intensity and magnitude mostly involve these compositions, or others very similar belonging to the same andesitic suite (see below the discussion related to Figure 1 at the end of section 3-1). The condition whereby a shallow, more chemically (dacitic in this case) evolved magma is invaded by more primitive (andesitic) magma reflects the typical occurrence at most volcanoes worldwide (references above). The ascent of magma between the two reservoirs, and from the shallow chamber to the surface, is modeled by a 1D, steady, isothermal, compressible, multiphase, multicomponent model. We show that accounting for the concomitant action of different magmatic reservoirs leads to a new world of volcanic eruption dynamics, and new understanding of several aspects of the observed eruption dynamics and evolutions. We first analyze the model through a global uncertainty and sensitivity analysis, leading to the identification of the major controls by the model input quantities. The controlling roles of such quantities are then explored through dedicated parametric studies. The efficiency of inter-chamber magma transfer emerges as a new, major quantity largely controlling the eruption dynamics and evolution, and the conditions for eruption termination.

FIGURE 1

FIGURE 1. Observed ($•$) vs. simulated ($\circ$) magnitude-intensity distribution. Observed data from the LaMEVE database (Crosweller et al. 2012; Brown et al. 2014; www2.bgs.ac.uk/vogripa/, accessed in February 2019). Magnitude and intensity as in Pyle (2015), also defined in the Methods. The space of variables has been sampled through a Latin Hypercube Sampling (LHS) technique. Among the simulated points, blue/red refers to cases for which eruption closure is due to the achievement of conditions of rock failure at shallow/deep chamber level.

2 Numerical Model

2.1 Modeling Approach

Figure 2 (left) shows the simulation setup. A deep, large reservoir hosting andesitic magma is connected through a vertical planar dyke to a shallow, smaller chamber hosting more chemically evolved dacite, connected to the surface through a cylindrical conduit. The input parameters (volume and depth of the chambers, conduit diameter, dyke thickness, volatile abundance in the two magmas) are varied in the ranges reported in Table 1, reflecting more frequently reported values for explosive eruptions in the literature (e.g., the references at the beginning of the Introduction section). Table 2 reports the composition of the employed dacitic and andesitic melts, which affects the major properties density and viscosity, and the multi-component volatile saturation surface. Water and carbon dioxide are considered as volatile species partitioning between the melt and gas phase. Magma properties and gas-melt volatile partition are locally defined on the basis of the $P-T$ composition conditions, and evolve together with system evolution. Our computations refer to only the quasi-steady phases during which large volcanic columns develop and evolve, with or without partial or total column collapses and generation of pyroclastic flows (Neri et al., 1998). Provided that the system as a whole evolves on a time scale sufficiently longer than the transit time of magma along the simulated domain, system evolution can be approximated as a series of discrete steady flow steps (Macedonio et al., 2005; Colucci et al., 2014), with boundary conditions at each step depending on the previous dynamics and reflecting system evolution. The distribution of flow variables and properties along the deeper dyke and shallow conduit during each discrete steady flow step is computed through the model in Papale (2001), summarized in the Appendix. Magma chamber conditions are determined through a lumped approach, and evolve as a response to progressive magma withdrawal (deep chamber) and balance between andesitic magma inflow and mixed magma outflow (shallow chamber). The computations are terminated when the pressure somewhere in the system falls below lithostatic by a quantity exceeding the rock tensile strength (Bower and Woods, 1998; Gudmundsson, 2007; Kabele et al., 2017) leading to rock failure and either eruption closure or transition to highly transient dynamics. When that happens at chamber level, caldera collapse is expected (Gudmundsson, 2007; Kabele et al., 2017).

FIGURE 2

FIGURE 2. Simulation setup. (A) Simulation domain. (B–D) Example calculations at different times.

TABLE 1

TABLE 1. Input parameters used in the simulations. Magma temperature is taken to be constant along the entire simulated domain, and equal to 950°C *Composition reported in Table (2).

TABLE 2

TABLE 2. Chemical composition of the employed andesitic and dacitic melts (wt%).

A point-like lumped system approximation for the two magma chambers keeps the overall complexity to manageable levels, at the same time eliminating further arbitrary variables and allowing the extraction of first order controls on the eruption dynamics by the composite nature of the plumbing system. Boundary conditions for dyke and conduit flow reflect magma withdrawal (deep chamber) and balance between magma input and output under instantaneous magma mixing (shallow chamber). The upflow conditions for dyke flow correspond to deep chamber conditions. Pressure in the shallow chamber represents both the downflow boundary condition for dyke flow, and the upflow boundary condition for conduit flow, and together with the evolution of magma composition upon mixing in the shallow chamber, it couples dyke and conduit flow dynamics, and ultimately results in consistent evolution of the entire simulated system from the deeper chamber to volcanic conduit. In details:

1. At each time step, the equations for the cylindrical upper conduit and for the planar dyke connecting the deep reservoir to the shallow chamber are solved, with boundary conditions given by pressure and composition in the two chambers.

2. At time zero, we set a lithostatic pressure boundary condition for both chambers. Initial pressurized conditions are also possible, but they are not considered here in order to keep to a minimum the number of independent variables characterizing the system. The conduit/dyke exit boundary conditions are determined by either choked flow or ambient pressure (atmopheric pressure for the conduit, and shallow chamber pressure for the dyke), as part of the solution of the simulated dynamics. Note that time zero does not correspond to the start of the eruption, which would be characterized by highly transient dynamics during which the volcanic conduit and eruption plume are established. These transient dynamics are not considered in our modeling. Accordingly, time zero corresponds to a reference time when steady flow conditions have been established along the entire simulated domain from the deep chamber to the surface.

3. For later times, composition and pressure inside each one of the two chambers are computed on the basis of mass inflow/outflow, and constitute the new boundary condition for the numerical calculations of dyke and conduit flow. For dyke flow, the third, not simulated dimension has been fixed to 100 m when converting from mass flow rate per unit length (computed through dyke ascent modeling) to mass flow rate, and from that, to chamber mass loss or gain (see below). In particular, after update of mixture density given by

${\rho }_m^{t+\Delta t}=\frac{m^{t+\Delta t}}{v^{t+\Delta t}},(1)$

where ${\rho }_m$ is mixture density, t is time, m is mass inside the chamber, v is chamber volume, the corresponding pressure is determined. At thermodynamic equilibrium assumed in the simulations, mixture density and pressure are univocally related via the real equations of state for the liquid and gas phases and the multi-component gas-melt equilibrium model employed (Papale et al., 2006). Such a relationship is non-linear, requiring a numerical solution for pressure. That implies determination of the multi-component volatile saturation surface and melt/gas densities at chamber conditions.

Under the quasi-steady flow assumption, the mass inside the chamber is approximated as

$m^{t+\Delta t}\approx m^t+\left({\dot{m}}_{in}^t-{\dot{m}}_{out}^t\right)\Delta t,(2)$

where $\dot{m}$ is mass flow rate. Chamber volume is approximated by

$v^{t+\Delta t}\approx v^t+\frac{v^t\Delta P^t}{\beta },(3)$

where β is the effective bulk modulus of the elastic walls surrounding the chambers. By repeating the same simulations with $\beta ={10}^{10}$ Pa (Bower and Woods, 1998) and $\beta =\infty$ (non-deformable chamber walls), we have found that inclusion of rock elasticity and associated chamber volume changes produce minor or negligible differences in the simulated dynamics (see the Supplementary Figure S1). Therefore, we refer to only the latter case ($\beta =\infty$, or non-deformable chamber walls) throughout this paper.

4. The time step $\Delta t$ used to advance along the sequence of steady states is determined at run time as twice the longest among the transit times in the conduit and dyke.

$\Delta t=2\ast t_{trans}.(4)$

The transit time $t_{trans}$ into the conduit/dyke is computed by space integration of the inverse of velocity.

5. The cycle stops when the local difference between the lithostatic and magmatic pressure anywhere in the computational domain exceeds the tensile strenght of rocks.

6. A-posteriori evaluation of the quasi-steady assumption is done by checking if the following criterion is satisfied:

$t_{\dot{m}}\gg t_{trans},(5)$

where $t_{\dot{m}}$ is the characteristic time of variation of the mass flow rate over $\Delta t$.

$t_{\dot{m}}=\frac{{\dot{m}}^{t+\Delta t}}{\left|{\dot{m}}^{t+\Delta t}-{\dot{m}}^t\right|}\Delta t.(6)$

In spite of its relative simplicity, such a set up captures the fundamental characteristic of many volcanic plumbing systems characterized by separated reservoirs located at different depth and hosting different magmas, being simultaneously activated during an eruption. Comparison with real eruptions mostly involves the quantities Magnitude and Intensity. Magnitude (M) is computed as in Pyle (2015) as $M=\log \left[m_T\left(kg\right)\right]-7$, where $m_T$ is the total mass erupted; in the simulations, $m_T$ is determined by integrating in time the computed mass flow rate. Intensity (I) refers to the eruption mass flow rate, and is defined as in Pyle (2015) as $I=\log \left[\dot{m}\left(kg/s\right)\right]+3$, where $\dot{m}$ is the average mass flow rate.

Point 5 above requires further discussion. As explained above, we stop a simulation run when the pressure somewhere in the simulated domain, extending from the deep chamber to the surface, falls below the critical threshold dictated by rock tensile strength. Beyond such conditions local rock failure is expected, terminating the sub-steady phase of the eruption, or the eruption itself. If that happens at magma chamber level, caldera collapse may occur. In the real world the processes are expected to be more complex. As an example, initiation of caldera roof collapse may cause the magmatic pressure to increase back to values sufficient to avoid further collapse, and a sub-steady phase may be restored until further magmatic pressure changes lead to a new instability. Similarly, if large pressure decrease occurs locally in the volcanic conduit (typically close to magma fragmentation where the pressure gradient is the highest, Papale and Dobran (1994)), that may lead to local conduit wall erosion and conduit shape changes (Macedonio et al., 1994; Aravena et al., 2017). Additional complexities may come from the mechanics of rock failure (Kabele et al., 2017), especially in layered environments (Gudmundsson, 2007), whereby the generation of ring faults and the occurrence of caldera collapse depend on a number of other factors including chamber shape. The development of a yield strength in a crystal-rich magma reservoir affects the size of the magma body and hence the evolution of pressure and eruption rate, leading to increasing eruptions rates through time (Karlstrom and Manga, 2012). Because we only refer to steady flow dynamics and do not include in the model fracturing and faulting, that are beyond the scopes of this work, we stop our simulations when the computed pressure changes do not guarantee further validity of the steady flow assumption. Such an approach is similar to the one employed in previous works: Bower and Woods (1998) assume a critical under-pressure threshold in the magma chamber of 35 MPa as a condition for termination of steady conduit flow, and with the same approach Macedonio et al. (2005) use the very large value of 80 MPa. In his comprehensive analysis of ring-fault formation in volcanic areas, Gudmundsson (2007) employs a critical under-pressure of only 5 MPa as a potential trigger for ring-fault formation above a magma chamber. More recently, Kabele et al. (2017) develop a finite element model for the mechanical interactions between a magma chamber and its surrounding rocks, employing a rock tensile strength of 8 MPa. Although the real dynamics of caldera formation may be complex, and they are beyond the scopes of this work, in the following we adopt a fixed threshold of 20 MPa, and refer to the cases where the simulation is terminated because of pressure decrease in one of the two magma chambers, as “shallow chamber collapse” or “deep chamber collapse” for the shallow and deep magma chamber, respectively.

The dyke thickness parameter in our modeling requires particular care. In fact, magmatic reservoirs located at different depths may communicate to each other through a far more complex system of interconnections, possibly involving multiple dykes and intermediate storage zones (Edmonds et al., 2019). We will show in the following that dyke thickness is a major factor controlling the eruption dynamics and evolution, its role being that of balancing between magma withdrawal from the deep chamber, and magma inflow and outflow at shallow chamber level. Accordingly, dyke thickness in our simulations should be taken as a proxy for the overall efficiency of magma transfer across the variably complex system connecting the magmatic reservoirs.

2.2 Uncertainty and Global Sensitivity Analysis

In order to explore the dynamics of the composite system represented in Figure 2, the model described above is applied to the range of input parameters reported in Table 1. An uncertainty and global sensitivity analysis is performed in order to extract the input parameters controlling the eruption dynamics.

Uncertainty and sensitivity analysis investigate how a subset of the model assumptions, in the present case constituted by the input parameters in Table 1, affects the model output. The uncertainty analysis consists in studying how the variability of the input propagates into the distribution of the output. This can be done by repeatedly running the model with different input values sampled over their plausible ranges, and studying the distribution of model output. The sensitivity analysis studies how the variability observed in model outputs has to be ascribed to the different input parameters, thus providing information on the relative importance of the input quantities (Saltelli et al., 2019). To perform a valid sensitivity analysis when a model features nonlinearities and interactions, as it is largely the case with the present model, a global approach is required. In the global sensitivity analysis, the variation of the output is obtained by varying all the input factors simultaneously, so that the sensitivity to each input factor considers the direct influence of that factor as well as the joint influence due to interactions. It has been noted (Saltelli et al., 2019) that in the Earth Science literature many sensitivity analyses are still based on changing the input factors one at a time, leaving the space of the input variables largely unexplored.

In the present analysis we refer to three model outputs, represented by the eruption magnitude (which measures how much magma is discharged), the eruption intensity (which measures how fast the magma is discharged), and the location (shallow vs. deep) where rock failure terminating the simulation occurs. The input parameters are reported in Table 1.

We first perform an uncertainty analysis through a latin hypercube sampling of the space of the input quantities using the Dakota software (Adams et al., 2014). Then, two types of global sensitivity analysis are performed: a regional sensitivity analysis and a novel density-based sensitivity analysis called PAWN (Pianosi and Wagener, 2018).

Regional Sensitivity Analysis (RSA, Pianosi et al. (2016)) aims at identifying regions in the input space that correspond to particular values of the output. By comparing the empirical cumulative distribution functions of an input factor associated with output values above or below a given threshold (for continuous output values such as magnitude and intensity), or associated with different boolean-type outputs (such as “shallow” or “deep” rock failure), we can determine regions of the input that have different influence on the considered output. The divergence between the two distributions is measured by the Kolmogorov-Smirnov statistics.

$KS=\underset{x_i}{\max }\left|F_{x_i|y_a}\left(x_i|y\in Y_a\right)-F_{x_i|y_b}\left(x_i|y\in Y_b\right)\right|,(7)$

where $F_{x_i|y_a}$ and $F_{x_i|y_b}$ are the cumulative distribution functions of the input parameter $x_i$ associated with output values, $y_a$, $y_b$ (above/below a given threshold, or behavioral/non-behavioural for boolean-type values). We estimate KS with a two-sample Kolmogorov-Smirnov test using the Matlab function kstest2. The larger the value of $KS$ (that is, the more divergent the two cumulative distributions), the more the specific input has an influence on the considered output (however, while a small or vanishing value of $KS$ is a necessary condition to approach insensitivity, it may not be sufficient because interactions between input factors cannot be detected).

The PAWN analysis (Pianosi and Wagener, 2015; Pianosi and Wagener, 2018) consists of a distribution-based method that uses the cumulative distribution function of the output and is particularly suitable when the output distribution is skewed or multi-modal, as in our case for the magnitude and intensity outputs. The PAWN sensitivity index for the ith input factor is defined as

$S_i=\underset{x_i}{stat}\underset{y}{\max }\left|F_y\left(y\right)-F_{y|x_i}\left(y|x_i\right)\right|,(8)$

where $F_y$ and $F_{y|x_i}$ are the unconditional and conditional cumulative distribution functions of the output y, and $stat$ is a statistic (e.g. maximum, median or mean). Notice that the inner maximum in Eq. 8 corresponds to the Kolmogorov-Smirnov (KS) statistic. However, it is worth noting that the RSA considers the KS statistics of the inputs (Eq. 7), while Eq. 8 considers the KS variations in the cumulative distribution of the output. The RSA studies how the input distributions vary when conditioning the output, while the PAWN analyses how the output distributions vary when conditioning the input. In this work we estimate the PAWN sensitivity indexes relative to the model outputs represented by eruption magnitude and intensity, following the strategy recommended in Pianosi and Wagener (2018), using the Matlab SAFE toolbox (https://www.safetoolbox.info/). According to Pianosi and Wagener (2018), to assess robustness and convergence of the present analysis, we repeat the calculations for different values of the sample size (N) and the tuning parameter n (number of conditioning intervals).

3 Simulations

The right panels in Figure 2 show one example of calculations (only a few times are reported). At each time step the flow variable distributions display well-known behaviors [e.g., Papale (2001)] reflecting the establishment of a friction-dominated bubbly flow region characterized by large gradients of the flow quantities culminating in magma fragmentation, followed by an inertia-dominated gas-particle region from fragmentation to conduit exit. Lumped pressure in the shallow magma chamber corresponds at any time to both dyke exit and conduit base pressure, coupling the dynamics in the entire domain. Mixing between compositionally different magmas in the shallow chamber causes a progressive shift in the composition, thus in the properties of the magma flowing along the conduit, and results in a jump of gas volume fraction, of (not shown) magma density, and by continuity, of ascent velocity. Progressive withdrawal of magma from the deep chamber, controlled by the evolution of quantities over the entire domain, causes progressive decrease of pressure at depth, in turn influencing the overall dynamics.

3.1 Uncertainty Analysis

Figure 1 shows the computed variability in the two major output quantities magnitude and intensity, plotted together with a large number of observed data points from the literature [black dots, from the LaMEVE database (Crosweller et al. (2012)]. The open circles are the results from the present numerical simulations, each circle originating from a sequence of magma ascent simulations such as the one exemplified in Figure 2. Red and blue circles correspond to cases for which eruption termination occurs as a consequence of excessive pressure drop in the deep and shallow magma chamber, respectively. Significant depressurization exceeding rock tensile strength at volcanic conduit level only emerges for narrow conduits <50 m, in agreement with the conclusions by Aravena et al. (2018) who suggest dominantly unsteady dynamics at such unstable conditions. Accordingly, we have limited our analysis to conduit diameters >50 m (see Table 1). Conversely, over the explored conditions in Table 1, pressure decrease along the dyke connecting the two magma chambers is never found to approach rock tensile strength. The correspondence between the trend given by the data and that from the results of the simulations is notable. There are 719 simulated cases in the figure, which are obtained by automatically sampling (through the Dakota software, as explained in the Methods section) the input quantities over the reasonable ranges in Table 1. The modeling approach illustrated above, that results in a satisfactory representation of the real magnitude-intensity space, suggests that our model is able to capture the first order aspects of real explosive eruptions, at least in terms of total erupted mass, mass discharge rates, and eruption duration. The data describe a low-intensity region $\left(I<10.5\right)$ toward the low-M side of the diagram, which is not reproduced by our model. For such low-intensity eruptions, a reported dominance of unsteady dynamics [e.g., Miller and Chouet (1994); Pioli et al. (2008); Pistolesi et al. (2016); Gudnason et al. (2018)] or other processes not accounted for by the modeling (e.g., permeable gas loss or interaction with external aquifers) may explain their absence in the simulated set. Eruptions not associated to the andesitic suite represent less than 15% of the data in Figure 1. Within the andesitic suite, the eruptions discharging andesitic, dacitic, or andesitic-dacitic magma constitute about 50% of the data. Eruptions with dominance of more (rhyolitic) or less (basaltic) chemically evolved compositions, or belonging to the alkaline series (trachy-andesites to trachytes and phonolites) may contribute to explain a larger spread of the data with respect to the simulations. Although with ample superposition, and in line with the general observations, deep chamber collapse is mostly associated to large magnitude eruptions, and characterizes most of the globally impacting eruptions with magnitude seven or larger.

3.2 Global Sensitivity Analysis

Figures 3–5 show the results of the regional sensitivity analysis (RSA) in relation to the three considered major model outputs outlined above, namely, location of the expected rock failure terminating the simulation (Figure 3), eruption magnitude (Figure 4), and eruption intensity (Figure 5). In the diagrams of these figures, the larger the distance, quantified by $KS$ (Eq. 7), between the two cumulative distributions represented by the blue and red lines, the more the output is sensitive to the input parameter under consideration. The quantity h reported on each panel is the logical test result, returned by the kstest2 Matlab function, on the null hypothesis that the two distribution functions are the same at the $5\text{\%}$ significance level.

FIGURE 3

FIGURE 3. Regional analysis results for occurrence of eruption closure by shallow (blue lines) or deep (red lines) chamber collapse. ks (Eq. 7) quantifies the difference between the two distribution functions; the larger the value of ks, the more influential the input parameter; h is the hypothesis test result (h = 1 indicates that the test rejects the null hypothesis that the two distribution functions are the same at the $5\text{\%}$ significance level).

FIGURE 4

FIGURE 4. Regional analysis results for eruption magnitude less than (blue lines) or larger than (red lines) the value of 6.5. See Figure 3 for the interpretation.

FIGURE 5

FIGURE 5. Regional analysis results for eruption intensity less than (blue lines) or larger than (red lines) the value of 11.5. See Figure 3 for the interpretation.

For each one of the major output quantities in Figures 3–5, dyke thickness turns out to be the input quantity associated with the largest $KS$ values (0.54–0.63), with larger dykes preferentially associated with the occurrence of rock failure at deep chamber level (Figure 3), larger eruption magnitude (Figure 4), and larger eruption intensity (Figure 5). In other words, dyke thickness emerges as the individual input quantity with the strongest control on the computed eruption dynamics and evolution. Over the large range of conditions evaluated, rock failure (Figure 3) occurs at shallow chamber level for dyke thicknesss less than about 10 m, and invariably at deep chamber level if the dyke thickness is larger than about 25 m. With reference to eruption magnitude and intensity (Figures 4, 5), the only other strong relationship $\left(KS>0.4\right)$ which emerges from the analysis is the well-known association of larger conduit diameters with larger eruption intensity, described since the early conduit flow models [e.g., Wilson et al. (1980)] and confirmed here as holding also when conduit flow is associated with multiple magma chamber withdrawal and inter-chamber magma transfer dynamics. Conversely, a number of input quantities other than dyke thickness emerge as exerting a control on the occurrence of shallow vs. deep rock failure (Figure 3), and among them, the most relevant appears to be the volume of the deep magma chamber $\left(KS=0.34\right)$.

Figure 6 shows the PAWN sensitivity indexes relative to the output quantities magnitude and intensity, calculated as the median $KS$ across conditioning intervals ($\text{stat}=\text{median}$ in Eq. 8). Bootstrapping is used to estimate the 95% confidence interval (vertical line) and mean value of the median (bar height) of each PAWN index. The red line shows the KS of the dummy parameter providing an indication of the reliability of the computed indexes: if an input quantity is associated with a sensitivity index significantly larger than the dummy sensitivity, one can confidently conclude that such an input quantity is indeed influential, otherwise its non-zero sensitivity may be due to approximation errors. As Figure 6 shows, dyke thickness emerges again as being highly influential on both eruption magnitude and intensity, and conduit diameter as highly influential on eruption intensity, confirming the RSA analysis at Figures 3–5.

FIGURE 6

FIGURE 6. PAWN sensitivity indices as the median KS across conditioning intervals (stat = median in Eq. 8), for sample size N = 719 and number of conditioning intervals n = 10. Vertical line represents the 95% confidence interval estimated by bootstrapping and the bar height indicates the mean value of each PAWN index. The red line is the KS of the dummy parameter (values corresponding to bar heights lower than the red line are unreliable). The analysis is obtained with a sample size N = 719, and setting a number of conditioning intervals n = 10 sufficient to reach the convergence of the PAWN sensitivity indices (see also Supplementary Figures S2–5).

Overall, the picture emerging from the RSA and PAWN analyses is that of a major control exerted by the dyke thickness parameter on the overall simulated eruption dynamics and evolution, revealed here in terms of eruption intensity and magnitude, and on simulation closure as due to excessive depressurization within the shallow vs. deep magmatic system. Large dyke thicknesses (or high efficiency of magma transfer between deep reservoirs) tend to be associated with higher eruption magnitudes and intensities (although conduit diameter plays a major control on the latter), and to eruption closure by rock failure at deep chamber level (with other quantities playing a role on the latter, and among them, the volume of the deep magma chamber). In the following we analyze the major controls outlined here by making use of parametric studies primarily on dyke thickness and magma chamber volumes and illustrate the mechanistic aspects of such controls, leading to the identification of first relevance features of explosive eruption dynamics that were previously overlooked or unknown.

3.3 Control by Deep Interconnections

Figure 7 shows the computed evolution of some relevant eruption quantities, in a parametric study where the dyke thickness is varied from one simulation to another, while the simulation set up including all other input quantities is let unchanged. With small dyke thickness magma transfer across reservoirs is poorly efficient, thus the eruption discharge rate is not compensated by deep magma arrival (Figure 7D, cold colors). As a consequence, shallow chamber pressure decays rapidly (Figure 7A), while deep chamber pressure changes only slightly (Figure 7C) reflecting minor magma withdrawal from depth. At such conditions eruption closure quickly follows from shallow chamber pressure drop after progressive decline of the eruption discharge rate, resulting in a low magnitude eruption (Figure 7B).

FIGURE 7

FIGURE 7. Computed pressure in the shallow and deep magma chambers, and mass flow rate along the conduit and dyke. Calculations corresponding to dyke thicknesses from 4 (cold colors) to 30 (hot colors) m, with regular dyke thickness increment by 2 m. Other simulation conditions: deep chamber volume = $200{\text{km}}^3$; shallow chamber volume = $3{\text{km}}^3$; depth of the shallow chamber = 3 km; conduit diameter = 90 m; total volatile contents: H₂O in dacite = 4.5 wt%, CO₂ in dacite = 1 wt%, H₂O in andesite = 5 wt%, CO₂ in andesite = 2 wt%. The quantity M in panel 3b is the computed eruption magnitude. The corresponding evolutions of magma density and gas volume fraction are reported in the Supplementary Figure S7; the corresponding erupted magma composition is reported in Figure 9.

With larger dyke thickness (10–12 m in Figure 7), the increased efficiency of deep mass transfer sustains the eruption discharge rate (Figures 3B,D, light blue/cyan curves), leading to longer sub-steady eruption phase (up to about 8 h in the figure) with nearly constant mass flow rate. Deep chamber pressure decrease becomes significant (Figure 7C) due to increased magma withdrawal, eventually reaching the rock collapse threshold (starting from the cyan curve, or dyke thickness of 12 m).

With further increase in dyke thickness, the mass flow rate along the dyke can overcome the one along the shallow conduit (hot-coloured lines in Figures 7B,D), causing shallow chamber pressurization (Figure 7A) and further deep chamber depressurization (Figure 7C). In turn, the former causes the eruption discharge rate to increase (Figure 7B). However, shallow pressure increase and deep pressure decrease rapidly lead to decreased mass flow rate along the dyke (Figure 7D), as the concomitant action of increasing downflow and decreasing upflow pressure driving dyke flow. Therefore, shallow chamber pressurization is self-buffered, as it induces changes in conduit and dyke flow ending up with less pressure increase. Accordingly, after initial shallow pressurization and intensity increase, the eruption evolves toward shallow pressure decrease (Figure 7A) and about constant or slightly decreasing mass flow rate (Figure 7B).

The eruption evolutions represented in Figure 7 correspond to a fixed deep chamber volume of $200{\text{km}}^3$. A similar figure for a deep chamber volume of $20\text{,}000{\text{km}}^3$ is provided in the Supplementary Figure S6.

3.4 Shallow Versus Deep Controlled Eruption Regimes

The dynamics illustrated above show that besides explaining the observed magnitude-intensity relationship (Figure 1), the interplay between magmatic reservoirs provides the physical framework explaining observed variabilities in the evolution of eruption intensity, with mass flow-rates that can decrease, oscillate, remain constant, or increase depending on such an interplay. Non-linear relationships dominate the dynamics: when the dyke thickness is small, its increase produces longer eruption duration; conversely, if the dyke thickness is sufficiently large, its further increase results in shorter eruption duration. The longest durations correspond to near-balance of conduit and dyke mass flow rates (Figures 7B,D). The eruption magnitude (reported in panel 7b) is also non-linearly related to dyke thickness: for small dyke thickness, the increase in dyke thickness results in larger eruption magnitude, whereas for dyke thicknesses sufficiently large to cause eruption termination by deep chamber collapse, the eruption magnitude becomes poorly dependent on further dyke thickness changes.

The above depicts two regimes for explosive volcanic eruptions (Figure 8): a shallow-controlled regime determined by low inter-chamber magma transfer efficiency (= small dyke thickness) and evolving to shallow chamber collapse, characterized by rapid shallow chamber pressure decrease, minor changes in the deep reservoir, and positive dependence of eruption magnitude and duration on deep magma transfer efficiency; a deep-controlled regime determined by large inter-chamber magma transfer efficiency and evolving to deep chamber collapse, characterized by initial shallow pressure and eruption mass flow rate increase, significant pressure decrease in the deep reservoir, negative dependence of eruption duration and poorly negative or no dependence of eruption magnitude on magma transfer efficiency.

FIGURE 8

FIGURE 8. Role of deep interconnections. Calculated eruption magnitude (top panel) and duration (bottom panel) as a function of dyke thickness, for different deep chamber volumes and other conditions as in Figure 7. The different symbols refer to different shallow chamber volumes: circles = $1\text{.}5{\text{km}}^3$; squares = $3{\text{km}}^3$; diamonds = $15{\text{km}}^3$; stars = $30{\text{km}}^3$. Symbol colors as in Figure 1: blue and red indicate eruption closure by shallow or deep chamber collapse, respectively.

The critical dyke thickness (or inter-chamber magma transfer efficiency) at regime transition is non-linearly related to other system conditions. For small deep chamber volume of $200{\text{km}}^3$ the critical thickness decreases with increasing shallow chamber volume, implying that the combination of small deep and large shallow chamber volume is most likely to result in deep chamber collapse. Furthermore, the extent of the deep controlled regime progressively shrinks when increasing the deep chamber volume, with the largest volume of $20\text{,}000{\text{km}}^3$ requiring dyke thicknesses, and comparable magma transfer efficiency, of as much as $\text{>}15\text{m}$. These elements may represent additional factors decreasing the relative abundance on Earth of very large calderas, compared to small-medium ones, besides the likely lower abundance of extremely large magma chambers that may originate them. As anticipated, in the deep-controlled regime and for sufficiently large deep chamber volume the eruption magnitude becomes independent from both shallow chamber volume and dyke thickness. For the largest deep chamber volume in Figure 8, corresponding to $20\text{,}000{\text{km}}^3$, the computed eruption magnitudes extend to approach the value of eight corresponding to volcanic super-eruptions, with corresponding eruption durations in the range 250–450 h (or roughly 10–20 days), and longest durations found to characterize eruptions close to regime transition where inflow and outflow with respect to the shallow chamber nearly compensate each other. For the smallest shallow chamber volume in Figure 8, corresponding to $1\text{.}5{\text{km}}^3$, the critical dyke thickness at regime transition remains constantly close to 15 m over a three orders of magnitude range of the deep chamber volume. Conversely, for the largest shallow chamber volume $30{\text{km}}^3$ the transition between shallow-controlled and deep-controlled regimes requires a dyke thickness of only about 7 m for a deep chamber volume of $200{\text{km}}^3$, increasing to about 15 m, largely independent from shallow chamber volume, for a deep chamber volume of $20\text{,}000{\text{km}}^3$.

4 Discussion

4.1 Relevance of Deep Magma Transfer Efficiency

This study illustrates how the rich variability of observed magnitude, intensity, duration, and evolution of explosive volcanic eruptions on Earth can follow from the interplay between different magmatic reservoirs, as they are generally constrained by the bulk volcanological and petrological reconstructions from the products of past eruptions. In our study we vary several parameters over a large portion of their reconstructed natural variability. The performed sensitivity analysis allows consideration of the variablity of the eruption dynamics and evolution over a wide range of possible system conditions. We have limited our analysis to the andesitic magmatic suite, which represents the vast majority of subaerial explosive volcanism on Earth, and to a representative situation where more chemically evolved magma, that we fixed to dacite, occupies the shallow portion of a volcanic plumbing system which hosts andesitic magma at larger depth. Use of slightly different compositions such as rhyolite or basaltic andesite is not expected to change appreciably the results as long as the more viscous, chemically evolved, volatile-depleted magma is placed at shallow level, as the totality of petrologic modeling and reconstructions suggest. Furthermore, large changes in many other eruption parameters, reported in Table 2 and resulting in more than 700 simulations, ensure that a comparably large proportion of the space of eruptive conditions is covered. Each individual simulation illustrates an eruption evolution path such as those represented by each individual curve in Figure 7, and each one produces a point in Figure 1 (or in Figure 8). In analyzing such a vast suite of numerical results, we concentrate on a few parameters that emerge from the sensitivity analysis as playing a more effective control on first order characteristics of the resulting eruptions. In particular, we concentrate on eruption intensity and magnitude, as well as on the occurrence of deep vs. shallow pressure drop, beyond a sufficiently large threshold assumed to result in local rock collapse and either eruption closure or termination of the sub-steady eruption phase. A new quantity emerges as playing a major control on the overall eruption dynamics and evolution, represented by the efficiency of deep interconnections between separated magmatic reservoirs, exemplified in our study by the dyke thickness parameter. By varying the interconnection efficiency, we have highlighted the existence of two eruption regimes mostly controlled by the shallow or by the deep magmatic system, characterized by markedly different trends and likely to terminate as a consequence of shallow vs. deep failure of the rock system surrounding the chamber. The analyzed conditions span a range extending toward the two end-member cases constituted by isolated, relatively small, shallow magma chamber on the side of minimum dyke thickness (minimum efficiency of inter-chamber magma transfer), and large, deep chamber directly connected with the surface on the opposite side of very efficient deep magma transfer, or very large dyke thickness. It is, however, in the large region between such extremes that the richest dynamics and evolutions emerge, reflecting the concomitant action, and interaction, between the deep and shallow portions of the volcanic plumbing system.

Before considering some relevant implications of the above results, it is worth discussing the delicate balance between crediting and interpreting the numerical results on one side, and considering the limitations of the study on the other side. Complex systems are usually approached by first isolating some relevant first order aspects, then progressively increasing complexity so to construct a progressively more robust reference for understanding and interpretations. That approach is not limited to numerical investigation, rather, it characterizes virtually any scientific approach, e.g., the experimental one, where individual experiments are similarly planned to simplify real systems and isolate, for as much as possible, specific elements of interest. Because the aim is to achieve conclusions that are as much general as possible, both analogue and numerical experiments are first designed to minimize unnecessary complications. That is usually done by referring to a sort of archetype system conditions; e.g., a vertical cylinder for a volcanic conduit, a vertical planar fracture for a dyke, a sphere or ellipsoid for a magma chamber, etc. Once the first order aspects of such archetype systems are understood, one can then move to consider further complications, e.g., variations in the geometry, inclination, number etc. of the domain elements, higher order aspects of the physics involved, etc., obviously subject to constraints by technical feasibility, including finite time and resources available.

In the present case, the objective is that of systematically exploring the first order aspects of volcanic eruption dynamics involving magma initially residing at different depths in the volcanic plumbing system, an occurrence which is revealed by petrology and geochemistry as being virtually ubiquitous, but that it is still heavily under-explored in terms of controls on eruption conditions and evolution. Accordingly, our archetype simulation domain has been designed to capture the first order aspects of a composite volcanic system on one side, and to minimize unnecessary complications on the other side. Besides satisfying the above concepts, avoidance of further complications in the simulated system also reduces the computational time to within manageable amounts, allowing extensive investigation of the space of variables for meaningful investigation as well as for sensitivity analysis. In particular, we take advantage of previous extensive experience in the numerical simulation of 1D, steady, multiphase flow dynamics along volcanic conduits and fissures, and combine it with a lumped description for the evolution of magma chamber conditions, within a consistent approach detailed in the Methods section. We remark that adding a magmastatic pressure profile description within the magma chamber, as in previous work [e.g., Bower and Woods (1998); Colucci et al. (2014)], would not be a major complication. However, Colucci et al. (2014) find that pressure differences as large as $\text{>}100\text{MPa}$ from chamber top to bottom do not play major roles in the eruption dynamics (in a one chamber plus conduit system). Furthermore, that would require additional assumptions on chamber geometry to relate chamber volume to height, or additional sensitivity analysis to explore possible chamber geometries. Finally, during the highly dynamic eruption phases dynamic pressure contributions, mechanical dissipation due to friction, generation of vortexes, etc. are expected to largely disrupt the magmastatic distribution, which would thus be an unjustified assumption. On the contrary, the magma chamber lumped approach adopted here does not introduce any further assumption while accounting for substantial complexities, such as the non-linear relationships between pressure and gas-melt partitioning of volatiles upon mixing between inflowing and outflowing magmas having different composition, volatile contents, and properties. We stress, again, the meaning of dyke thickness in our approach, which is that of a proxy for the efficiency of deep, inter-chamber magma transfer through a system which is expected to be complex, e.g., multiple, possibly intersecting dykes crossing intermediate sills or minor magma storage regions. Accordingly, in the following discussion we mostly refer to dyke thickness as the “deep magma transfer efficiency”.

The primary message from this study is the existence of eruption regimes having different evolution and termination, mostly controlled by the deep magma transfer efficiency, and resulting in a major control on the eruption evolution by either the shallow or the deep region of the volcanic plumbing system (Figure 8). Factors other than the deep magma transfer efficiency have an influence in determining the onset of the shallow or deep-controlled eruption regime, and among them, the size of the deep magmatic reservoir is important (Figure 3). However, no other quantity emerges so strikingly as the deep magma transfer efficiency. Its central role in the eruption dynamics is that of balancing fast magma withdrawal from the shallow magmatic system during explosive eruptions. If the mass flow-rate across reservoirs does not compensate for the eruption mass flow-rate, then the pressure is destined to decrease in the shallow magmatic system (Figure 7) (in turn increasing the deep mass flow-rate, thus triggering non-linear effects accounted for in the present modeling). The eruption intensity progressively declines until the difference between lithostatic and magmatic pressure exceeds the rock tensile strength and shallow rock failure occurs, after a minor amount of deep magma has been erupted together with the shallow, more evolved magma making up nearly the totality of the volcanic products (Figure 9). If, instead, the deep magma transfer efficiency over-compensates for the eruption mass flow-rate, then the shallow system undergoes an initial phase of pressurization and quick increase of the eruption intensity (that is, increase of mass flow-rate and volcanic plume height), before entering a long phase where the eruption mass flow-rate (and plume height) can oscillate, remain constant, or slowly decrease (Figure 7), until a sufficient amount of deep magma is extracted from the deep system, and a portion of it is erupted to the surface (Figure 9), causing sufficient pressure decrease in the deep magmatic system and triggering unstable conditions at such deep level.

FIGURE 9

FIGURE 9. Calculated evolution of the average composition of the erupted magma, for the same simulation cases as in Figure 7.

4.2 Comparison With Real Eruptions

The overall picture summarized above is confirmed by the close correspondence between the intensity-magnitude relationship computed in this study by just varying the simulation conditions over a-priori defined, reasonable ranges from the literature (Table 1), and that from observed or reconstructed eruptions worldwide (Figure 1). The eruption magnitudes and durations resulting from the two regimes can largely overlap (Figures 1, 8), making the attribution of real eruptions to one or the other regime not straightforward. However, in the analyzed range of conditions the largest magnitudes and durations, approaching those of super-eruptions, are found for the deep-controlled regime (Figures 1, 4, 8). Similarly, if the size of the collapsing chamber is taken as a proxy for the width of a possibly associated caldera, caldera sizes can largely overlap between the two regimes, but the largest ones, of order tens of km, are found in association to the largest simulated eruption magnitudes in the deep-controlled regime (Figure 3). It seems plausible, therefore, that the known largest volcanic super-eruptions developed under conditions corresponding with those along the red curves in Figures 7–9: if intermediate magma accumulation zones existed between the large, deeper reservoir and the surface, the efficiency of the deep magma interconnections was however large enough to avoid quick shallow pressure decrease and rapid eruption closure. It may be speculated that if, instead, the deep magma transfer efficiency was not sufficiently large, then the volcanic system might have entered a phase of frequent low magnitude eruptions, each one terminating as a consequence of shallow system control, and each one quickly consuming part of the deep chamber overpressure. In other words, the simulation results seem to suggest that huge accumulation of eruptible magma at some depth below a volcano is necessary but not sufficient for the occurrence of a very large, globally impacting eruption. Such an extreme event appears to equally require the establishment of unusual conditions allowing high rates of magma transfer at depth, across an effective, highly efficient network of interconnections linking the large reservoir to its shallow, smaller hypophysis closer to the Earth surface. Figure 8 shows that for the conditions of the simulations in the figure, the minimum dyke thickness for deep-controlled eruption regime nearly doubles when the deep chamber volume changes from 200 to 20,000 km³. The establishment of such unusual conditions may in turn provide a justification to the rareness of volcanic super-eruptions [VEI eight eruptions represent only one-two out of 10 million volcanic eruptions on Earth, (Papale. 2018)], independent from the existence of sufficiently large magmatic reservoirs, that may exist but not be able to produce a huge eruption. Geophysical imaging suggests the presence of huge amounts of magma beneath an increasing number of volcanoes: at Yellowstone (WY, United States) a shallow (upper crust) rhyolitic magma storage zone is estimated to have a volume of about $10\text{,}000{\text{km}}^3$, plus additional $46\text{,}000{\text{km}}^3$ of deeper, more mafic magma (Husen et al., 2004; Reuber et al., 2018); at Campi Flegrei (Italy) a magmatic source with composition variable with depth has an estimated volume of $40\text{,}000{\text{km}}^3$ and it includes 30% of molten material (Fedi et al., 2018); and in a number of other cases, including Mount St. Heles, WA (Waite and Moran, 2009) and Nevado del Ruiz, Colombia (Londoño and Kumagai, 2018) large magma accumulations are highlighted down to the largest depths (7–10 km) allowed by the required resolution. These and other similar cases led previous investigators to wonder how can magma accumulate to the amounts required for the occurrence of a volcanic super-eruption, without being discharged to the surface [e.g., Miller and Wark (2008)]. The results of this study offer a possible explanation: without the establishment of a system of interconnections allowing extremely efficient transfer of magma from deep reservoirs, no matter how large, to their shallow counterparts, magmatic giants cannot give rise to volcanic giants. With reference to the $20\text{,}000{\text{km}}^3$ deep chamber case in Figure 8, poorly efficient interconnections would correspond instead to eruptions over the left end portion of the diagram, characterized by the occurrence of medium size eruptions controlled by the shallow magmatic system. The situation at Campi Flegrei looks exemplary: here the two major 39 ka Campanian Ignimbrite and 15 ka Yellow Tuff eruption had a magnitude of 7.7–7.8 and 6.8, respectively (Albert et al., 2019; Silleni et al., 2020), contributing to shape a large caldera over a minimum area of 100 ${\text{km}}^2$. Although a minimum volume of 10,000 ${\text{km}}^3$ of trachy-andesitic to trachytic melt is likely available to-date in a huge reservoir below 8 km depth (Astbury et al., 2018; Fedi et al., 2018), tens of explosive eruptions during last 15 ka have produced volcanic materials each in the range $\text{<}100{\text{Mm}}^3-3{\text{km}}^3$ (De Vita et al., 1999; Orsi et al., 2009) (magnitude in the range <4–5.5), resulting in a constellation of craters and smaller calderas with area typically less than $4{\text{km}}^2$ which punctuate the large caldera floor. Such small to medium eruptions were fed by chemically evolved magma that resided at shallow, 2–4 km depth, and that is invariably mixed up with minor amounts of less chemically evolved magma with deeper provenance (Arienzo et al., 2010; Di Renzo et al., 2011; Arienzo et al., 2015; Astbury et al., 2018). Although the magmas erupted at Campi Flegrei belong to a different compositional suite with respect to the andesitic suite investigated here, and that characterizes plate convergence regions where most active volcanoes concentrate, the general trends are expected to be qualitatively similar. Thus, although the amount of melt sufficient to feed a globally impacting VEI seven or eight eruption may be available at Campi Flegrei, it appears that the system of deep interconnections during last 15 ka has not been sufficiently developed to sustain large mass flow rates, resulting in shallow-controlled conditions such as those toward the left end side in the diagrams of Figure 8, and leading to eruption closure by rock failure at shallow level and to the myriad of small craters and calderas within the large caldera, remnant of eruptions that were instead controlled from depth.

4.3 Implications for Volcanic Hazard Forecasts

The present study suggests that the magnitude, thus the impacts, of a volcanic eruption largely depends on the dyke thickness parameter, a proxy for the efficiency of deep magmatic interconnections. That is an unfortunate situation for our forecasting capabilities. In fact, we do not have, currently, any viable means of estimating such an elusive quantity at any real volcano. As an example, what conditions should we assume in case of escalating unrest at Yellowstone, which has never erupted during last 70 ka? And, would it be safe to assume that the next eruption at Campi Flegrei, where the last eruption about 500 years ago was preceded by about 3.4 ka of quiescence, would still reflect low efficient deep magma interconnections as for the previous millennia? While we have progressed enormously in revealing even the smallest unrest signals at well monitored volcanoes, and we expect therefore to recognize the early signs of escalating dynamics, our capability to anticipate the size of an impending eruption has practically not evolved during last decades. Considerations based on the global size distribution of volcanic eruptions on Earth suggest that such an uncomfortable situation may result from exceedingly non-linear eruption dynamics, translating into the observed power law distribution of the frequency of dominantly explosive eruptions (Papale, 2018; Papale et al., 2021). This study shows that modeling the eruption dynamics leads to the similar conclusions that finding confident relationships between unrest dynamics and size of the impending eruption may continue to be, at least for explosive eruptions, a very hard task, if not a hopeless undertaking. The situation may be more favorable after the eruption has begun. In fact, conditions of high deep magma transfer efficiency can be associated with an initial phase of increasing eruption mass flow-rate, as it is shown in Figure 7 (and in the Supplementary Figure S6). For the conditions in such figures, and assuming no mass loss from the volcanic plume due to partial or total column collapses, the volcanic column would be seen to increase by about 2–>10 km (Costa et al., 2018) over the first one to several hours from the establishment of sub-steady conditions. Therefore, an initially progressively increasing (sub-steady) column height may be suggestive of deep-controlled eruption regime and associated high to very high eruption duration and magnitude, with the potential to evolve into a large caldera-forming eruption. A number of well-studied Plinian fallout deposits from the largest caldera-forming eruptions on Earth exhibit basal reverse grading; e.g., the Plinian fallout preceding the Campanian Ignimbrite eruption, mentioned above, shows basal reverse grading for more than its half height (Rosi et al., 1999); and the Plinian phase of the 760 ka M8.3 Bishop Tuff eruption, Long Valley, United States, shows reverse grading throughout its entirety, reflecting nearly continuous eruption column height thus mass flow-rate increase before pyroclastic flows leading to the gigantic ignimbrite formation (Gardner et al., 1991). In practice, however, referring to time-wize mass flow-rate increase as a syn-eruptive early sign of deep-controlled eruption regime thus of potentially more dangerous, more impacting conditions, can be difficult. In fact, mass flow-rate increase during initial eruption stages may reflect other factors and processes, like progressive wall rock erosion and establishment of the volcanic conduit (Aravena et al., 2017), progressive tapping of more water-rich, less chemically evolved (thus less viscous) magma (Papale et al., 1998), and others that can affect to various degrees the stages preceding and accompanying full development of a Plinian volcanic column. Many of these processes can be identified and reconstructed from the analysis of the associated volcanic deposits [e.g., Cioni et al. (2021)], but are difficult to constrain while the eruption is on-going. With the above limitations, in light of the present results we suggest that for those cases where the eruption mass flow-rate is seen to increase over hours, then it may be wise to consider an increased likelihood of established deep-controlled regime conditions, as well as the possibility of a long eruption duration of order days to weeks with comparably large accumulation of volcanic ash, potentially ending up with large caldera collapse and generation of large pyroclastic flows.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

SC and PP conceived the present idea and the algorithm, discussed the results, and wrote the manuscript. SC implemented the algorithm and performed the numerical simulations.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2021.681083/full#supplementary-material

Appendix

The fundamental equations describing transport of mass and momentum along a conduit or dyke assume 1D, steady, isothermal, multiphase conditions (Papale, 2001). The mass conservation for the gas phase (bubbles below fragmentation, continuous above it) is given by

$\frac{d}{dz}\left(\epsilon {\rho }_G{u}_G\right)={\dot{m}}_A\frac{d{w}_G}{dz},(9)$

and the equation for the dense phase (liquid + crystals below fragmentation, particles/droplets above it) by

$\frac{d}{dz}\left[\left(1-\epsilon \right){\rho }_D{u}_D\right]=-{\dot{m}}_A\frac{d{w}_G}{dz}.(10)$

The next two equations are for momentum balance. The equation for the gas phase is

$\epsilon {\rho }_G{u}_G\frac{d{u}_G}{dz}=-\epsilon \frac{dP}{dz}-\epsilon {\rho }_{G}g-F_{wG}$

$-\xi \left(u_G-u_D\right)-\delta \left(u_G-u_D\right){\dot{m}}_A\frac{d{w}_G}{dz},(11)$

and the equation for the dense phase is

$\left(1-\epsilon \right){\rho }_D{u}_D\frac{d{u}_D}{dz}=-\left(1-\epsilon \right)\frac{dP}{dz}$

$-\left(1-\epsilon \right){\rho }_{D}g-F_{wD}+\xi \left(u_G-u_D\right)$

$-\left(1-\delta \right)\left(u_G-u_D\right){\dot{m}}_A\frac{d{w}_G}{dz}-\Pi \frac{d{\sigma }_s}{dz}.(12)$

In the above equations, ε is gas volume fraction, ρ is density, u is velocity, w is mass fraction, z is the vertical coordinate upward positive, P is pressure, g is gravity acceleration, $F_{wG}$ and $F_{wD}$ are gas-wall and dense phase-wall friction terms, respectively, ξ is the interphase drag coefficient, δ is the coefficient of interphase momentum transfer due to mass transfer, ${\dot{m}}_A$ is the mass flow rate per unit area, ${\sigma }_s$ is the normal stress due to particle-particle collisions ($\Pi =0$ for bubbly flow, $\Pi =1$ for gas-particle flow), and the subscripts G and D refer to the gas and dense phase, respectively.

The interphase mass transfer term $d{w}_G/dz$ in the above equations is given by

$\frac{d{w}_G}{dz}=-\frac{\left(1-w_{TC}\right)\left(1-w_V^o\right)}{{\left(1-w_{DV}\right)}^2}\frac{d{w}_{DV}}{dP}\frac{dP}{dz},(13)$

where $w_{TC}$ is the total mass fraction of crystals, $w_V^o$ is the total volatile content (water and carbon dioxide in the present modeling) with respect to the crystal-free magma, and $w_{DV}$ is the mass fraction of dissolved volatiles with respect to the crystal-free magma. In turn, the quantities $w_{DV}$ and $d{w}_{DV}/dP$ in Eq. 13 are computed from the real multi-component melt-gas equilibrium model in Papale et al. (2006).

Magma fragmentation is determined according to local conditions in terms of brittle/ductile response of the magmatic mixture. Fragmentation occurs when the following criterion is satisfied (Papale, 1999):

$\frac{d{u}_z}{dz}>k\frac{G_{\infty }}{{\mu }_m},(14)$

where $k=0.01$ and $G_{\infty }$ is magma elastic modulus at infinite frequency, and ${\mu }_m$ is the viscosity of the multiphase magma.

The above equations are the same for conduit and dyke flow. However, the two cases differ for the friction term, which expresses in 1D the effective mechanical energy dissipation due to the internal friction related to the effective velocity distribution in the plane perpendicular to flow direction. The friction term $F_{wG}$ and $F_{wD}$ at Eqs 11, 12 for the bubbly flow regime (below fragmentation, subscript BF) and gas-particle flow regime (above fragmentation, subscript GPF) are given by:

$F_{wG,BF}=0,(15)$

$F_{wD,BF},F_{wG,GPF}=\frac{2f{\dot{m}}_A^2}{D_H{\rho }_m},(16)$

$F_{wD,GPF}=\frac{2f_s\left(1-\epsilon \right){\rho }_D{u}_D^2}{D_H},(17)$

$f=\frac{B_1}{\mathrm{Re}}+B_2,\mathrm{Re}=\frac{{\dot{m}}_A{D}_H}{{\mu }_m},(18)$

where $D_H$ is the hydraulic diameter given by

$D_H=\frac{4A}{\chi },(19)$

where A is the cross-sectional area and χ is the wetted perimeter. For a cylindrical conduit and a planar dyke Eq. 19 gives

$\text{conduit}:D_H=D,(20)$

$\text{dyke}:D_H=2W,(21)$

where D is conduit diameter and W is dyke thickness.

The quantity f in Eqs 18 is the friction coefficient, which embeds the internal dissipation associated with the effective velocity profile. For a Newtonian flow and with reference to the first of Eq. 18, the following holds:

$\text{conduit}:B_1=16,(22)$

$\text{dyke}:B_1=24,(23)$

B₂$B_2$ takes the value of 0.01 in both cases. Finally, $f_s$ at Eq. 17 is particle friction factor given by:

$f_s=\frac{0.0285{\left(g{d}_D\right)}^{1/2}}{u_D},(24)$

where $d_D$ is the average diameter of particles flowing above fragmentation.

Other equations describing the additional terms in Eqs 9–12 are reported in Papale (2001).

References

Adams, B., Ebeida, M. S., Eldred, M. S., Jakeman, J. D., Swiler, L. P., Stephens, J. A., et al. (2014). Dakota, a Multilevel Parallel Object-Oriented Framework for Design Optimization, Parameter Estimation, Uncertainty Quantification, and Sensitivity Analysis: Version 6.0 User’s Manual. New Mexico: Sandia Technical Report SAND, 2014–4633. doi:10.2172/1177048

CrossRef Full Text

Albert, P. G., Giaccio, B., Isaia, R., Costa, A., Niespolo, E. M., Nomade, S., et al. (2019). Evidence for a Large-Magnitude Eruption from Campi Flegrei Caldera (Italy) at 29 Ka. Geology 47, 595–599. doi:10.1130/G45805.1

CrossRef Full Text | Google Scholar

Aravena, Á., Cioni, R., de’Michieli Vitturi, M., and Neri, A. (2018). Conduit Stability Effects on Intensity and Steadiness of Explosive Eruptions. Sci. Rep. 8, 4125. doi:10.1038/s41598-018-22539-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Aravena, Á., de' Michieli Vitturi, M., Cioni, R., and Neri, A. (2017). Stability of Volcanic Conduits during Explosive Eruptions. J. Volcanology Geothermal Res. 339, 52 – 62. doi:10.1016/j.jvolgeores.2017.05.003

CrossRef Full Text | Google Scholar

Araya, N., Nakamura, M., Yasuda, A., Okumura, S., Sato, T., Iguchi, M., et al. (2019). Shallow Magma Pre-charge during Repeated Plinian Eruptions at Sakurajima Volcano. Sci. Rep. 9, 1979. doi:10.1038/s41598-019-38494-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Arienzo, I., D’Antonio, M., Di Renzo, V., Tonarini, S., Minolfi, G., Orsi, G., et al. (2015). Isotopic Microanalysis Sheds Light on the Magmatic Endmembers Feeding Volcanic Eruptions: the Astroni 6 Case Study (Campi Flegrei, italy). J. Volcanology Geothermal Res. 304, 24–37. doi:10.1016/j.jvolgeores.2015.08.003

CrossRef Full Text | Google Scholar

Arienzo, I., Moretti, R., Civetta, L., Orsi, G., and Papale, P. (2010). The Feeding System of Agnano-Monte Spina Eruption (Campi Flegrei, Italy): Dragging the Past into Present Activity and Future Scenarios. Chem. Geology. 270, 135–147. doi:10.1016/j.chemgeo.2009.11.012

CrossRef Full Text | Google Scholar

Astbury, R. L., Petrelli, M., Ubide, T., Stock, M. J., Arienzo, I., D'Antonio, M., et al. (2018). Tracking Plumbing System Dynamics at the Campi Flegrei Caldera, italy: High-Resolution Trace Element Mapping of the Astroni crystal Cargo. Lithos 318-319, 464–477. doi:10.1016/j.lithos.2018.08.033

CrossRef Full Text | Google Scholar

Bonaccorso, A., Calvari, S., Linde, A., and Sacks, S. (2014). Eruptive Processes Leading to the Most Explosive Lava fountain at etna Volcano: The 23 November 2013 Episode. Geophys. Res. Lett. 41, 4912–4919. doi:10.1002/2014GL060623

CrossRef Full Text | Google Scholar

Bower, S. M., and Woods, A. W. (1998). On the Influence of Magma chambers in Controlling the Evolution of Explosive Volcanic Eruptions. J. Volcanology Geothermal Res. 86, 67–78. doi:10.1016/s0377-0273(98)00081-x

CrossRef Full Text | Google Scholar

Brown, S. K., Crosweller, H. S., Sparks, R. S. J., Cottrell, E., Deligne, N. I., Guerrero, N. O., et al. (2014). Characterisation of the Quaternary Eruption Record: Analysis of the Large Magnitude Explosive Volcanic Eruptions (LaMEVE) Database. J. Appl. Volcanol. 3. doi:10.1186/2191-5040-3-5

CrossRef Full Text | Google Scholar

Cioni, R., Pistolesi, M., Pompilio, M., and Scaillet, B. (2021). Understanding Volcanic Systems and Their Dynamics Combining Field and Physical Volcanology with Petrology studiesForecasting and Planning for Volcanic Hazards, Risks, and Disasters. Elsevier, 285–328. doi:10.1016/b978-0-12-818082-2.00007-x

CrossRef Full Text | Google Scholar

Colucci, S., de' Michieli Vitturi, M., Neri, A., and Palladino, D. M. (2014). An Integrated Model of Magma Chamber, Conduit and Column for the Analysis of Sustained Explosive Eruptions. Earth Planet. Sci. Lett. 404, 98 – 110. doi:10.1016/j.epsl.2014.07.034

CrossRef Full Text | Google Scholar

Costa, A., Suzuki, Y. J., and Koyaguchi, T. (2018). Understanding the Plume Dynamics of Explosive Super-eruptions. Nat. Commun. 9, 1–6. doi:10.1038/s41467-018-05545-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Crosweller, H. S., Arora, B., Brown, S. K., Cottrell, E., Deligne, N. I., Guerrero, N. O., et al. (2012). Global Database on Large Magnitude Explosive Volcanic Eruptions (Lameve). J. Appl. Volcanol. 1. doi:10.1186/2191-5040-1-4

CrossRef Full Text | Google Scholar

De Vita, S., Orsi, G., Civetta, L., Carandente, A., D'Antonio, M., Deino, A., et al. (1999). The Agnano-Monte Spina Eruption (4100 Years BP) in the Restless Campi Flegrei Caldera (Italy). J. Volcanology Geothermal Res. 91, 269–301. doi:10.1016/s0377-0273(99)00039-6

CrossRef Full Text | Google Scholar

Di Renzo, V., Arienzo, I., Civetta, L., D'Antonio, M., Tonarini, S., Di Vito, M. A., et al. (2011). The Magmatic Feeding System of the Campi Flegrei Caldera: Architecture and Temporal Evolution. Chem. Geology 281, 227–241. doi:10.1016/j.chemgeo.2010.12.010

CrossRef Full Text | Google Scholar

Edmonds, M., Cashman, K. V., Holness, M., and Jackson, M. (2019). Architecture and Dynamics of Magma Reservoirs. Phil. Trans. R. Soc. A. 377, 20180298. doi:10.1098/rsta.2018.0298

PubMed Abstract | CrossRef Full Text | Google Scholar

Fedi, M., Cella, F., D’Antonio, M., Florio, G., Paoletti, V., and Morra, V. (2018). Gravity Modeling Finds a Large Magma Body in the Deep Crust below the Gulf of naples, italy. Sci. Rep. 8. doi:10.1038/s41598-018-26346-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Folch, A., and Martí, J. (2009). Time-dependent Chamber and Vent Conditions during Explosive Caldera-Forming Eruptions. Earth Planet. Sci. Lett. 280, 246 – 253. doi:10.1016/j.epsl.2009.01.035

CrossRef Full Text | Google Scholar

Folch, A., and Martı́, J. (1998). The Generation of Overpressure in Felsic Magma chambers by Replenishment. Earth Planet. Sci. Lett. 163, 301 – 314. doi:10.1016/S0012-821X(98)00196-4

CrossRef Full Text | Google Scholar

Gardner, J. E., Sigurdsson, H., and Carey, S. N. (1991). Eruption Dynamics and Magma Withdrawal during the Plinian Phase of the bishop Tuff Eruption, long valley Caldera. J. Geophys. Res. 96, 8097–8111. doi:10.1029/91jb00257

CrossRef Full Text | Google Scholar

Gudmundsson, A. (2007). Conceptual and Numerical Models of Ring-Fault Formation. J. Volcanology Geothermal Res. 164, 142 – 160. doi:10.1016/j.jvolgeores.2007.04.018

CrossRef Full Text | Google Scholar

Gudnason, J., Thordarson, T., Houghton, B. F., and Larsen, G. (2018). The 1845 Hekla Eruption: Grain-Size Characteristics of a Tephra Layer. J. Volcanology Geothermal Res. 350, 33 – 46. doi:10.1016/j.jvolgeores.2017.11.025

CrossRef Full Text | Google Scholar

Husen, S., Smith, R. B., and Waite, G. P. (2004). Evidence for Gas and Magmatic Sources beneath the Yellowstone Volcanic Field from Seismic Tomographic Imaging. J. Volcanology Geothermal Res. 131, 397–410. doi:10.1016/S0377-0273(03)00416-5

CrossRef Full Text | Google Scholar

Kabele, P., Žák, J., and Somr, M. (2017). Finite-element Modeling of Magma Chamber-Host Rock Interactions Prior to Caldera Collapse. Geophys. J. Int. 209, 1851–1865. doi:10.1093/gji/ggx121

CrossRef Full Text | Google Scholar

Karlstrom, L., Rudolph, M. L., and Manga, M. (2012). Caldera Size Modulated by the Yield Stress within a crystal-rich Magma Reservoir. Nat. Geosci 5, 402–405. doi:10.1038/ngeo1453

CrossRef Full Text | Google Scholar

Keiding, J. K., and Sigmarsson, O. (2012). Geothermobarometry of the 2010 Eyjafjallajökull Eruption: New Constraints on Icelandic Magma Plumbing Systems. J. Geophys. Res. 117. doi:10.1029/2011JB008829

CrossRef Full Text | Google Scholar

Londoño, J. M., and Kumagai, H. (2018). 4D seismic tomography of Nevado del Ruiz Volcano, Colombia, 2000-2016. J. Volcanology Geothermal Res. 358, 105–123. doi:10.1016/j.jvolgeores.2018.02.015

CrossRef Full Text | Google Scholar

Macedonio, G., Dobran, F., and Neri, A. (1994). Erosion Processes in Volcanic Conduits and Application to the Ad 79 Eruption of Vesuvius. Earth Planet. Sci. Lett. 121, 137 – 152. doi:10.1016/0012-821X(94)90037-X

CrossRef Full Text | Google Scholar

Macedonio, G., Neri, A., Martì, J., and Folch, A. (2005). Temporal Evolution of Flow Conditions in Sustained Magmatic Explosive Eruptions. J. Volcanology Geothermal Res. 143, 153 – 172. doi:10.1016/j.jvolgeores.2004.09.015

CrossRef Full Text | Google Scholar

Mandeville, C. W., Carey, S., and Sigurdsson, H. (1996). Magma Mixing, Fractional Crystallization and Volatile Degassing during the 1883 Eruption of Krakatau Volcano, indonesia. J. Volcanology Geothermal Res. 74, 243–274. doi:10.1016/S0377-0273(96)00060-1

CrossRef Full Text | Google Scholar

Marianelli, P., Sbrana, A., and Proto, M. (2006). Magma Chamber of the Campi Flegrei Supervolcano at the Time of Eruption of the Campanian Ignimbrite. Geol 34, 937–940. doi:10.1130/G22807A.1

CrossRef Full Text | Google Scholar

Martí, J., Folch, A., Neri, A., and Macedonio, G. (2000). Pressure Evolution during Explosive Caldera-Forming Eruptions. Earth Planet. Sci. Lett. 175, 275 – 287. doi:10.1016/S0012-821X(99)00296-4

CrossRef Full Text | Google Scholar

Melnik, O., and Costa, A. (2014). Chapter 3 Dual-Chamber-Conduit Models of Non-linear Dynamics Behaviour at Soufrière Hills Volcano, Montserrat. Geol. Soc. Lond. Mem. 39, 61–69. doi:10.1144/m39.3

CrossRef Full Text | Google Scholar

Miller, C. F., and Wark, D. A. (2008). SUPERVOLCANOES AND THEIR EXPLOSIVE SUPERERUPTIONS. Elements 4, 11–15. doi:10.2113/GSELEMENTS.4.1.11

CrossRef Full Text | Google Scholar

Miller, T. P., and Chouet, B. A. (1994). The 1989-1990 Eruptions of Redoubt Volcano: an Introduction. J. Volcanology Geothermal Res. 62, 1 – 10. doi:10.1016/0377-0273(94)90025-6

CrossRef Full Text | Google Scholar

Murphy, M. D., Sparks, R. S. J., Barclay, J., Carroll, M. R., and Brewer, T. S. (2000). Remobilization of Andesite Magma by Intrusion of Mafic Magma at the Soufriere hills Volcano, montserrat, West Indies. J. Petrology 41, 21–42. doi:10.1093/petrology/41.1.21

CrossRef Full Text | Google Scholar

Neri, A., Papale, P., and Macedonio, G. (1998). The Role of Magma Composition and Water Content in Explosive Eruptions: 2. Pyroclastic Dispersion Dynamics. J. Volcanology Geothermal Res. 87, 95 – 115. doi:10.1016/S0377-0273(98)00102-4

CrossRef Full Text | Google Scholar

Orsi, G., Di Vito, M. A., Selva, J., and Marzocchi, W. (2009). Long-term Forecast of Eruption Style and Size at Campi Flegrei Caldera (italy). Earth Planet. Sci. Lett. 287, 265–276. doi:10.1016/j.epsl.2009.08.013

CrossRef Full Text | Google Scholar

Pallister, J. S., Hoblitt, R. P., and Reyes, A. G. (1992). A basalt Trigger for the 1991 Eruptions of Pinatubo Volcano? Nature 356, 426–428. doi:10.1038/356426a0

CrossRef Full Text | Google Scholar

Papale, P., and Dobran, F. (1994). Magma Flow along the Volcanic Conduit during the Plinian and Pyroclastic Flow Phases of the May 18, 1980, Mount St. Helens Eruption. J. Geophys. Res. 99, 4355–4373. doi:10.1029/93JB02972

CrossRef Full Text | Google Scholar

Papale, P. (2001). Dynamics of Magma Flow in Volcanic Conduits with Variable Fragmentation Efficiency and Nonequilibrium Pumice Degassing. J. Geophys. Res. 106, 11043–11065. doi:10.1029/2000JB900428

CrossRef Full Text | Google Scholar

Papale, P. (2018). Global Time-Size Distribution of Volcanic Eruptions on Earth. Sci. Rep. 8, 6838. doi:10.1038/s41598-018-25286-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Papale, P., Marzocchi, W., and Garg, D. (2021). Global Volume Distribution for continental Volcanism on Earth. Under Review. J. Geophys. Res. doi:10.1029/2021jb021763

CrossRef Full Text | Google Scholar

Papale, P., Moretti, R., and Barbato, D. (2006). The Compositional Dependence of the Saturation Surface of H2o+co2 Fluids in Silicate Melts. Chem. Geology 229, 78 – 95. doi:10.1016/j.chemgeo.2006.01.013

CrossRef Full Text | Google Scholar

Papale, P., Neri, A., and Macedonio, G. (1998). The Role of Magma Composition and Water Content in Explosive Eruptions. J. Volcanology Geothermal Res. 87, 75–93. doi:10.1016/s0377-0273(98)00101-2

CrossRef Full Text | Google Scholar

Papale, P. (1999). Strain-induced Magma Fragmentation in Explosive Eruptions. Nature 397, 425–428. doi:10.1038/17109

PubMed Abstract | CrossRef Full Text | Google Scholar

Pianosi, F., Beven, K., Freer, J., Hall, J. W., Rougier, J., Stephenson, D. B., et al. (2016). Sensitivity Analysis of Environmental Models: A Systematic Review with Practical Workflow. Environ. Model. Softw. 79, 214–232. doi:10.1016/j.envsoft.2016.02.008

CrossRef Full Text | Google Scholar

Pianosi, F., and Wagener, T. (2015). A Simple and Efficient Method for Global Sensitivity Analysis Based on Cumulative Distribution Functions. Environ. Model. Softw. 67, 1–11. doi:10.1016/j.envsoft.2015.01.004

CrossRef Full Text | Google Scholar

Pianosi, F., and Wagener, T. (2018). Distribution-based Sensitivity Analysis from a Generic Input-Output Sample. Environ. Model. Softw. 108, 197–207. doi:10.1016/j.envsoft.2018.07.019

CrossRef Full Text | Google Scholar

Pioli, L., Erlund, E., Johnson, E., Cashman, K., Wallace, P., Rosi, M., et al. (2008). Explosive Dynamics of Violent Strombolian Eruptions: The Eruption of Parícutin Volcano 1943-1952 (Mexico). Earth Planet. Sci. Lett. 271, 359 – 368. doi:10.1016/j.epsl.2008.04.026

CrossRef Full Text | Google Scholar

Pistolesi, M., Isaia, R., Marianelli, P., Bertagnini, A., Fourmentraux, C., Albert, P. G., et al. (2016). Simultaneous Eruptions from Multiple Vents at Campi Flegrei (Italy) Highlight New Eruption Processes at Calderas. Geology 44, 487–490. doi:10.1130/G37870.1

CrossRef Full Text | Google Scholar

Pyle, D. M. (2015). Sizes of Volcanic Eruptions. In The Encyclopedia of Volcanoes (Second Edition), ed. H. Sigurdsson (Amsterdam: Academic Press). Second edition edn., 257 – 264. doi:10.1016/B978-0-12-385938-9.00013-4

CrossRef Full Text | Google Scholar

Reuber, G. S., Kaus, B. J. P., Popov, A. A., and Baumann, T. S. (2018). Unraveling the Physics of the Yellowstone Magmatic System Using Geodynamic Simulations. Front. Earth Sci. 6, 117. doi:10.3389/feart.2018.00117

CrossRef Full Text | Google Scholar

Reverso, T., Vandemeulebrouck, J., Jouanne, F., Pinel, V., Villemin, T., Sturkell, E., et al. (2014). A Two-Magma Chamber Model as a Source of Deformation at Grímsvötn Volcano, Iceland. J. Geophys. Res. Solid Earth 119, 4666–4683. doi:10.1002/2013jb010569

CrossRef Full Text | Google Scholar

Saltelli, A., Aleksankina, K., Becker, W., Fennell, P., Ferretti, F., Holst, N., et al. (2019). Why So many Published Sensitivity Analyses Are False: A Systematic Review of Sensitivity Analysis Practices. Environ. Model. Softw. 114, 29–39. doi:10.1016/j.envsoft.2019.01.012

CrossRef Full Text | Google Scholar

Samaniego, P., Pennec, L., Robin, L., and Hidalgo, S. (2011). Magma Mixing: a Mechanism for Triggering Acid Explosive Eruption. Nature 257, 315–318. doi:10.1038/267315a0

CrossRef Full Text | Google Scholar

Silleni, A., Giordano, G., Isaia, R., and Ort, M. (2020). The Magnitude of the 39.8 Ka Campanian Ignimbrite Eruption, italy: Method, Uncertainties and Errors. Front. Earth Sci.

Google Scholar

Sparks, S. R. J., Sigurdsson, H., and Wilson, L. (1977). Magma Mixing: a Mechanism for Triggering Acid Explosive Eruptions. Nature 267, 315–318. doi:10.1038/267315a0

CrossRef Full Text | Google Scholar

Stix, J., Layne, G. D., and Williams, S. N. (2003). Mechanisms of degassing at nevado del ruiz volcano, colombia. J. Geol. Soc. 160, 507–521. doi:10.1144/0016-764902-028

CrossRef Full Text | Google Scholar

Tepley, F. J., Tilling, R. I., Arth, J. G., and Arth, J. G. (2000). Magma Mixing, Recharge and Eruption Histories Recorded in Plagioclase Phenocrysts from El Chichón Volcano, Mexico. J. Petrology 41, 1397–1411. doi:10.1093/petrology/41.9.1397

CrossRef Full Text | Google Scholar

Venezky, D. Y., and Rutherford, M. J. (1997). Preeruption Conditions and Timing of Dacite-Andesite Magma Mixing in the 2.2 Ka Eruption at mount rainier. J. Geophys. Res. 102, 20069–20086. doi:10.1029/97JB01590

CrossRef Full Text | Google Scholar

Waite, G. P., and Moran, S. C. (2009). VP Structure of Mount St. Helens, Washington, USA, Imaged with Local Earthquake Tomography. J. Volcanology Geothermal Res. 182, 113–122. doi:10.1016/j.jvolgeores.2009.02.009

CrossRef Full Text | Google Scholar

Wilson, L., Sparks, R. S. J., and Walker, G. P. L. (1980). Explosive Volcanic Eruptions -- IV. The Control of Magma Properties and Conduit Geometry on Eruption Column Behaviour. Geophys. J. Int. 63, 117–148. doi:10.1111/j.1365-246X.1980.tb02613.x

CrossRef Full Text | Google Scholar