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 imply a limit to eruption size forecasts based on observations and measurements during volcanic unrest.


Introduction
Explosive eruptions are characterised by a variety of mass flow rate evolutions ([22] and references therein) and erupted mass spanning several orders of magnitude [23]. The activation of multiple reservoirs during individual eruptions is extensively documented (e.g., [2] and references therein) and accepted as a general reference model [8]. Still, dynamic models of volcanic eruptions consider at most one single, progressively emptied magma reservoir [3,5,9,10,14,15]. Such models strictly associate larger magma chambers to larger eruption intensity and longer eruption duration, and invariably predict progressive chamber pressure and mass flow-rate decrease with time. On the contrary, real eruptions show virtually any mass flow-rate (or intensity) evolution with time, requiring therefore a different approach to explain such first order observations. By considering multiple magma chambers mutually interacting during volcanic eruptions, an entirely new world appears, explaining the observed eruption magnitudes and intensities and the large variability of mass flow rate evolutions. A new fundamental quantity emerges, represented by the efficiency of magma transfer across the magmatic reservoirs. That quantity largely controls eruption evolution and the causes of its termination, and ultimately the impacts on the surrounding environment, with substantial implications for volcanic hazards.
1 Numerical simulations Figure 1 (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 width, volatile abundance in the two magmas) are varied in the range reported in the figure caption, reflecting more frequently reported values for explosive eruptions in a wide literature. 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 [17]. 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 [6,14,18], with boundary conditions at each step depending on the previous dynamics and reflecting system evolution. The computations are terminated when the pressure somewhere in the system falls below lithostatic by a quantity exceeding the rock tensile strength [3,11,13] leading to rock collapse and either eruption closure or transition to highly transient dynamics. When that happens at chamber level, caldera collapse is expected [11,13].
The right panels in Figure 1 show one example of calculations (only a few times are reported). At each time step the flow variable distributions display well-known behaviours (e.g., [18]) reflecting a large gradient frictiondominated region leading to magma fragmentation, followed by an inertiadominated region with lower gradients close to conduit exit. Lumped pressure in the shallow chamber corresponds to both dyke exit and conduit base pressure, whereas magma mixing in the chamber results in a jump of gas volume fraction and (not shown) magma density, and by continuity, magma ascent velocity. Figure 2, obtained using the Dakota software [1] through a latin hypercube sampling within the range of employed system conditions, reports the two computed major quantities (defined in the Methods section) magnitude (a measure of the mass of discharged magma) and intensity (a measure of mass flow-rate ). For eruptions with intensity < 10.5, a reported dominance of unsteady dynamics (e.g., [12,16,20,21]) likely explains their absence in the simulated set. Apart from such low intensity eruptions, the observed magnitude-intensity relationship is well reproduced from relatively small, frequent M4 to globally impacting, destructive, rare M8 eruptions, suggesting that our simulations capture, on a first order, the major factors controlling the evolution of explosive volcanic eruptions. The colored symbols refer to cases for which the termination of the simulation is due to deep (red) or shallow (blue) chamber collapse (termination at volcanic conduit level is rare in our computations, and none for the deep dyke). Although with ample superposition, and in line with observations, deep chamber collapse is mostly associated to large magnitude eruptions, and characterises about all of the globally impacting M7+ eruptions.

Control by deep interconnections
The extent of magma chamber interconnections determines the transfer capability of magma across the reservoirs, exemplified in our simulations by the dyke width parameter. Dyke width is therefore a proxy for the effective mass flow rate across a complex system of interconnections possibly involving multiple dykes and intermediate storage zones [8]. Accordingly, parameterization of dyke width means parameterizing the overall interconnection efficiency. With small dyke width magma transfer across reservoirs is poorly efficient, thus the eruption discharge rate is not compensated by deep magma arrival (Fig. 3b,d, cold colours). As a consequence, shallow chamber pressure decays rapidly (Fig. 3a), while deep chamber pressure changes only slightly ( Fig. 3c) 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.
With larger dyke width (10-12 m in Fig. 3), the increased efficiency of deep mass transfer sustains the eruption discharge rate (Fig. 3b, d, light blue/cyan curves), leading to longer sub-steady eruption phase (up to about 8 hours in the figure with nearly constant mass flow rate). Deep chamber pressure decrease becomes significant (Fig. 3c) due to increased magma withdrawal, eventually reaching the rock collapse threshold (starting from the cyan curve, or dyke width of 12 m).
With further increase in dyke width, the mass flow rate along the dyke can overcome the one along the shallow conduit (hot-coloured lines in Fig.  3b,d), causing deep chamber depressurization (Fig. 3c) and shallow chamber pressurization (Fig. 3a). In turn, the latter causes the eruption discharge rate to increase (Fig. 3b). However, shallow pressure increase and deep pressure decrease rapidly lead to decreased of mass flow rate along the dyke (Fig. 3d), 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 towards shallow pressure decrease (Fig. 3a) and about constant or slightly decreasing mass flow rate (Fig. 3b).

Shallow versus deep controlled eruption regimes
The dynamics illustrated above show that besides explaining the observed magnitude-intensity relationship, 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, obscillate, remain constant, or increase depending on such an interplay. Non-linear relationships dominate the dynamics: with sufficiently small dyke width, its increase produces longer eruption duration; conversely, with larger dyke width, its further increase results in shorter eruption duration. The longest durations correspond to near-balance of conduit and dyke mass flow rates ( Fig.  3b,d). Eruption magnitude (reported in panel 3b) is also non-linearly related to dyke width: for small dyke width, its increase results in larger eruption magnitude, whereas for dyke widths sufficiently large to cause eruption termination by deep chamber collapse, the eruption magnitude becomes poorly or not dependent on further dyke width changes.
This situation depicts the existence of two regimes for explosive volcanic eruptions (Fig. 4): a shallow-dominated regime determined by low interchamber magma transfer efficiency and evolving to shallow chamber collapse, characterised 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; and a deep-dominated regime determined by large inter-chamber magma transfer efficiency and evolving to deep chamber collapse, characterised by initial shallow pressure and eruption mass flow rate increase, significant pressure decrease in the deep reservoir, negative dependence of eruption duration and poor or no dependence of eruption magnitude on magma transfer efficiency.

Implications for volcanic hazard
We have shown that the efficiency of magma transfer between reservoirs, controlled by deep interconnections, plays a central role in determining the eruption magnitude and evolution, and ultimately the impacts on the surrounding environment. Unfortunately, we do not have yet any means to image with sufficient accuracy the complexity of volcanic plumbing systems.In fact, we know little or nothing about the deep system of interconnections extending down towards deeper magmatic reservoirs. Our results show that without a reliable and accurate estimate of the corresponding magma transfer efficiency, we may never be able to forecast the eruption magnitude, intensity and duration, or its evolution towards large caldera collapse, no matter the type, intensity and evolution of the observed signals. As a matter of fact, no robust, generally accepted method to forecast the size of an impending eruption has been produced to-date. Previous analysis [19] suggests that the highly non-linear nature of volcanic processes limits deterministic predictions of volcanic eruption size. This study shows that even an the deterministic approach employed here leads to conclude that finding confident relationships between unrest dynamics and size of the impending eruption may continue to be a very hard task, if not a hopeless undertaking.   . Among the simulated points, blue refers to cases for which eruption closure is due to shallow chamber collapse, and red to cases for which eruption closure is due to deep chamber collapse.    [17] A. Neri, P. Papale, and G. Macedonio

Methods
Our computations refer to the sub-steady eruption phases for which the time scale of system evolution is significantly longer than the transit time of magma in the simulated domain. 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. The fundamental equations describing transport of mass and momentum along a conduit or dyke under 1D, steady, isothermal, multiphase conditions are reported in [8]. A 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. The computation is terminated when the magmatic pressure somewhere in the simulated domain falls below the local lithostatic pressure by a quantity exceeding the rock tensile strength, here fixed to an average value of 20 MPa [2,4,5]. 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 characterising 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 characterised by highly transient dynamics during which the volcanic conduit and eruption plume are established. These transient dynamics are not considered in our modelling. 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 is 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 modelling) to mass flow rate, and from that, to chamber mass loss or gain (see below). In particular, after update of mixture density given by where ρ m is mixture density, t is time, m is mass inside the chamber, v is chamber volume, the corresponding pressure is calculated. In fact, 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 [10]. 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 whereṁ is mass flow rate. Chamber volume is approximated by where β is the effective bulk modulus of the elastic walls surrounding the chambers. By repeating the same simulations with β = 10 10 Pa [2] and β = ∞ (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 S3). Therefore, we refer to only the latter case (β = ∞, non-deformable chamber walls) throughout this paper.
4. The time step ∆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 ∆t = 2 * t trans .
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: where tṁ is the characteristic time of variation of the mass flow rate over ∆t In spite of its relative simplicity, such a set up captures the fundamental characteristic of many volcanic plumbing systems characterised 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 (see Fig. 2). Magnitude (M) is computed as in [11] as M = log[m T (kg)] − 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, see panel b in Fig. 3. Intensity (I) refers to the eruption mass flow rate, and is defined as in [11] as I = log[ṁ(kg/s)] + 3, whereṁ 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 collapse is expected to occur, terminating the sub-steady phase of the eruption. When such a condition is found in correspondence of one of the magma chambers, a caldera collapse may occur. In the real world the processes may be more complex: in some cases the 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 new instabilities. Similarly, if large pressure decrease occurs locally in the volcanic conduit (typically close to magma fragmentation where the pressure gradient is the highest [9]), that may lead to local conduit wall erosion and conduit shape changes [1,6] and to transient eruption dynamics. Additional complexities relate to the mechanics of rock failure [5], especially in layered environments [4], whereby the generation of ring faults and the occurrence of caldera collapse depend on a number of other factors including chamber shape. Because we only refer to steady flow dynamics and do not solve the rock mechanics associated to 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 that employed by previous authors [2,3,7] who considered magma withdrawal from one single chamber.
The CONDUIT4 code [8] used for the simulations of magma flow along the conduit and dyke is available upon request to the authors. The SOLW-CAD code [10] for multi-component volatile saturation is embedded in the CONDUIT4 code, and openly available from www.pi.ingv.it/progetti/ eurovolc.  Figure S1: Time series for density and gas volume. Calculated evolution of magma density (a and c, referred to the gas-melt mixture) and gas volume (b and d) inside the shallow and deep chamber, for the same simulation cases as in Fig. 3.   Figure S3: The effects of chamber wall elasticity on computed eruption magnitude and duration. Open symbols as in Fig. 4 correspond to assuming rigid walls. Solid symbols assume instead elastic walls with bulk modulus equal to 10 10 Pa. Note that most of the solid symbols are superimposed to the open ones.