Evolution of the magma feeding system during a Plinian eruption: The case of Pomici di Avellino eruption of Somma–Vesuvius, Italy

Earth and Planetary Science Letters 482 (2018) 545–555

Contents lists available at ScienceDirect

Earth and Planetary Science Letters www.elsevier.com/locate/epsl

Evolution of the magma feeding system during a Plinian eruption: The case of Pomici di Avellino eruption of Somma–Vesuvius, Italy S. Massaro a , A. Costa b , R. Sulpizio a,∗ a b

Dipartimento di Scienze della Terra e Geoambientali, Via Orabona 4, 70125, Bari, Italy Istituto Nazionale di Geofisica e Vulcanologia, Via Donato Creti 12, 40128, Bologna, Italy

a r t i c l e

i n f o

Article history: Received 18 March 2017 Received in revised form 11 November 2017 Accepted 14 November 2017 Available online xxxx Editor: T.A. Mather Keywords: conduit/feeding system sustained explosive eruption eruption behaviour driving pressure partition lithic erosion

*

a b s t r a c t The current paradigm for volcanic eruptions is that magma erupts from a deep magma reservoir through a volcanic conduit, typically modelled with fixed rigid geometries such as cylinders. This simplistic view of a volcanic eruption does not account for the complex dynamics that usually characterise a large explosive event. Numerical simulations of magma flow in a conduit combined with volcanological and geological data, allow for the first description of a physics-based model of the feeding system evolution during a sustained phase of an explosive eruption. The method was applied to the Plinian phase of the Pomici di Avellino eruption (PdA, 3945 ±10 cal yr BP) from Somma–Vesuvius (Italy). Information available from volcanology, petrology, and lithology studies was used as input data and as constraints for the model. In particular, Mass Discharge Rates (MDRs) assessed from volcanological methods were used as target values for numerical simulations. The model solutions, which are non-unique, were constrained using geological and volcanological data, such as volume estimates and types of lithic components in the fall deposits. Three stable geometric configurations of the feeding system (described assuming elliptical cross-section of variable dimensions) were assessed for the Eruptive Units 2 and 3 (EU2, EU3), which form the magmatic Plinian phase of PdA eruption. They describe the conduit system geometry at time of deposition of EU2 base, EU2 top, and EU3. A 7-km deep dyke (length 2a = 200–400 m, width 2b = 10–12 m), connecting the magma chamber to the surface, characterised the feeding system at the onset of the Plinian phase (EU2 base). The feeding system rapidly evolved into hybrid geometric configuration, with a deeper dyke (length 2a = 600–800 m, width 2b = 50 m) and a shallower cylindrical conduit (diameter D = 50 m, dyke-to-cylinder transition depth ∼2100 m), during the eruption of the EU2 top. The deeper dyke reached the dimensions of 2a = 2000 m and 2b = 60 m at EU3 peak MDR, when the shallower cylinder had enlarged to a diameter of 60 m and a transition depth of 3000 m. The changes in feeding system geometry indicate a partitioning of the driving pressure of the eruption, which affected both magma movement to the surface and dyke growth. This implies that a significant portion of the magma injected from the magma chamber filled the enlarging dyke before it erupted to the surface. In this model, the lower dyke acted as a sort of magma “capacitor” in which the magma was stored briefly before accelerating to the cylindrical conduit and erupting. The capacitor effect of the lower dyke implies longer times of transit for the erupting magma, which also underwent several steps of decompression. On the other hand, the decompression of magma within the capacitor provided the driving pressure to maintain the flow into the upper cylindrical conduit, even as the base of the dyke started to close due to the drop in driving pressure from progressive emptying of the magma chamber. The shallower cylindrical conduit was shaped through the erosion of conduit wall rocks at and above the fragmentation level. Using the lithic volume and duration of EU3, the erosion rate of shallower cylindrical conduit was calculated at ∼5 × 103 m3 /s. The outcomes of this work represent an important baseline for further petrologic and geophysical studies devoted to the comprehension of processes driving volcanic eruptions. © 2017 Elsevier B.V. All rights reserved.

Corresponding author. E-mail address: [email protected] (R. Sulpizio). https://doi.org/10.1016/j.epsl.2017.11.030 0012-821X/© 2017 Elsevier B.V. All rights reserved.

546

S. Massaro et al. / Earth and Planetary Science Letters 482 (2018) 545–555

1. Introduction During volcanic eruptions, magma erupts from a chamber through a volcanic conduit, which is the path through the Earth’s crust that connects the chamber to the surface. The geometric configuration of the coupled magma chamber–conduit system, which represents the feeding system of an eruption, is controlled by the interaction of magma with host rocks. Contrarily to the common assumption of a rigid unaltered cylindrical conduit, the feeding system can be modified during the evolution of the eruption, affecting the dynamics of magma ascent (e.g. Macedonio et al., 1994; Costa et al., 2007a, 2009). The time evolution of conduit/feeding systems during a volcanic eruption are still far from being understood. In particular, understanding how the ascending magma interacts mechanically and thermally with the surrounding rocks (e.g., Campbell et al., 2013) is of paramount importance for our knowledge of the physical processes that drive volcanic eruptions. An understanding of this interaction will shed light on some pivotal questions in modern volcanology research, such as: i) how do volcanic conduits evolve, and, ii) how do the interactions of magma with host rocks control eruption dynamics? The effect of conduit walls erosion on eruption dynamics was investigated in some pioneering works (Varekamp, 1993; Macedonio et al., 1994), while Campbell et al. (2013) investigated the thermal milling effect of pyroclastic mixtures on conduit walls. Costa et al. (2007a, 2007b, 2007c, 2009, 2011) described the effects of elastic variations of the conduit cross section on eruption dynamics, while de’ Michieli Vitturi et al. (2008) investigated the effects of changing radius for cylindrical conduits. The changes in fragmentation level during an explosive eruption were deduced from the study of lithic fragments in pyroclastic deposits, which represent one of the main sources of information about the nature of the crust underlying a volcano (Barberi et al., 1989). Because lithic fragments are assumed to be formed mainly at or above the fragmentation level (Barberi et al., 1989; Macedonio et al., 1994), changes in lithic type in pyroclastic successions have been interpreted as an evidence for changes in the depth of the fragmentation front during an explosive eruption (Barberi et al., 1989; Varekamp, 1993; Cioni et al., 2000; Sulpizio et al., 2010). Discrete lithic-enriched horizons can also reflect intermittent episodes of magma–water interaction due to increasing explosivity or flashing of liquid water to vapour in fluidhosting country rocks (e.g. Macedonio et al., 1994; Jurado-Chichay and Walker, 2001), or onset of caldera-forming phases (e.g. Cioni et al., 1999; Brown and Gardner, 2004). To date, erosion rates of country rocks during explosive eruptions are not available, and eroded lithics (preserved in pyroclastic deposits) have not been used to constrain numerical models of physical processes driving and governing explosive eruptions (physics-based models). Indeed, physics-based models have only been recently used in volcanology, and specifically for the inversion of large amounts of geophysical data relating to lava dome eruptions through the adoption of physics-based models with the sophisticated Bayesian approach (Anderson and Segall, 2011, 2013). Inversion of geological and volcanological data characterizing explosive eruptions is more challenging, due to the paucity of available data sets and the large uncertainties in the stratigraphic data. For these reasons, here we apply a physics-based model of volcanic conduit named CPIUD (Costa et al., 2009, 2011; Costa and Marti, 2016) with a less sophisticated approach as a first step in this direction, which consists of a grid search over a few geometry configurations describing the feeding system. The uncertainty of the inversion approach was reduced through all available geological and volcanological data either as constraints and target outputs for numerical simulations. This has allowed us to create

a comprehensive physics-based model linking magma properties, conduit evolution and eruption dynamics. This approach was applied to the evolution of the conduit system during the Plinian phase of the Pomici di Avellino (PdA) eruption (Somma–Vesuvius, Italy), using geological (rock type; Sulpizio et al., 2010), and volcanological (mass discharge rate (MDR), lithic volume, median terminal velocity of particles in fall deposits; Sulpizio et al., 2010; Hanson, 2016) data to constrain outputs from physical modelling (CPIUD code; Costa et al., 2009). The PdA Plinian phase (Cioni et al., 2000; Sulpizio et al., 2010; Sevink et al., 2011) represents an excellent opportunity to reconstruct the evolution of the conduit/feeding system during a large explosive event. This is because it has the following advantages: i) it was a simple eruptive phase, with establishment of a sustained column that drove continuous deposition of pyroclastic fall deposits (only one, very small pyroclastic density current (PDC) is interbedded; Sulpizio et al., 2010); ii) it was driven by purely magmatic fragmentation; iii) detailed studies described the stratigraphy and the physical volcanology of the eruption (Cioni et al., 2000; Sulpizio et al., 2010); iv) the subsurface stratigraphy of the country rocks is well known (Fig. 1; Santacroce, 1987); and v) the pre-eruptive conditions of the magma are well constrained (Santacroce et al., 2008; Scaillet et al., 2008; Pappalardo and Mastrolorenzo, 2010). The CPIUD model can simulate steady flows of magma inside elastic conduits that have elliptical cross-sections with major and minor semi-axes that change with depth. This also allows the description of a dyke-shaped conduit (a subvertical magma-filled crack) by approximating it as an elongated ellipse. Quasi-static elastic deformation of the dyke is accounted for by an analytical solution that couples cross-sectional area with the magma pressure. The model allows input into the conduit bottom of either liquid magma or pyroclastic particulate into the conduit bottom, and solves the mass and momentum conservation equations. Input parameters include the depth of magma chamber top (i.e. conduit length), overpressure at the top of the magma chamber, rheological properties of magma (e.g. melt and magma viscosity), and geometry of the conduit/feeding system (pure dyke or deeper dyke passing into a shallower cylinder). The MDRs calculated with volcanological methods (Sulpizio et al., 2010) were used as target parameters for the numerical simulations. The amount of lithic ejecta and their provenance depths (indicative of fragmentation level) were used to constrain the possible ranges of model solutions as functions of different geometries of the conduit/feeding system. The model outputs provided snapshots of the geometry of the conduit/feeding system during the Plinian phase of the PdA eruption, highlighting the passage from a pure dyke to a more hybrid geometry composed of a deeper dyke that evolves into a shallower cylinder. The different geometric configurations are also discussed in the light of dyke propagation due to magma injection from the chamber. Combining the results from different snapshots, the emerging model depicts, for the first time, the evolution of a feeding/conduit system during a sustained Plinian eruption. In particular, the relations between conduit/feeding system evolution and eruption dynamics highlight how the lower dyke acts as a sort of magma capacitor, which can store the magma for a while before it is conveyed to the fragmentation level around the base of the shallow cylinder conduit. The latter is shaped through time by lithic erosion, and modulates the MDR. The capacitor effect of the dyke could sustain the eruption even after the connection to the magma chamber is interrupted through the closure of the dyke base by fall of the overpressure below the lithostatic load. The paper is organised as a summary of the PdA eruption, a description of input data and constraints, a report of results, discussion and conclusive remarks. Three Appendixes provide more de-

S. Massaro et al. / Earth and Planetary Science Letters 482 (2018) 545–555

547

tails on PdA stratigraphy (Appendix A), CPIUD code (Appendix B), and model sensitivity to different input parameters (Appendix C). 2. The Avellino eruption The Pomici di Avellino eruption (PdA) is the third Plinian eruption produced by Somma–Vesuvius, dated at 3945 ± 10 cal yr BP (Sevink et al., 2011). First recognised by Johnston Lavis (1884), it occurred during the Early Bronze Age, as indicated in numerous associated archaeological sites (Cioni et al., 2000; Di Vito et al. 2009). Sulpizio et al. (2010) described the stratigraphic succession formed by three main eruptive phases (opening, magmatic Plinian, and phreatomagmatic), which contain five Eruptive Units (sensu Fisher and Schminke, 1984; EU1–EU5; Appendix A). Short-lived sustained columns occurred twice during the opening phase (EU1a, b; Ht of 13 and 21.5 km, respectively; Sulpizio et al., 2010), and dispersed thin fall deposits together with small PDCs onto the volcano slopes. A thin ash bed separates the EU1 and the beginning of the Plinian phase, representing a pause in eruptive activity that allowed finer material to settle out. The magmatic Plinian phase produced the main volume of the erupted deposits, emplacing white phonolitic (EU2) and grey tephri-phonolitic (EU3) pumice fallouts. During EU2, the eruption intensity increased with time, resulting in a reversely graded deposit (Sulpizio et al., 2010). The estimated maximum column height is of ∼23 km, with a MDR of ∼5.7 × 107 kg/s (Sulpizio et al., 2010). The fall deposits have a total estimated volume of 0.3 km3 (Sulpizio et al., 2010). Hanson (2016) calculated the amount of erupted lithic clasts during EU2 to be ∼0.003 km3 . The passage from EU2 to EU3 is indicated in the field by a thin layer of white-grey banded pumice clasts, indicative of syn-eruptive mingling between the EU2 and EU3 magma sources. Column height reached ∼31 km in EU3, corresponding to a MDR of ∼1.7 × 108 kg/s (Sulpizio et al., 2010). The EU3 fall deposits have an estimated total volume of ∼1.0 km3 . Total volume of erupted lithic clasts during EU3 was assessed at ∼0.03 km3 (Hanson, 2016). Only one small PDC deposit, due to partial column collapse was recognised in the EU3 succession (Sulpizio et al., 2010). A short eruptive pause occurred after the end of the EU3, which allowed a thin ash bed to be deposited. The eruption resumed with the establishment of a short-lived sustained column that emplaced the EU4 fall deposits (Sulpizio et al., 2010). The PdA eruption ended with a phreatomagmatic phase (EU5, Sulpizio et al., 2010), characterised by extensive generation of PDCs, with total erupted volume estimated at ∼1 km3 (Gurioli et al., 2010). Experimental petrology data indicate a magma storage pressure of 200 MPa (Scaillet et al., 2008; Balcone-Boissard et al., 2012). This pressure value corresponds to a chamber top located at about 7 km, similarly to those estimated for other Plinian eruptions of Somma–Vesuvius (Scaillet et al., 2008). 3. Input data and constraints A detailed summary of CPIUD code (Costa et al., 2009) is provided in Appendix B, while a sensitivity analysis of the code outputs to variable input parameters is contained in Appendix C. The input parameters for the CPIUD code were compiled from available geophysical, petrological, and stratigraphic datasets for both EU2 and EU3, and are summarised in Table 1. The total length of the conduit system (conduit depth) was set at 7 km based on experimental petrology data (Scaillet et al., 2008; Table 2; Fig. 1). The magma chamber was assumed to be an ellipse of 2 km (height) × 4 km (width), and confined within limestones (Fig. 1), as known from the lithological subsurface data (Santacroce, 1987). According to experimental petrology (Scaillet

Fig. 1. Sketch of the subsurface stratigraphy of the Somma–Vesuvius area as inferred from Trecase borehole log data (from Santacroce, 1987). The schematic magma chamber is not to scale and its top is placed at 7 km below Somma–Vesuvius summit (Scaillet et al., 2008). The schematic conduit shows the position of fragmentation level during EU2 base, EU intermediate, EU top, and EU3 inferred from lithic type (Sulpizio et al., 2010). The black arrows indicate the shallowing of fragmentation level in the early stage of the Plinian phase and its deepening starting at the middle of EU2.

et al., 2008) the initial pressure at the top of the magma chamber was set at 200 (±20) MPa, lowering to 175 MPa in the final phase of each simulation. Magma temperatures for the phonolitic (EU2) and tephra-phonolitic (EU3) magma were assumed at 820 ◦ C (1093 K) and 850 ◦ C (1123 K), respectively (Pappalardo and Mastrolorenzo, 2010; Table 1). The viscosity of the melt was calculated with the Giordano et al. (2008) model, while the effective magma viscosity (melt + crystals) was calculated with a crystal mass fraction in the chamber less than 20% (Sulpizio et al., 2010; Pappalardo and Mastrolorenzo, 2010). The assessed effective magma viscosity is in the order of 107 Pa s (EU2) and 106 Pa s (EU3). We estimated the host rock rigidity of limestone using the empirical relationship between static and dynamic ranges of the Young Modulus (Wang, 2000), and a Poisson’s ratio ν equal to 0.3 (Table 1). We need to highlight that the solutions of the inverse problem are non-unique (e.g. Anderson and Segall, 2011, 2013). Possible solutions were selected using target values of MDR and constraints from volcanological and stratigraphic data. The target peak MDRs are 5.7 × 107 kg/s for the EU2 top and 1.7 × 108 kg/s for the EU3 (Sulpizio et al., 2010, values obtained using the method of Sparks, 1986). No estimate of the MDR was available for the initial phase of the EU2, but it is assumed to be the same as that of the opening phase of PdA (EU1a–b, average MDR of ∼1 × 107 kg/s; Sulpizio et al., 2010). This assumption is justified because of the similar grain size of EU1 and the base of EU2 (Sulpizio et al., 2010; Table 2). An uncertainty in MDR estimates of about ±50% was considered (Bonadonna et al., 2015; Costa et al., 2016; Table 2; Figs. 3–5), and our preferred matches

548

S. Massaro et al. / Earth and Planetary Science Letters 482 (2018) 545–555

Table 1 Input parameters used in the simulations. The data variability is reported in brackets. xtot , and T from Pappalardo and Mastrolorenzo (2010); xc , L, ρl0 , and ρc0 from Sulpizio et al. (2010); ρr from Santacroce (1987); μ obtained using Giordano et al. (2005); P ch from Scaillet et al. (2008); s, n, and β from Costa et al. (2009); G and ν from Touloukian et al. (1989). Notation

Description

xtot xc L

Values

Concentration of dissolved gas (wt.%) Crystal fraction (wt.%) Conduit length (m) Density of the melt phase (kg m−3 ) Density of the crystal (kg m−3 ) Host rock density (kg m−3 ) Magma temperature (K) Effective magma viscosity (Pa s) Initial magma chamber pressure (MPa) Final magma chamber pressure (MPa) Solubility coefficient (Pa−1/2 ) Solubility exponent Static host rock rigidity (GPa) Poisson ratio Bulk modulus of melt/crystal (GPa)

ρl0 ρc0 ρr T

μ Pi Pf s n G

ν β

EU2

EU3

2 (1.9 ± 0.2) 0.14 (0.11–0.17) 7000 2410 ± 30 2600 2600 1093 (1023–1259) 107 200 ± 20 175 4.1 × 10−6 0.5 13 0.3 10

1.4 (1.3 ± 0.1) 0.13 (0.11–0.15) 7000 2700 ± 30 2700 2600 1123 (1023–1259) 106 200 ± 20 175 4.1 × 10−6 0.5 13 0.3 10

Table 2 Matching values used for defining acceptable solutions from the CPIUD code. Bulk lithic volumes from Hanson (2016), obtained using thickness from B-spline contouring method (Engwell et al., 2015) and exponential decay (Pyle, 1989), Mass Discharge Rates from Sulpizio et al. (2010) using the method of Sparks (1986). Eruptive units

Bulk lithic volume (km3 )

Height (km)

MDR (kg/s)

EU3 EU2top EU2beg

0.02 (±30%) 0.003 (±30%) 0

31 23 21.5

1.7 × 108 (0.9–2.6 × 108 ) 5.7 ×107 (2.9–8.6 × 108 ) ∼107 (0.5–1.5 × 108 )

Table 3 Outputs of CPIUD simulations for EU2 and EU3 phases. MDR = mass discharge rate; FL = depth of fragmentation level; 2a = dyke length; 2b = dyke width; D = upper cylinder diameter. EUs

MDR (kg/s × 107 )

FL (m)

2a (m)

2b (m)

D (m)

EU3 EU2top EU2beg

17 (25.5–8.5) 5.7 (8.55–2.85) 1 (1.5–0.5)

3000 2100 5100

2000 600–800 200–400

120 50 10–12

120 50 0

(i.e. the solutions of the adopted physics-based model that best match the volcanological data) of MDRs obtained from the numerical simulations were chosen within these typical uncertainty ranges (Bonadonna et al., 2015; Costa et al., 2016). The constraints used for the CPIUD simulations derive from stratigraphic and volcanological data are summarised in Table 2 (Sulpizio et al., 2010; Hanson, 2016). The total volume of lithic fragments erupted during the Plinian phase of the PdA eruption was recently assessed at 0.003 ± 0.001 km3 and 0.03 ± 0.01 km3 for EU2 and EU3, respectively (Table 2; Hanson, 2016). Detailed variation of the relative abundance of different rock types (mainly lavas and carbonates) versus stratigraphic height was previously reported in Sulpizio et al. (2010), and it is summarised in Fig. 2. They are used as a proxy for depth of the fragmentation level, under the assumption that they form at or above the fragmentation level itself (Macedonio et al., 1994). The geometric constraints for the feeding system include the assumption that the Plinian phase started due to the arrival at the surface of a dyke connected with the magma chamber (Sulpizio et al., 2010). This was followed by the progressive development of a hybrid dyke/cylinder feeding system, whose geometry was shaped by erosion of wall rocks due to shock waves generation at the fragmentation depth (e.g., Macedonio et al., 1994) and the abrasive effect of the py-

Fig. 2. Variation of median value of terminal velocity (Mdφ TV ), abundance of lithic clasts and lava/carbonate ratio in the lithic clasts vs. normalised stratigraphic height in EU2 and EU3 deposits. More negative values of Mdφ TV indicate higher terminal velocities. Red data are from a composite proximal (∼6 km from the vent) stratigraphic section, while the blue data are from the medial (∼15 km from the vent) composite stratigraphic section (from Sulpizio et al., 2010, modified). The peak MDR value is reported for both EU and EU3. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

roclastic mixture above fragmentation depth (e.g., Campbell et al., 2013). For the sake of simplicity, all simulations of the hybrid feeding system were carried out assuming a transition zone from dyke to cylinder at different depths (from 2000 m up to 4500 m). The diameter of the cylinder region D was set equal to the dyke width 2b. 4. Results As previously mentioned, the proposed solutions of conduit/ feeding system geometry are non-unique, selecting the best compromise between the evidence obtained from the stratigraphic data (Sulpizio et al., 2010; Hanson, 2016) and the physical parameters used in the governing equations for flow in the conduit system (Costa et al., 2009). An iterative approach was used for achieving the best solutions for dyke/cylinder widths and lengths, chosen by

S. Massaro et al. / Earth and Planetary Science Letters 482 (2018) 545–555

549

Fig. 3. Diagram showing the calculated MDR at varying dyke cross-sectional axes (2a and 2b) for the pure dyke geometry assumed at EU2 beginning. The b/a ratio is a proxy for eccentricity of dyke cross-sectional area. A value greater than 1/20 is considered unrealistic based on observation of dyke length and width in nature (Rivalta et al., 2015). Solid symbols indicate MDR values respecting the b/a ratio constraint, while the open symbols fall below the acceptable b/a value and are not further considered. The dashed lines indicate the assumed range of MDR, calculated from simulated MDR ±50%. The green dashed lines indicate the minimum acceptable value of 2a and the corresponding MDR, taking into account the b/a constraint. The blue value of 2a indicates the minimum value of acceptable 2a. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

matching solutions against the target range of MDR and the lithic provenance depths (Table 3). 4.1. Base of EU2 Fig. 3 illustrates the results of the simulations carried out using a pure dyke geometry. The runs were carried out using the dyke length (2a) and dyke width (2b) as independent variables, and the MDR as output. The initial dyke width (2b) range was set from 5 m to 50 m (considering 5 m, 10 m, 20 m and 50 m). Runs with a circular plan shape (2a = 2b) to the conduit were not taken into account. It is worth noting that runs with 2b = 5 m have no stable solutions (Costa et al., 2009) for 2a > 200 m, while runs with 2b = 10 m have no solutions for 2a > 500 m. The target range of MDR (∼1 × 107 kg/s ± 50%; Table 2) only intersects curves corresponding to 2b equal to 20 m and 10 m (Fig. 3). Taking into account the common length-thickness ratio of dykes (in the range 60 and 250 for non-feeder dykes and even higher for feeder dykes, Kusumoto et al., 2013) all possible solutions were found excluding very unlikely values of 2b/2a > 1/20 (grey area), which then limits the acceptable values to the light blue area observed in Fig. 3. The limits of such an area correspond to the 2b range between 10–12 m and 2a range between 200–400 m. This geometrical combination corresponds to a fragmentation depth at 5100 m, a value consistent with the deposit componentry and stratigraphic evidence (Fig. 1). 4.2. Top of EU2 The simulations for the top of EU2 were performed allowing for hybrid dyke/cylinder geometries. The dyke length and width were

Fig. 4. CPIUD simulations for the EU2 top using hybrid (dyke plus cylinder) geometry. a) Diagram showing calculated MDR at varying dyke length (2a) and width (2b; squares = 10 m; diamonds = 20 m; crosses = 50 m; triangles = 100 m). The dykecylinder transitions are fixed at 4500 m (black lines), 3500 m (blue lines), 2000 m (orange lines). The grey area corresponds to the acceptable MDR range considering an uncertainty of ±50% (Bonadonna et al., 2015). b) Diagram showing fragmentation level (FL) vs. dyke length (2a) corresponding to four different dyke-cylinder transitions (triangles = 4500 m; crosses = 3500 m; squares = 2000 m). The best match solutions, indicated in green, have dyke length 2a between 600–800 m and dyke width 2b of 50 m, with dyke-cylinder transition at 2000 m. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

the same used for the calculations at the beginning of EU2, while the diameter of the cylinder D was set equal to 2b. The dykecylinder transition depth was varied from 2000 m, to 3500 m and 4500 m. Fig. 4 shows how only results of the runs with 2b = 50 m intersect the MDR of EU2 top (Sulpizio et al., 2010; Table 2), which was expanded by ±50% to take into account the uncertainty and variability in the MDR estimate (Bonadonna et al., 2015; Costa et al., 2016). However, Fig. 4a does not discriminate between the different transition depths, as all the curves intersecting the acceptable values of MDR. In order to determine the best solution for dyke/cylinder dimensions and the transition depth, we can consider the constraints coming from the fragmentation level

550

S. Massaro et al. / Earth and Planetary Science Letters 482 (2018) 545–555

a narrow range between 4.2 and 4.5 × 107 kg/s and dyke length 2a between 600 and 800 m (green shaded areas in Fig. 4a). The best match occurs with a fragmentation depth at 2100 m, corresponding to a MDR of 4.2 × 107 kg/s. The erosion of a cylinder starts from a dyke width 2b = 10 m and length 2a ∼ 200–400 m and evolves towards 2b = 50 m and 2a between 600 and 800 m, producing 0.003 km3 of lithic fragments, in agreement with calculations based on field data (Hanson, 2016; Table 2). This geometric configuration also accounts for erosion of carbonates, and satisfies, at least semi-quantitatively, the observed lava/carbonate ratio in the lithic clasts (Fig. 2). 4.3. EU3

Fig. 5. CPIUD simulations for the EU3 using hybrid (dyke plus cylinder) geometry. a) Diagram showing MDR at varying dyke length (2a) and width (2b; squares = 60 m; diamonds = 80 m; crosses = 100 m; triangles = 120 m). The dyke-cylinder transitions are fixed at 4000 m (black line), 3000 m (blue line), 2000 m (orange line). The grey area corresponds to the acceptable MDR considering an uncertainty of ±50% (Bonadonna et al., 2015). b) Diagram showing fragmentation level (FL) vs. dyke length 2a. The green dashed line indicates the best match solution with dyke length 2a of 2000 m, dyke width 2b of 120 m, and a dyke-cylinder transition of 3000 m. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

(Fig. 4b). Volcanological data show that the fragmentation depth moves from shallower to deeper levels in the conduit between intermediate and top of EU2 (Fig. 1). Assuming that for a sustained phase of the eruption the cylindrical part of the conduit forms around the fragmentation level due to the erosion/disruption of dyke walls (Costa et al., 2009), we choose the solution that had a fragmentation level around that transition depth. The analysis of Fig. 4b shows that the only transition curve that intersects the depth of the fragmentation level is the 2000 m curve, whereas all the others yield fragmentation levels shallower than the transition itself. The combination of Figs. 4a and 4b determines the possible range of MDR, fragmentation depth, dyke length (2a), dyke width (2b) and cylinder diameter values for a transition dyke/cylinder at 2000 m. In particular, the MDR has

The CPIUD runs for the EU3 take into account that the volcanological data indicate a deepening of the fragmentation level and an enlargement of the conduit cross-section (increase in MDR) with respect to the preceding EU2 top. The dyke width (2b) was explored using values of 60 m, 80 m, 100 m, and 120 m. The diameter of the cylinder D was set equal to the underlying dyke width 2b. The dyke-cylinder transition depths were considered at 2000 m, 3000 m, and 4000 m. Using the same approach illustrated for EU2 top, Fig. 5a shows that only the runs with transitions at 3000 m and 4000 m are compatible with the MDR estimated by Sulpizio et al. (2010) (expanded by ±50%, Table 2). In particular, the MDR range intersects the dyke width 2b values of 100 m and 120 m for a level of transition from dyke to cylinder (TL) at 4000 m, and 2b = 120 m for the TL at 3000 m. The analysis of Fig. 5b indicates the 3000-m transition curve is the best match because it is the only one at which the fragmentation and transition level have the same depth. The combination of graphs in Figs. 5a and 5b shows the best match of MDR, fragmentation depth, dyke length 2a, dyke width 2b and cylinder diameter values are found with a transition dyke-cylinder at 3000 m. In particular, the MDR coincides with the published one at 1.7 × 108 kg/s, where dyke length 2a = 2000 m (dyke width, 2b, and cylinder diameter, D, equal to 120 m). The fragmentation depth is at 3000 m, coinciding with the transition depth. The enlargement of the cylinder from a diameter of 50 m to 120 m, and its deepening from 2000 m to 3000 m are in good agreement with a lithic clast production of 0.03 km3 , as calculated from field data (Hanson, 2016; Table 2). This geometric configuration also accounts for erosion of carbonates, and satisfies the observed lava/carbonate ratio of the lithic clasts (Fig. 2). 5. Discussion Our preferred matches (i.e. the CPIUD solutions that best reproduce the volcanological data) for EU2 and EU3 highlight how the changes in the geometry of the feeding system due to dyke propagation and enlargement, and widening of the upper cylinder due to wall-rock erosion, has a strong effect on MDR and fragmentation level depth. In particular, it is important to note how the geometry of the feeding system changes during the different phases of the eruption (Fig. 6). For the sake of clarity, the evolution of the deeper dyke and shallower cylindrical conduit will be discussed separately in the following sections. It is also important to bear in mind that each geometrical setting illustrated in Fig. 6 is a “snapshot” capturing an established phase of the eruption, given that the results are based on the assumption of steady-state conditions (Costa et al., 2009). 5.1. Deeper dyke evolution The mechanisms driving the deeper dyke evolution must be analysed in the context of the partitioning of the driving pressure,

S. Massaro et al. / Earth and Planetary Science Letters 482 (2018) 545–555

551

Fig. 6. Cartoon showing the three “snapshots” of the PdA feeding system. On the right side, evolution of Mass Discharge Rate (MDR, red line), volume of eroded lithic clasts (black line), fragmentation level (FL, purple line), depth of transition level (TL, orange line), dyke length (2a, blue line), and dyke width (2b, green line) are shown. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

here defined as the magma pressure available for doing mechanical work. The CPIUD outputs earlier described indicate the crosssectional area of the feeding dyke increased from 3.8 × 103 m2 at the onset of EU2 to 1.9 × 105 m2 during EU3. This growth of cross-sectional area reflects an increase of MDR, which passes from ∼107 kg/s to ∼108 kg/s from the beginning of EU2 to the climax of EU3. The enlargement of the cross-sectional area occurred steadily throughout EU2, as testified by the continuous increase of Mdφ TV (assumed to be a good proxy of MDR variation) during EU2 (Fig. 2). This is also in agreement with the mass flow conservation required to maintain the sustained column conditions during EU2 and EU3. It follows that the driving pressure ( P d ) not only main-

tained mass flow through the conduit/feeding system, but part of it was also used for increasing the dyke cross-sectional area. The pressure partition during the eruption may be expressed by modifying the equilibrium equation (Gudmundsson, 2012; Sulpizio and Massaro, 2017):

P d = P open + P f + P η + P ρ + σ3

(1)

where P d is the driving pressure, P open is the pressure required to keep the dyke open, P f is the pressure spent fracturing the rocks, P η the viscous pressure, P ρ = (ρr − ρm ) gh is the buoyancy term, and σ3 is the minimum component of the stress in the crust.

552

S. Massaro et al. / Earth and Planetary Science Letters 482 (2018) 545–555

Once the dyke is open to the surface and the mass flux established, the fracturing pressure acts only laterally to increase the dyke length 2a. This allows for grouping of P open and P f as the pressure component acting on the horizontal plane ( P a–b ), while P η and P ρ define the vertical component of the driving pressure ( P z ):

P d = P a–b + P z + σ3

(2)

The partition between horizontal and vertical driving pressure may change during the eruption, based on changes in the dyke geometry. Considering σ3 is invariant at the scale of the eruption, the analysis of the two other terms of the driving pressure yields some interesting clues about the eruption dynamics. The vertical component of driving pressure P z controls the vertical flow of magma in the dyke, provided that the following condition is satisfied:

Pρ Pη

>1

(3)

The analysis of the horizontal partition of the driving pressure ( P a–b ) provides additional insights into the cross-section evolution. P a–b is composed of P open and P f , where these two terms can be analysed separately. P open relates to b (minor semi-axis of the dyke) following the equation (e.g. Gonnermann and Taisne, 2015):

P open =

E

b

2(1 − ν 2 ) a

(4)

where ν is the Poisson’s number, E the Young’s modulus, and a and b are the major and minor semi-axes of the dyke. Equation (4) states that a and b vary accordingly, and that the pressure needed to keep the dyke open is proportional to the b/a ratio. The b/a ratio remains nearly constant in our preferred solutions for the PdA Plinian phase equating to 0.05 at beginning of EU2, 0.08 at EU2 top, and 0.06 in the EU3. This emphasizes that the pressure within the dyke experiences little variation, and qualitatively mirrors the variation in fragmentation level depth. The fracturing pressure ( P f ) can be written as (e.g. Gonnermann and Taisne, 2015):

Kc Pf ≈ √ a

(5)

where K c is the fracture toughness, or the minimum stress at which the fracture occurs or a crack propagates (value between 1 and 100 MPa; Rivalta et al., 2015; Gonnermann and Taisne, 2015), which depends on pressure gradient within the dyke, and on the scale of the fracture (e.g., Rivalta et al., 2015). Although K c is not constant, its variation is minimal with respect to a. The equation (5) states that P f is high when a is short (beginning of dyke expansion) and decreases following a power law with increasing a. It follows that the increase in cross-sectional area cannot be infinite, and the b/a ratio of equation (5) will approach some limit value that will stop the dyke propagation. This general behaviour is in agreement with that inferred during the EU2 eruption, which shows an increase of MDR and the contemporaneous enlargement of the dyke cross-section, until it approaches a nearly constant value in EU3. 5.2. Upper conduit evolution CPIUD simulations for the EU2 base were carried out using simplified dyke geometries. This is a simplification of reality, because geological evidence highlights that many near-surface volcanic conduits are approximately cylindrical and become elongated dykes with depth (Costa et al., 2007a, 2009, 2011; Melnik et al., 2008; Melnik and Costa, 2014; Rivalta et al., 2015). Therefore, it

is not possible to exclude the presence of a short cylindrical conduit on top of the main dyke feeding the onset of the PdA Plinian phase, but this does not significantly change the simulation results. Once the sustained explosive activity was established, erosion of the original dyke proceeded through mechanical breaking of country rocks near the fragmentation level or by erosion from magma and pyroclasts along the walls (Macedonio et al., 1994; Costa et al., 2009; Campbell et al., 2013). In the case of an explosive eruption such as the PdA, the continuum magma flow in the deeper conduit passes into a pyroclastic mixture above the fragmentation level. Irrespective if it is due to magma or pyroclastic mixture, shear traction in an elliptical conduit (dyke) is not uniform on all of its walls, but rather changes with the position, being maximum at the minor axis tips and minimum at the major axis terminations. Assuming the erosion rate of the wall of the conduit wall is proportional to shear traction the conduit cross section will become more circular with time (Dragoni and Santini, 2007). Wallrock erosion rates have been calculated to be on the order of 10−7 m/s for Kilauean lava flows (about 2.5 m/month; Hulme, 1982), which closely approximates the erosion rates expected from the fluid magma flow in the deeper conduit. These estimates may be increased if viscous heating occurs at the magma–rock interface (Costa and Macedonio, 2005; Costa et al., 2007a). Considering the typical time frame of a sustained explosive eruption (on the order of few hours), it is clear that magma flow erosion rates are too slow below the fragmentation level to account for the observed amount of lithic fragments found in the EU2–EU3 deposits. It is instead in agreement with the amount of xenoliths found within pumice clasts, which is about 1 wt.% of the total lithic content (Hanson, 2016). It follows that the great majority of lithic fragments observed in EU2 and EU3 deposits are derived from erosion of country rocks due to fragmentation and motion of the pyroclastic mixture. At the fragmentation level, the pressure difference between the pyroclastic mixture and the lithostatic load of the surrounding rocks (Costa et al., 2009), combined with the formation of shock waves due to rupture of magma and the subsequent release of pressurised gas, may induce mechanical breaking of wall rocks (Macedonio et al., 1994). At the same time, the motion of the pyroclastic mixture induces shear traction and erosion on shallower conduit walls (Campbell et al., 2013). Theoretically, it is possible to calculate the total mass eroded in a volcanic conduit during a given phase of an explosive volcanic eruption as long as the erosion rate is known (Macedonio et al., 1994). Nevertheless, the erosion rate is not easy to estimate theoretically, being dependent upon many poorly constrained parameters such as the mechanical properties of wall rocks, the concentration of solids in the pyroclastic mixture, the stress state (pressurisation or under-pressurisation with respect to lithostatic pressure) of wall rocks, and the rock types in the volcanic edifice. In the case of the PdA eruption, erosion rates for EU2–EU3 can be obtained empirically using the published volcanological data (Sulpizio et al., 2010; Hanson, 2016). To this aim, the EU3 is particularly useful for calculating the erosion rate, having a known erupted mass (∼1012 kg) and MDR (1.7 × 108 kg/s; Sulpizio et al., 2010). The ratio between the mass of the EU3 and its MDR yields an eruption duration on the order of ∼7 × 103 s (∼2 h). It thus follows that the approximate mean erosion rate for EU3 is ∼5 × 103 m3 /s by dividing the EU3 bulk lithic volume (∼3 × 107 m3 ; Hanson, 2016) by the EU3 duration. This is the first derivation of an erosion rate during a phase of a sustained explosive eruption, and can be of paramount importance for the assessment of conduit shape, obtained using both numerical simulations and volcanological methods.

S. Massaro et al. / Earth and Planetary Science Letters 482 (2018) 545–555

553

Fig. 7. Relationships among magma flux and feeding system geometries. a) Diagram showing the upper cylinder conduit length (C l ) vs. magma volume in the dyke (V md ). The data fits with an exponential law, having α = 2.9 and a β = 0.0018. b) Diagram showing the cylinder volume (C v ) vs. magma volume in the dyke (V md ). The data fit with a linear trend, having χ = 7.1 and δ = 22. c) Diagram showing the mass discharge rate (MDR) vs. magma volume in the dyke (V md ). The data fit with a linear trend, having ε = 15 and φ = 0.21. d) Diagram showing the variation of erupted magma volume (V e ) vs. magma volume in the dyke (V md ). Solid line illustrates the cumulative erupted volume, while the dashed line shows the volume erupted in each phase (EU2 and EU3). Both regression lines show a linear trend, having γ = 2.3, η = 0.6, ι = 6.5, and λ = 0.5.

5.3. Geometric evolution of the feeding system and eruption dynamics The onset of the PdA eruption was characterised by two small explosive episodes witnessed in the EU1a and b deposits (Sulpizio et al., 2010), which indicate the initial difficulty to establish continuous mass flow through the feeding dyke. This behaviour can be explained using the findings of Costa et al. (2009), which demonstrated that a magma supply system with a simple dyke geometry is unlikely to sustain an explosive eruption for a long time, especially for evolved magma compositions, which tend to produce deep fragmentation levels. This is mainly because the underpressure generated at the fragmentation level tends to close the dyke and fracture the conduit walls (Costa et al., 2009). The restriction of the dyke cross-section and the wall collapse may drive a periodic blockage of the mass flux, which results in an intermittent discharge, as inferred for the EU1 case. The transition from pure dyke to cylinder geometry passing from EU2 base to EU2 top introduces significant differences in terms of deformation, fragmentation depth and flow rate. In the case of a dyke that evolves to a cylinder, the fragmentation level becomes shallower, which corresponds to what is interpreted during EU2, where the fragmentation depth rose into the volcanic pile (Fig. 1). After the formation of the shallower cylindrical conduit

the fragmentation level started to deepen and the cylinder became/continued to enlarge. The concomitant enlargement of the deeper dyke and shallower cylinder accounts for the deepening of the fragmentation level, required to maintain mass flow through the conduit system. Figs. 7a–d show how evolution of the deeper dyke and shallower cylinder are intimately related. In particular, the dyke and cylinder volumes have a first-order correlation (linear regression, Fig. 7b), where the dyke volume grew at a rate 22 times greater than the cylinder during the EU2–EU3 phases. Fig. 7c shows that the MDR increases with the volume of magma stored in the dyke. This relationship suggests that the dyke acted as a “capacitor” for magma rising from the chamber. The capacitor behaviour of the dyke is also confirmed by diagram of Fig. 7d, which shows how the cumulative erupted material increased at a slower rate with respect to the volume of magma in the dyke (Fig. 7d), which requires some form of temporary storage in the magma supply system. Moreover, the erupted volume in EU2 and EU3 (dashed line) is equal to half of the magma stored in the dyke. This result has important general consequences for the thermodynamics of an erupting magma, which would then experience longer decompression paths that might influence groundmass crystallization and vesiculation.

554

S. Massaro et al. / Earth and Planetary Science Letters 482 (2018) 545–555

5.4. What ended the Plinian phase? To fully described the evolution of the feeding system during the Plinian phase of PdA eruption, it is important to discuss what processes drove the sustained mass flux to cease Assuming no magma influx into the chamber during the eruption (in accordance with petrological evidence), it is important to bear in mind that any injection of magma volume into the feeding dyke resulted in a pressure drop in the chamber. The pressure drop is provided by the equation (e.g. Rivalta et al., 2015):

P =

M

1

M (βc + βm )

(6)

where  M is the mass of withdrawn magma, M is the original mass of the magma chamber, βc is the magma chamber compressibility and βm is the magma compressibility controlled mainly by magma bubbles. Both parameters can vary by an order of magnitude, depending on chamber size and depth (Melnik and Sparks, 2005). With time, the pressure loss drove the P d (equation (1)) below the lithostatic pressure, with consequent closure of the feeding dyke, starting from its base. The closure of the dyke feeding the EU2–EU3 is also signalled by the appearance of deep-seated rocks (skarn, marbles, cumulates) in the lithic components (Fig. 2; Sulpizio et al., 2010), because of the mechanical breaking of rocks at the base of the dyke and top of the magma chamber. The deepseated rocks appear in the lithic components in the upper 3/5 of the normalised stratigraphic height of EU3, suggesting that the sustained explosive phase was fed for a considerable time (around 1.5 h) after the onset of the closure of the dyke base. This is not surprising as long as the condition of equation (3) is maintained ( P ρ / P η >1). The buoyancy term P ρ is negative if we consider the density of EU3 magma (dense rock), which is around 2700 kg/m3 (Sulpizio et al., 2010; Table 1) and greater than those of carbonates and lavas forming the Somma–Vesuvius subsurface. Nevertheless, the decompression undergone by EU3 as it filled the feeding dyke induced gas exsolution and vesiculation, which significantly reduced the density of the magma and provided the necessary buoyancy pressure to maintain the mass flow rate through the upper conduit. This inference is supported by the componentry data, which shows a prevalence of vesicular pumice fragments throughout the EU3 stratigraphy (Sulpizio et al., 2010). The progressive emptying of the dyke probably induced generalised conduit wall collapse, which blocked the flow for a while and terminated the Plinian phase (EU2–EU3) of the PdA eruption. The temporary re-pressurisation of the system led to the establishment of a transient mass flow that produced the thin fall deposits of the EU4 (Sulpizio et al., 2010). 5.5. Some perspectives The physics-based model used for the description of the of the Plinian phase of PdA eruption can also be applied to other eruptions, provided the assumption of steady state conditions is achieved and enough constraints are available. The latter point is critical, because the solutions of the numerical model are nonunique and need to be constrained by geological and volcanological datasets. Further, the reconstruction of other geometric conditions of conduit/feeding systems may provide a better picture of physical processes controlling eruption dynamics. This is of particular importance because additional geometric constraints will help not only to explain variations of physical parameters at the vent, but can offer clues for interpreting geophysical data collected during an explosive eruption. The physics-based model for the PdA Plinian phase offers also some interesting hints for explaining observed pyroclastic rock texture. For example, the capacitor effect of the lower dyke promotes

prolonged the transit times for the magma to the surface, inducing a stepwise decompression pathway. This would result in a longer time available for groundmass crystallization and for bubble nucleation and growth. Reinterpretation of textures observed in pyroclastic rocks, as well as new studies aimed at understanding crystal and bubble growth kinetics should be compared with the results from these and other numerical simulations. This combined approach would allow for a comprehensive description of the physical and chemical processes occurring during an explosive eruption. Additionally, further research is required to characterise the physical mechanics of different rock types, allowing for increased understanding of erosion rates and the possible contribution of magma–rock chemical reactions to the dynamics of eruption. 6. Conclusions A first-order numerical approach combined with volcanological and geological data allowed, for the first time, the reconstruction of the evolution of conduit/feeding system during a Plinian eruption. The constrained numerical outputs depicted three different geometrical configurations of the conduit/feeding system at different times of the Plinian phase of the PdA eruption. The onset of the Plinian phase was reconstructed through a pure dyke geometry, which evolved into hybrid geometry (EU2 top) composed of a deeper dyke having width 2b = 50 m, length 2a ∼ 600–800 m, and a shallower cylindrical conduit (D = 2b) from 2000 m depth upward. The geometry of the feeding system changed until it reached the peak MDR in EU3, with a dyke width 2b = 120 m, length 2a = 2000 m, and a cylindrical conduit (D = 2b) that started at 3000 m depth. The enlargement of the dyke from EU2 base to EU3 indicates a dynamic, inter-dependent eruptive system composed of the magma chamber, a deeper dyke, and a shallower cylindrical conduit. In particular, the volume of magma within the lower dyke shows first-order relations with MDR and erupted volume. This indicates that the deeper dyke could act as a sort of capacitor for erupting magma, whose MDR was regulated by the development and growth of the upper cylindrical conduit. The capacitor effect of the lower dyke also guaranteed the continuous mass flux necessary for the sustained eruption even when the base of the feeding dyke started to close due to the lowering of driving pressure in the magma chamber. These results depict new ideas on how the feeding system evolves during a sustained explosive eruption, which may represent the baseline for developing new models of volcanic eruptions and for the comprehension of their driving processes. Acknowledgements All authors contributed to the conception and design of the study, acquisition of data, and analysis and interpretation of data, SM performed the numerical simulations. We thank Pierfrancesco Dellino for the revision of the draft of the manuscript. Associate Editor Tamsin Mather is acknowledged for handling the manuscript and two anonymous reviewers are thanked for a thorough revision that significantly improved the quality of the manuscript. Jonathan Hanson kindly provided access to data of his PhD thesis. Michael Ort and Madison Myers are warmly thanked for editing the English language. Appendix A. Supplementary material Supplementary material related to this article can be found online at https://doi.org/10.1016/j.epsl.2017.11.030.

S. Massaro et al. / Earth and Planetary Science Letters 482 (2018) 545–555

References Anderson, K., Segall, P., 2011. Physics-based models of ground deformation and extrusion rate at effusively erupting volcanoes. J. Geophys. Res. 116 (B7). https:// doi.org/10.1029/2010JB007939. Anderson, K., Segall, P., 2013. Bayesian inversion of data from effusive volcanic eruptions using physics-based models: application to Mount St. Helens 2004–2008. J. Geophys. Res. 118 (B5), 2017–2037. Balcone-Boissard, H., Boudon, G., Ucciani, G., Villemant, B., Cioni, R., Civetta, L., Orsi, G., 2012. Magma degassing and eruption dynamics of the Avellino pumice Plinian eruption of Somma–Vesuvius (Italy). Comparison with Pompeii eruption. Earth Planet. Sci. Lett. 331 (332), 257–268. Barberi, F., Cioni, R., Santacroce, R., Sbrana, A., Vecci, R., 1989. Magmatic and phreatomagmatic phases in explosive eruptions of Vesuvius as deduced by grain-size and component analysis of the pyroclastic deposits. J. Volcanol. Geotherm. Res. 38, 287–307. Bonadonna, C., Biass, S., Costa, A., 2015. Physical characterization of explosive volcanic eruptions based on tephra deposits: propagation of uncertainties and sensitivity analysis. J. Volcanol. Geotherm. Res. 296, 80–100. Browne, B.L., Gardner, J.L., 2004. The nature and timing of caldera collapse as indicated by accidental lithic fragments from the AD V1000 eruption of Volcan Ceboruco, Mexico. J. Volcanol. Geotherm. Res. 130, 93–105. Campbell, M.E., Russell, J., Porrit, L.A., 2013. Thermomechanical milling of accessory lithics in volcanic conduits. Earth Planet. Sci. Lett. 377–378, 276–286. Cioni, R., Santacroce, R., Sbrana, A., 1999. Pyroclastic deposits as a guide for reconstructing the multi-stage evolution of the Somma–Vesuvius Caldera. Bull. Volcanol. 60, 207–222. Cioni, R., Levi, S., Sulpizio, R., 2000. Apulian Bronze Age pottery as a long-distance indicator of the Avellino Pumice eruption (Vesuvius, Italy). Geol. Soc. (Lond.) Spec. Publ. 171, 159–177. Costa, A., Macedonio, G., 2005. Viscous heating effects in fluids with temperaturedependent viscosity: triggering of secondary flows. J. Fluid Mech. 540, 21–38. Costa, A., Melnik, O., Vedeneeva, E., 2007a. Thermal effects during magma ascent in conduits. J. Geophys. 112 (B12). Costa, A., Melnik, O., Sparks, R.S.J., Voight, B., 2007b. The control of magma flow in dykes on cyclic lava dome extrusion. Geophys. Res. Lett. 34, L02303. Costa, A., Melnik, O., Sparks, R.S.J., 2007c. Controls of conduit geometry and wallrock elasticity on lava dome. Earth Planet. Sci. Lett. 260 (1–2), 137–151. Costa, A., Sparks, R.S.J., Macedonio, G., Melnik, O., 2009. Effects of wall–rock elasticity on magma flow in dykes during explosive eruptions. Earth Planet. Sci. Lett. 288, 455–462. https://doi.org/10.1016/j.epsl.2007.05.024. Costa, A., Gottsmann, J., Melnik, O., Sparks, R.S.J., 2011. A stress-controlled mechanism for the intensity of very large magnitude explosive eruptions. Earth Planet. Sci. Lett. 310, 161–166. Costa, A., Marti, J., 2016. Stress field control during large caldera-forming eruptions. Front. Earth Sci. 4, 92. Costa, A., Suzuki, Y.J., Cerminara, M., Devenish, B., Esposti Ongaro, T., Herzog, M., Van Eaton, A., Denby, L.C., Bursik, M., de Michieli Vitturi, M., Engwell, S., Neri, A., Barsotti, S., Folch, A., Macedonio, G., Girault, F., Carazzo, G., Tait, S., Kaminski, E., Mastin, L., Woodhouse, M., Phillips, J.C., Hogg, A.J., Degruyter, W.J., Bonadonna, C., 2016. Results of the eruption column model inter-comparison study. J. Volcanol. Geotherm. Res. 326, 2–25. https://doi.org/10.1016/j.jvolgeores. 2016.01.017. Di Vito, M.A., Zanella, E., Gurioli, L., Lanza, R., Sulpizio, R., Bishop, J., Tema, E., Bonzi, G., La Forgia, E., 2009. The Afragola settlement near Vesuvius, Italy: the destruction and abandonment of a Bronze Age village revealed by archaeology, volcanology and rock-magnetism. Earth Planet. Sci. Lett. 277, 408–421. de ’Michieli Vitturi, M., Clarke, A.B., Neri, A., Voight, B., 2008. Effects of conduit geometry on magma ascent dynamics in dome-forming eruptions. Earth Planet. Sci. Lett. 272, 567–578. Dragoni, M., Santini, S., 2007. Lava flow in tubes with elliptical cross sections. J. Volcanol. Geotherm. Res. 160, 239–248. Engwell, S.L., Aspinall, W.P., Sparks, R.S.J., 2015. An objective method for the production of isopach maps and implications for the estimation of tephra deposit volumes and their uncertainties. Bull. Volcanol. 77, 61. Fisher, R.V., Schmincke, H.U., 1984. Pyroclastic Rocks. Springer, Berlin. Gurioli, L., Sulpizio, R., Cioni, R., Sbrana, A., Santacroce, R., Luperini, W., Andronico, D., 2010. Pyroclastic flow hazard assessment at Somma–Vesuvius based on the geological record. Bull. Volcanol. 72, 1021–1038. Giordano, D., Nichols, A.R., Dingwell, D.B., 2005. Glass transition temperatures of natural hydrous melts: a relationship with shear viscosity and implications for the welding process. J. Volcanol. Geotherm. Res. 142 (1), 105–118. Giordano, D., Russell, J.K., Dingwell, D.B., 2008. Viscosity of magmatic liquids: a model. Earth Planet. Sci. Lett. 271, 123–134. https://doi.org/10.1016/j.epsl. 2008.03.038.

555

Gonnermann, H., Taisne, B., 2015. Magma Transport in Dikes. In: The Encyclopedia of Volcanoes, second edition, pp. 215–224. Gudmundsson, A., 2012. Magma chambers: formation, local stresses, excess pressures, and compartments. J. Volcanol. Geotherm. Res. 237–238, 19–41. Hanson, J., 2016. Conduit Evolution During the Avellino Eruption at Vesuvius Based on Lithic Abundance and Componentry. PhD Thesis. University of Bristol, Bristol, UK. Hulme, G., 1982. A review of lava flow processes related to the formation of lunar sinuous rilles. Geophys. Surv. 5, 245–279. Johnston-Lavis, H.J., 1884. The geology of the Mt. Somma and Vesuvius: being a study of volcanology. Q. J. Geol. Soc. Lond. 40 (35), 149. Jurado-Chichay, Z., Walker, G.P.L., 2001. Variability of plinian fall deposits: examples from Okataina volcanic centre, New Zealand. J. Volcanol. Geotherm. Res. 111, 239–263. Kusumoto, S., Geshi, N., Gudmundsson, A., 2013. Aspect ratios and magma overpressures of non-feeder dikes observed in the Miyake-jima volcano (Japan), and fracture toughness of its upper part. Geophys. Res. Lett. 40, 1065–1068. Macedonio, G., Dobran, F., Neri, A., 1994. Erosion processes in volcanic conduits and application to the AD79 eruption of Vesuvius. Earth Planet. Sci. Lett. 121, 137–152. Melnik, O., Sparks, R.S.J., 2005. Controls on conduit magma flow dynamics during lava dome building eruptions. J. Geophys. Res. 110 (B2). https://doi.org/10.1029/ 2004JB003183. Melnik, O., Sparks, R.S.J., Costa, A., Barmin, A.A., 2008. Volcanic Eruptions: Cyclicity during Lava Dome Growth. In: Meyers: Encyclopedia of Complexity and Systems Science. Melnik, O., Costa, A., 2014. Dual-chamber–conduit models of non-linear dynamics behavior at Soufriere Hills Volcano, Montserrat. In: The Eruption of Soufriere Hills Volcano, Montserrat from 2000 to 2010, vol. 39. Geological Society, London Memories, pp. 61–69. Pappalardo, L., Mastrolorenzo, G., 2010. Short residence times for alkaline Vesuvius magmas in a multi-depth supply system: evidence from geochemical and textural studies. Earth Planet. Sci. Lett. 296, 133–143. Pyle, D.M., 1989. The thickness, volume and grainsize of tephra fall deposits. Bull. Volcanol. 51, 1–15. Rivalta, E., Taisne, B., Bunger, A.P., Katz, R.F., 2015. A review of mechanical models of dike propagation: schools of thought, results and future directions. Tectonophysics 638, 1–42. Santacroce, R., 1987. Somma Vesuvius. In: CNR Quaderni Ricerca Sci., vol. 114, p. 251. Santacroce, R., Cioni, R., Marianelli, P., Sbrana, A., Sulpizio, R., Zanchetta, G., Donahue, D.J., Joron, J.-L., 2008. Age and whole rock-glass composition of proximal pyroclastics from the major explosive eruptions of Somma–Vesuvius: a review as a tool for distal tephrostratigraphy. J. Volcanol. Geotherm. Res. 177, 1–18. Scaillet, B., Pichavant, M., Cioni, R., 2008. Upward migration of Vesuvius magma chamber over the past 20,000 years. Nature 45, 11. https://doi.org/10.1038/ nature07232. Sevink, J., van Bergen, M.J., van der Plicht, J., Feiken, H., Anastasia, C., Huizinga, A., 2011. Robust date fro the Bronze Age Avellino eruption (Somma–Vesuvius): 3945 ± 10 cal BP(1995 ± 10 cal BC). Quat. Sci. Rev. 30, 1035–1046. Sparks, R.S.J., 1986. The dimensions and dynamics of volcanic eruption columns. Bull. Volcanol. 48, 3–15. Sulpizio, R., Bonasia, R., Dellino, P., Mele, D., Di Vito, M.A., La Volpe, L., 2010. The Pomici di Avellino eruption of Somma–Vesuvius (3.9 ka BP) part II: sedimentology and physical volcanology of pyroclastic density current deposits. Bull. Volcanol 72, 559–577. https://doi.org/10.1007/s00445-009-0340-4. Sulpizio, R., Massaro, S., 2017. Influence of stress field changes on eruption initiation and dynamics: a review. Front. Earth Sci. 5, 18. Touloukian, Y.S., Judd, W.R., Roy, R.F., 1989. Physical Properties of Rocks and Minerals, vol. 548. Hemisphere, New York. Varekamp, J.C., 1993. Some remarks on volcanic vent evolution during plinian eruptions. J. Volcanol. Geotherm. Res. 54, 309–318. Wang, Z., 2000. Dynamic versus static elastic properties of reservoir rocks. In: SEG Books, vol. 19, pp. 531–539.

Further reading Macedonio, G., Neri, A., Martì, J., Folch, A., 2005. Temporal evolution of flow conditions in sustained magmatic explosive eruptions. J. Volcanol. Geotherm. Res. 143, 153–172. Muskhelishvili, N., 1963. Some Basic Problems in the Mathematical Theory of Elasticity. Noordhof, Leiden, The Netherlands. Sparks, R.S.J., 1978. The dynamics of bubble formation and growth in magmas: a review and analysis. J. Volcanol. Geotherm. Res. 3, 1–37.