Understanding the plume dynamics of explosive super-eruptions

Corrected: Author correction ARTICLE DOI: 10.1038/s41467-018-02901-0 OPEN Understanding the plume dynamics of explosive super-eruptions 1234567890():,; Antonio Costa 1, Yujiro J. Suzuki2 & Takehiro Koyaguchi2 Explosive super-eruptions can erupt up to thousands of km3 of magma with extremely high mass flow rates (MFR). The plume dynamics of these super-eruptions are still poorly understood. To understand the processes operating in these plumes we used a fluiddynamical model to simulate what happens at a range of MFR, from values generating intense Plinian columns, as did the 1991 Pinatubo eruption, to upper end-members resulting in coignimbrite plumes like Toba super-eruption. Here, we show that simple extrapolations of integral models for Plinian columns to those of super-eruption plumes are not valid and their dynamics diverge from current ideas of how volcanic plumes operate. The different regimes of air entrainment lead to different shaped plumes. For the upper end-members can generate local up-lifts above the main plume (over-plumes). These over-plumes can extend up to the mesosphere. Injecting volatiles into such heights would amplify their impact on Earth climate and ecosystems. 1 Istituto Nazionale di Geofisica e Vulcanologia, Bologna 40128, Italy. 2 Earthquake Research Institute, University of Tokyo, 1-1-1 Yayoi, Bunkyo-ku, Tokyo 1130032, Japan. Correspondence and requests for materials should be addressed to A.C. (email: [email protected]) NATURE COMMUNICATIONS | (2018)9:654 | DOI: 10.1038/s41467-018-02901-0 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-02901-0 E xplosive super-eruptions eject from several hundreds to thousands of km3 of magma at extremely high flow rates1. Many of these eruptions have had significant impacts to the climate and ecosystems2–5. Explosive super-eruptions cover areas within hundreds km from the vent with thick pyroclastic flows, blanket continent-size regions with ash, and inject large quantities of aerosols into the atmosphere4. Volatiles injected into the stratosphere can alter the Earth climate on a global scale even causing a volcanic winter that can persist for years to decades2,3. On the other hand, tephra layers associated with these catastrophic events are invaluable chronological markers across the affected regions6,7. The mass erupted during super-eruptions is orders of magnitude larger than the biggest eruptions experienced in historic times1,4,8–10. Estimates of mass flow rates (MFRs) during these super-eruptions, obtained from different independent approaches, suggest that they are extremely high, ranging from 109 to 1011 kg/s4,8,10–12. Such large MFRs require multiple vents or continuous emission along dykes13,14. Plume dynamics of explosive super-eruptions are not well understood as such large events have not been witnessed. In order to understand how such volumes of material are ejected and dispersed we rely on field evidence and models that can produce the observed deposits. Our current understanding on how plumes of explosive super-eruptions behave is from extrapolations of simple integral models developed for describing columns generated from small MFR. These simple models15–17 are based on the Buoyant Plume Theory (BPT) but the similarity assumption behind has been shown not to be valid for large MFRs18. Large explosive eruptions produce Plinian columns when the erupted mixtures of fragmented hot magma and gas entrain air, which heats up and expands making the plume buoyant. Above a critical MFR the eruption column becomes unstable19 and collapses, producing pyroclastic flows that spreads laterally on the ground. At high MFR, the dilute parts of the hot pyroclastic flows can also become buoyant as they also entrain air, forming a coignimbrite eruption plume that can rise up to the stratosphere carrying massive quantities of elutriated fine ash and volatiles. Results Fluid-dynamical regimes of eruptive plumes for large MFRs. Here, in order to avoid making unrealistic assumptions, we investigate the plume dynamics using a three-dimensional computational fluid-dynamical code (see Methods) designed to describe the evolution of volcanic plumes and umbrella clouds19. The code simulates the injection of a well-coupled mixture of solid pyroclasts (ash) and volcanic gas (assumed to be water vapour) from vents of different shapes above a flat surface into a stratified atmosphere. The model does not consider particle sedimentation and particle decoupling20 but captures the plume dynamics (see Methods). In this study we will not consider the effects of the rotation of the Earth on the plume dynamics, which can be very significant for very large eruptions, affecting, among other things, their spreading and shape of the plume21,22. For this reason and the for the sake of simplicity, here we focus on eruptions occurring in the equatorial belt, where these effects are negligible4,23 and consider tropical windless atmospheric conditions only (see Methods). Considering the input parameters reported in Table 1 and atmospheric properties described in Methods, we explored the effects of variable MFRs for different vent geometries, such as a single circular vent, fissure, and vents at different distances. For the sake of simplicity, we focus on the results from a circular vent but these are rather general. The fluid dynamics of large Plinian columns fed by a MFR~109 kg/s have been described in several studies19,23,24. In these columns, the fountain-like structure (radially suspended flow24,25) generated 2 NATURE COMMUNICATIONS | (2018)9:654 Table 1 Common input parameters and constants used for the three-dimensional simulations Variable Exit velocity Exit temperature Exit water fraction Exit density Gas constant of volcanic gas Gas constant of atmospheric air Specific heat of solid pyroclasts Specific heat of volcanic gas at constant volume Specific heat of air at constant volume Gravity body force Value 256 m/s 1053 K 0.06 3.5 kg/m3 462 J/kg/K 287 J/kg/K 1100 J/kg/K 1348 J/kg/K 717 J/kg/K 9.81 m/s2 in the lower part of the column is characterised by a highconcentration of erupted mixture and it is denser than the ambient air. In this region, the erupted mixture mixes with the air in largescale vortexes and this mixture becomes rapidly buoyant (see Fig. 1a and Supplementary Movie 1). Approaching the critical MFR at ~109.5 kg/s (see Methods), the radially suspended flow24,25 becomes unstable producing partial collapses, but the main plume still survives. The lower central part of the plume is a Negatively Buoyant Region (NBR), while the area around it is fed with relatively pure air that maintains its buoyancy and efficiently transports the mixture up to the stratosphere (see Fig. 1b and Supplementary Movie 2). The highest velocity remains in the central region generating a mushroom shape plume, with very high mass fractions in the central part up to the top of the plume (see Fig. 1b and Supplementary Movie 2). Increasing MFR up to 1010 kg/s produces a total collapse of the radially suspended flow, which generates continuous fountaining to the ground, feeding pyroclastic density currents and increasing the radius of the hotter NBR, resulting in a basal region (~60 km diameter) from where the large co-ignimbrite plume will develop (see Fig. 1c). In this case, because of the vigorous rising velocities in the periphery owing to the more effective local air entrainment, the upper central portion of such a large plume has a relatively low mass fraction compared with the outer region (Fig. 1c and Supplementary Movie 3). The resulting plume still maintains a mushroom shape but the plume top has flat, rather than umbonate, cap (Fig. 1c and Supplementary Movie 3). A further increase of MFR up to 1011 kg/s will produce a larger co-ignimbrite plume (>150 km in diameter). In this case, the vortices entrain air mainly at the periphery of the co-ignimbrite plume, without affecting the area around the hotter NBR (Fig. 1d and Supplementary Movies 4 and 5). This allows the periphery region of the co-ignimbrite plume to become much more buoyant and increase its velocity. Because of mass conservation, the vertical velocity in the inner part of the plume decreases. This regime results in the formation of a sort of toroid umbrella (donut-like shape), giving to the plume a depressed-cap mushroom shape (i.e., two separate lobes in a 2D cross-section). Our simulations show that vent geometry has a strong control on the dynamics and stability of the Plinian columns but once coignimbrite plumes are generated the processes are predominantly controlled by the diameter of the co-ignimbrite plume, which is typically larger that the vent area, i.e., longer than the fissures or the distance between multiple vents (Supplementary Movies 1–6). For these reasons, the described features apply to all vent geometries despite the fact that a circular vent was used for the simulations (see Supplementary Movie 6). Implications for the assessment of eruption parameters. MFRs are typically estimated from the total plume heights assuming BPT is valid. A similar approach has been extended to co| DOI: 10.1038/s41467-018-02901-0 | www.nature.com/naturecommunications ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-02901-0 b MFR = 109 kg/s MFR = 109.5 kg/s 50 50 40 40 Z (km) Z (km) a 30 20 30 20 10 10 0 −100 −50 0 50 0 −100 100 −50 0 50 0.0 0.2 0.4 0.6 0.8 1.0 mfmgm(−) 50 40 40 Z (km) Z (km) 0.0 0.2 0.4 0.6 0.8 1.0 mfmgm(−) 50 30 20 30 20 10 10 0 −100 −50 0 50 0 −100 100 −50 0 X (km) 0.5 1.0 50 50 40 40 30 20 –40 0 X (km) 40 80 –80 –40 80 40 0 1.0 20 Y (km) d MFR = 1010 kg/s 40 40 Z (km) 50 30 20 –40 0 X (km) 80 –80 40 –40 80 40 0 Y (km) MFR = 1011 kg/s 50 30 20 10 0 −100 −50 0 50 0 −100 100 −50 0 50 100 X (km) X (km) 0.0 0.2 0.4 0.6 0.8 1.0 mfmgm(−) 0.0 0.2 0.4 0.6 0.8 1.0 mfmgm(−) 50 50 40 40 Z (km) Z (km) 0.5 30 0 –80 10 30 20 30 20 10 10 0 −100 −50 0 50 0 −100 100 −50 0 −1.0 −0.5 0.0 ∆ 0.5 −1.0 1.0 50 40 40 Z (km) 50 30 20 10 –40 0 40 50 100 X (km) X (km) Z (km) 0.0 ∆ 10 0 –80 0 –80 100 −1.0 −0.5 Z (km) Z (km) 0.0 ∆ 10 Z (km) 50 X (km) −1.0 −0.5 c 100 X (km) X (km) 80 –80 –40 0 40 80 Y (km) X (km) −0.5 0.0 ∆ 0.5 1.0 30 20 10 0 –120 –60 0 X (km) 60 120 –120 –60 60 0 Y (km) 120 Fig. 1 Simulation results for the plumes with different MFRs. a 109 kg/s, b 109.5 kg/s, c 1010 kg/s, d 1011 kg/s. The figures show snapshots at t = 800 s after the eruption initiation: vertical cross-sections of the mass fraction of the erupted mixture is 0.01 (a–d, upper panels); density difference relative to the atmospheric density (a–d, middle panels); three-dimensional isosurface where the mass fraction of the erupted mixture is 0.01 (a–d, lower panels) ignimbrite plume15 and is still largely used in the volcanological community16,17. However, the dynamics of large co-ignimbrite plumes are markedly different from BPT as their horizontal extension is typically much larger than their height. The simulations indicate that, above a critical MFR, the co-ignimbrite plume rises from a source pffiffiffiwith ffiffiffiffiffiffiffi radius RCI , increasing with MFR as by RCI
2:8  104 MFR (with the co-ignimbrite plume radius, RCI , expressed in km and MFR in kg/s, see Fig. 2); this NATURE COMMUNICATIONS | (2018)9:654 power-law dependence of the MFR with run-out distance was predicted in simple models of pyroclastic flows26. The mechanisms of air entrainment from such broad sources are profoundly different from, and invalidate the similarity assumption used in, BPT. The difference in the scale of the horizontal extension of plume also affects the behaviour of the upper part of the plume, including the column height and the dynamics in the umbrella cloud. | DOI: 10.1038/s41467-018-02901-0 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-02901-0 100 1000 Radius (km) RCI (km) 70 30 10 9 10 1010 MFR (kg/s) 70 Ht 50 7 0.0 Height (km) 0.10 40 0.20 0.30 30 0.12 0.30 20 HMSL 0.10 0.20 0.30 1.00 HNBL 10 0 109 1010 MFR (kg/s) 1011 Fig. 3 Characteristic plume heights as a function of the MFR. Red bars represent the height over the whole plume (Ht). Blue squares represent the altitude of maximum level of radial spreading of the umbrella cloud (HMSL). Green bars represent the altitude of NBL (HNBL). The characteristic heights from the BPT model are shown as dashed and solid curves. In this diagram, the total heights are compared with the heights where the upward velocity becomes zero in the BPT calculations, whereas the MSLs and NBLs are compared with the NBL calculated by the BPT model (see also Table 2). Values on the curves are the effective entrainment coefficients that explain the 3D results. In the calculations of co-ignimbrite ash clouds (the MFR > 1010 kg/s) using the BPT model, the initial temperature and mass fraction of the entrained air were estimated so that the initial density is equal to the atmospheric density near the ground (see Table 3). According to this mass fraction, the effective MFR at the ground can be estimated. The initial radii were calculated using the scaling law reported in Fig. 2 Here, we focus on the dependence of maximum plume height with MFR (Fig. 3 and Fig. 4 for the dynamics of umbrella cloud). The result shows that from 109 to 1010 kg/s the maximum plume heights remain similar and are between 40 and 60 km. Coignimbrite plumes appear steadier than Plinian columns, which show a more oscillatory behaviour (see Supplementary Movies 1–5). 4 100 t 2/3 10 1011 Fig. 2 Co-ignimbrite plume radius as function of MFR. Cross denotes values estimated from the 3D simulations and the line the empirical relationship given by RCI
a  MFRb with a ffi 2:8  104 km kgb sb and b ffi 0:5. Note that the model simulates a dusty-gas mixture, whereas the pyroclastic flows from such super-eruptions can be formed of concentrated dispersions, which depend on imperfect coupling between gas and particles10 60 109 kg/s 109.5 kg/s 1010 kg/s 1010.5 kg/s 1011 kg/s NATURE COMMUNICATIONS | (2018)9:654 1 100 1000 5000 Time (s) Fig. 4 Radii of the spreading umbrella clouds as a function of the time for the different simulations. Dashed lines correspond to the relationships between radius, R, and time, t: RðtÞ ¼ At2=3 , where A is the proportionality constant. This dependence is commonly used to estimate MFR23,35,36. The plot shows that with respect to the Plinian case, co-ignimbrite plumes with MFR > 1010 kg/s have a steeper spreading rate, with the power-law exponent23,35,36 passing from ~0.66 to 0.96. These results validate the robustness of simple box models for describing the umbrella region35,36. However, the fact that the spreading of the umbrella cloud with time is faster for co-ignimbrite plumes with respect Plinian plumes suggests caution when we use relationships obtained for the latter to estimate MFR above the critical value The maximum column height and highest mass fraction in the umbrella region is reached for MFR around the critical value, i.e., 109.5 kg/s. Remarkably, the Neutral Buoyancy Level (NBL) remains at ~20 km for all the simulations with MFRs above the critical MFR (see Table 2). For MFR of ~1011 kg/s, maximum plume height is in the peripheral region rather than in the centre, owing to the more efficient entrainment of air from the border of the plume. The maximum height for the bulk mass is at ~50 km but local up-lifts, having a diameter of ~30–40 km, develop above the umbrella region (see Supplementary Figure 1) and keep rising up to the mesosphere (60–70 km). In this case, two different effective heights should be considered, one for the bulk mass spreading around the umbrella region and one for the maximum height reached by the local up-lifts (we call them “local over-plume” hereafter). These local over-plumes develop from the base of the periphery of the coignimbrite plumes, because of local heterogeneities in the efficiency of air entrainment, and are characterised by higher velocities and larger mass fractions (see Supplementary Figure 1). The complex relationship between plume height and MFR in Fig. 3 suggests that for large eruption intensities (MFR > 109 kg/s) we cannot use column height estimations to assess the value of MFR15,27. To compare our results from the three-dimensional simulations with those of the simple BPT integral models15,27, we estimated the main mean variables20, such as mixture density difference, Δρ, vertical velocity, U, temperature, T, and mass fraction ξ (see Supplementary Figure 2) and extracted optimal parameter values (see Fig. 3 and Table 3). These BPT models do not adequately describe co-ignimbrite plumes but if they are used as extrapolations the effective entrainment coefficient, k, should be properly tuned and not assumed as an invariant. Since variations of HNBL with MFR from sustained Plinian column to fully co-ignimbrite plume are almost negligible (HNBL ~15–20 km, see Table 2), accordingly to BPT, this implies that the effective entrainment coefficient should increase with MFR and have significantly different values for the two regimes. | DOI: 10.1038/s41467-018-02901-0 | www.nature.com/naturecommunications ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-02901-0 Table 2 Typical plume heights for different MFRs Plume dynamics MFR (kg/s) 109 12–19 23 40–45 Low HNBL ðkmÞa HMSL ðkmÞb Ht (km)c Mean ξ in the umbrella a b c 109.5 19–23 25 50–60 Very high 1010 15–18 22 40 Medium 1010.5 16–20 26 50–60 High 1011 18–21 29 60–70 Very high The lower value corresponds to the column region, while the upper value to the entire umbrella region The radial maximum spreading level was estimated as the height containing the greatest mass of airborne particles34 The lower value corresponds to dilute part of the plume (ξ = 10−3), while the upper value to the concentrated part of the plume (ξ = 10−1) Table 3 Input parameters for BPT model inferred from three-dimensional simulations Parameter Plume radius (km) Exit velocity (m/s) Exit temperature (K) Exit gas fraction MFR (kg/s) 109 0.60 256 1053 0.06 109.5 1.06 256 1053 0.06 1010 28.0 4.87 858 0.33 1010.5 50.0 4.87 858 0.33 1011 88.5 4.87 858 0.33 The total column height, Ht, varies with MFR (see Table 2), rising from Ht ffi 40  45 km at MFR ¼ 109 kg/s to Ht ffi 50  60 km at MFR ¼ 109:5 kg/s, then decreasing Ht ffi 40 km at MFR ¼ 1010 kg/s and increasing again to Ht ffi 50  60 km at MFR ¼ 1010:5 kg/s. To estimate plume height ranges, we considered a fraction mass of ξ ¼ 103 for the dilute upper region and ξ ¼ 101 for the lower values, which refer to a more concentrated plume region. This implies that the behaviour of the plumes above NBL and the plume height is controlled by an effective entrainment coefficient, kU, which is different from the entrainment coefficient, kL, that governs the air entrainment below NBL (see Fig. 3). The total heights obtained from threedimensional simulations cannot be described by the BPT model (for the optimal input values see Table 3) even if the value of kU is empirically tuned. The maximum height is reached by the oscillating local over-plume at MFR > 1010.5 kg/s. It is inferred that the total column height (reported in Fig. 3) is generated from a large amplitude gravity wave that is excited by the intensive plume. This again confirms the difficulty of simple integral models to capture such complex plume dynamics. Discussion Our results have enormous implications for the assessment of the dynamics of super-eruptions. For the most extreme MFRs13, ~1011 kg/s, similar to those estimated for the Young Toba Tuff (YTT) eruption4, the total plume height increases beyond the stratosphere (up to ~60–70 km) owing to the development of local over-plumes above the umbrella region, even though the NBL remains at ~20 km. The bulk mass spreads in the umbrella region between ~20 and ~50 km, and the local over-plumes develop above the main umbrella region (see Supplementary Figure 1). All previous studies relied on the results of the Woods and Wohletz model15 and used a relatively low plume height (~30–40 km) but our simulations indicate that volatiles and fine ash can be transported up to 70 km. This has important implications for the effects of eruptions, like YTT, on the Earth’s climate2,3,28. Climate models are sensitive to the injection height because atmosphere stratification and availability of H2O influence the conversion of SO2 into sulphate, which governs the climate response. Injections into the stratosphere affect the albedo of the atmosphere on the order of decades2 but the longevity of SO2 in the mesosphere could be considerably longer. NATURE COMMUNICATIONS | (2018)9:654 Our simulations also show that partially collapsing plumes, generated slightly above the critical MFR value, can reach heights of up to 50–60 km and efficiently transport the mixture up to the stratosphere (see Fig. 1b and Supplementary Movie 2). This results in effective upward transport of a large mass fraction of fine ash that is generated from the pyroclastic flows, enriching the fine ash content of the umbrella cloud with respect to the groundhugging pyroclastic flows. This can be the case for eruptions similar to the Campanian Ignimbrite event for which a MFR of ~109 kg/s was empirically estimated for the initial Plinian phase and ~2–5×109 kg/s for the co-ignimbrite phase12,29. Methods Computational fluid-dynamical model. To simulate fluid dynamics of volcanic plumes, we used the pseudo gas model by Suzuki et al.19, in which the momentum and heat exchanges between the volcanic ash and gas are assumed to be so rapid that the velocity and temperature are the same for all phases. The fluid dynamics model solves a set of partial differential equations describing the conservation of mass, momentum, and energy, and a set of constitutive equations describing the thermodynamics state of the mixture of volcanic ash, gas, and air. The computational model was designed to reproduce the injection of a mixture of volcanic ash and volcanic gas from a vent into a stratified atmosphere. For this study typical tropical atmosphere conditions were used as initial atmospheric conditions (see Supplementary Figure 3). For given temperature gradients of atmosphere, the atmospheric density and pressure were calculated from the hydrostatic relationship. Initial wind velocity was set to be zero in the whole computational domain. Free slip conditions were applied at the ground boundary, and fixed inflow conditions at the vent. At the other boundaries, mass, momentum, and energy fluxes were assumed to be continuous (i.e., free outflow/inflow conditions). The use of the latter conditions can avoid the wave reflection at the boundaries even when the initial atmospheric conditions, such as hydrostatic pressure and density, are perturbed by the eruption. The governing equations were solved numerically by the Roe scheme30 with MUSCL (Monotone Upstream-centred Scheme for Conservation Laws) interpolation31 for spatial integration and time splitting method for time integration. The present model with these schemes has a third-order accuracy in space and second-order accuracy in time. We used a generalised coordinate system in the computational domain. Computational settings and simulation strategy. The grid size near the vent was set to be sufficiently smaller than the vent diameter (i.e., D0/20, where D0 is the vent diameter) in order to resolve the flow structures and turbulent mixing at a low altitude, whereas it was increased at a constant rate up to 300 m/grid with the distance from the vent. We confirmed that the flow regime and plume height are basically unaffected by the numerical resolutions by changing the grid sizes (see Supplementary Figure 4). The computations were carried out on the Earth Simulator (NEC SX-ACE; 64 GFLOPS/core) at the JAMSTEC and also on the Fujitsu PRIMERGY CX400 (22 GFLOPS/core) at the Research Institute, Kyushu University. Each simulation took 52–280 h with 512 cores. In order to investigate how variation of MFR affects the plume dynamics, we carried out a set of simulations with variable vent size, following the method described by Suzuki et al.25. The vent radius ranged from 600 to 6000 m. The pressure, density, and velocity at the vent were kept fixed among the simulations. The pressure at the vent was assumed to be in equilibrium with the atmospheric pressure. We assumed a magmatic temperature of 1053 K and water content of 6 wt.%, similar to the magmatic properties in the Pinatubo 1991 eruption23. Considering these values of pressure, temperature, and density the equation of state gives a density at the vent of 3.5 kg/m3. The exit velocity was assumed to be 256 m/s, corresponding to a fixed Mach number of 1.5. As a result, the MFR ranges from 109 to 1011 kg/s. It is important to stress that in this study we kept fixed the magmatic properties such as magma temperature and water content. The variation of these properties | DOI: 10.1038/s41467-018-02901-0 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-02901-0 can lead to different critical conditions between the flow regimes described in the main text. However, the qualitative features of each flow regimes would be same in the possible ranges of magmatic temperature and water content for magmatic eruptions. We also assumed the equilibrium of pressure and the supersonic flow as the conditions at the vent. The disequilibrium and sub/supersonic flow can change the flow structures near the vent, which results in the change of the final distance of PDC and therefore the change of transition between the flow regimes32,33. In addition, we ignore the non-equilibrium effects between the volcanic ash and gas phases. However, the non-equilibrium effects are less relevant in the strong eruptions rather than in the weak eruptions20. Data availability. The authors declare that all data supporting the findings of this study are available in the article and in Supplementary Information. Additional information is available from the corresponding author upon request. Received: 29 July 2017 Accepted: 8 January 2018 References 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 6 Self, S. The effects and consequences of very large explosive volcanic eruptions. Philos. Trans. R. Soc. A 364, 2073–2097 (2006). Robock, A. et al. Did the Toba volcanic eruption of ~74 ka B.P. produce widespread glaciation? J. Geophys. Res. 114, D10107 (2009). Timmreck, C. et al. Climate response to the Toba super-eruption: regional changes. Quat. Int. 258, 30–44 (2012). Costa, A., Smith, V., Macedonio, G. & Matthews, N. The magnitude and impact of the Youngest Toba Tuff super-eruption. Front. Earth Sci. 2, 16 (2014). Black, B. A., Neely, R. R. & Manga, M. Campanian Ignimbrite volcanism, climate, and the final decline of the Neanderthals. Geology 43, 411–414 (2015). Lowe, J. J. et al. Volcanic ash layers illuminate the resilience of Neanderthals and early modern humans to natural hazards. Proc. Natl. Acad. Sci. USA 109, 13532–13537 (2012). Smith, V. C., Isaia, R., Engwell, S. & Albert, P. G. Tephra dispersal during the Campanian Ignimbrite (Italy) eruption: implications for ultra-distal ash transport during the large caldera-forming eruption. Bull. Volcanol. 78, 45 (2016). Wilson, C. J. N. & Hildreth, W. The Bishop Tuff: new insights from eruptive stratigraphy. J. Geol. 105, 407–439 (1997). Mastin, L. G., Van Eaton A. R. & Lowenstern J. B., Modeling ash fall distribution from a Yellowstone supereruption. Geochem. Geophys. Geosyst. 15, 3459–3475 (2014). Roche, O., Buesch, D. C. & Valentine, G. A. Slow-moving and far- travelled dense pyroclastic flows during the Peach Spring super-eruption. Nat. Commun. 7, 10890 (2016). Wilson, C. J. N. & Walker, G. P. L. in Tephra Studies, Proceedings of the NATO Advanced Study Institutes Series, vol. 75 (eds Self, S. & Sparks, R. S. J.) 441–448 (Springer, Dordrecht, 1981). Marti, A., Folch, A., Costa, A. & Engwell, A. Reconstructing the plinian and co-ignimbrite sources of large volcanic eruptions: a novel approach for the Campanian Ignimbrite. Sci. Rep. 6, 21220 (2016). Costa, A. & Marti, J. Stress field control during large caldera-forming eruptions. Front. Earth Sci. 4, 92 (2016). Costa, A., Gottsmann, J., Melnik, O. & Sparks, R. S. J. A stress-controlled mechanism for the intensity of very large magnitude explosive eruptions. Earth. Planet. Sci. Lett. 310, 161–166 (2011). Woods, A. W. & Wohletz, K. The dimensions and dynamics of coignimbrite eruption columns. Nature 350, 225–227 (1991). Engwell, S. L., de’ Michieli Vitturi, M., Esposti Ongaro, T. & Neri, A. Insights into the formation and dynamics of coignimbrite plumes from onedimensional models. J. Geophys. Res. 121, 4211–4231 (2016). Jessop, D. E., Gilchrist, J., Jellinek, A. M. & Roche, O. Are eruptions from linear fissures and caldera ring dykes more likely to produce pyroclastic flows? Earth. Planet. Sci. Lett. 454, 142–153 (2016). Costa, A. et al. Results of the eruption column model inter-comparison study. J. Volcanol. Geotherm. Res. 326, 2–25 (2016). Suzuki, Y. J., Koyaguchi, T., Ogawa, M. & Hachisu, I. A numerical study of turbulent mixing in eruption clouds using a three-dimensional fluid dynamics model. J. Geophys. Res. 110, B08201 (2005). Cerminara, M., Esposti Ongaro, T. & Neri, A. Large eddy simulation of gas–particle kinematic decoupling and turbulent entrainment in volcanic plumes. J. Volcanol. Geotherm. Res. 326, 143–171 (2016). Baines, P. & Sparks, R. S. J. Dynamics of giant volcanic ash clouds from supervolcanic eruptions. Geophys. Res. Lett. 32, L24808 (2005). Baines, P. G., Jones, M. T. & Sparks, R. S. J. The variation of large- magnitude volcanic ash cloud formation with source latitude. J. Geophys. Res. 113, D21204 (2008). NATURE COMMUNICATIONS | (2018)9:654 23. Suzuki, Y. J. & Koyaguchi, T. A three-dimensional numerical simulation of spreading umbrella clouds. J. Geophys. Res. 114, B03209 (2009). 24. Neri, A. & Dobran, F. Influence of eruption parameters on the thermofluid dynamics of collapsing volcanic columns. J. Geophys. Res. 99, 11833–11857 (1994). 25. Suzuki, Y. J., Costa, A. & Koyaguchi, T. On the relationship between eruption intensity and volcanic plume height: insights from three-dimensional numerical simulations. J. Volcanol. Geotherm. Res. 326, 120–126 (2016). 26. Bursik, M. I. & Woods, A. W. The dynamics and thermodynamics of large ash flows. Bull. Volcanol. 58, 175–193 (1996). 27. Sparks, R. S. J. The dimensions and dynamics of volcanic eruption columns. Bull. Volcanol. 48, 3–15 (1986). 28. Rampino, M. R. & Self, S. Volcanic winter and accelerated gla- ciation following the Toba super-eruption. Nature 359, 50–52 (1992). 29. Costa, A. et al. Quantifying volcanic ash dispersal and impact of the Campanian Ignimbrite super-eruption. Geophys. Res. Lett. 39, L10310 (2012). 30. Roe, P. L. Approximate Riemann solvers, parameter vectors, and difference schemes. J. Comput. Phys. 43, 357–372 (1981). 31. van Leer, B. Towards the ultimate conservative difference scheme III. Upstream-centered difference schemes for ideal compressible flow. J. Comput. Phys. 23, 263–275 (1977). 32. Ogden, D., Glatzmaier, G. A. & Wohletz, K. H. Effects of vent overpressure on buoyant eruption columns: implications for plume stability. Earth Planet. Sci. Lett. 268, 283–292 (2008). 33. Carcano, S., Bonaventura, L., Esposti Ongaro, T. & Neri, A. A semi-implicit, second-order-accurate numerical model for multiphase underexpanded volcanic jets. Geosci. Model. Dev. 6, 1905–1924 (2013). 34. Suzuki, Y. J. et al. Inter-comparison of three-dimensional models of volcanic plumes. J. Volcanol. Geotherm. Res. 326, 26–42 (2016b). 35. Woods, A. & Kienle, J. The dynamics and thermodynamics of volcanic clouds: theory and observations from the April 15 and April 21, 1990 eruptions of Redoubt Volcano, Alaska. J. Volcanol. Geotherm. Res. 62, 273–299 (1994). 36. Costa, A., Folch, A. & Macedonio, G. Density-driven transport in the umbrella region of volcanic clouds: implications for tephra dispersion models. Geophys. Res. Lett. 40, 4823–4827 (2013). Acknowledgements A.C. was partially supported by a grant of the International Research Promotion Office Earthquake Research Institute, the University of Tokyo. Y.J.S. and T.K. were partially supported by KAKENHI (grand nos. 25750142 and 17K01323). We warmly thank V.C. Smith for revising the English and very helpful suggestions. Author contributions A.C. and Y.J.S. designed the simulation set. Y.J.S. performed the runs and processed the results. A.C., Y.J.S., and T.K. analysed and interpreted the results. A.C. wrote the manuscript with input from all the co-authors. Additional information Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467018-02901-0. Competing interests: The authors declare no competing financial interests. Reprints and permission information is available online at http://npg.nature.com/ reprintsandpermissions/ Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. © The Author(s) 2018 | DOI: 10.1038/s41467-018-02901-0 | www.nature.com/naturecommunications 

DOI: 10.1038/s41467-018-05545-2 OPEN Author Correction: Understanding the plume dynamics of explosive super-eruptions Antonio Costa 1, Yujiro J. Suzuki2 & Takehiro Koyaguchi2 We became aware of a mistake in the data displayed in the original version of Fig. 3. Specifically, the lines showing the relationship between column height and MFR for MFR larger than 1010 kg/s were based on simulations in which the exit gas fraction was assumed to be an unrealistic value of 0.76 rather than the correct value of 0.33. The correct version of Fig. 3 is: 70 Ht 60 50 7 0.0 0.10 Height (km) 1234567890():,; Correction to: Nature Communications https://doi.org/10.1038/s41467-018-02901-0; published online 12 Feb 2018 40 0.20 0.30 30 0.12 0.30 20 HMSL 0.10 0.20 0.30 1.00 HNBL 10 0 109 1010 MFR (kg/s) 1011 Fig. 3 1 Istituto Nazionale di Geofisica e Vulcanologia, 40128 Bologna, Italy. 2 Earthquake Research Institute, University of Tokyo, 1-1-1 Yayoi, Bunkyo-ku, Tokyo 1130032, Japan. Correspondence and requests for materials should be addressed to A.C. (email: [email protected]) NATURE COMMUNICATIONS | (2018)9:3471 | DOI: 10.1038/s41467-018-05545-2 | www.nature.com/naturecommunications 1 AUTHOR CORRECTION NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-05545-2 which replaces the previous incorrect version: 70 0 0.1 Ht 60 Height (km) 50 7 0.0 40 30 H MSL 0.12 0.30 20 0.30 H NBL 10 0 0.20 109 1010 MFR (kg/s) 1011 This has been corrected in both the PDF and the HTML versions of the Article. The text was written on the basis of the correct plots, and so this error does not affect the original discussion or conclusions of the Article. The authors apologize for the confusion caused by this mistake. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. © The Author(s) 2018 2 NATURE COMMUNICATIONS | (2018)9:3471 | DOI: 10.1038/s41467-018-05545-2 | www.nature.com/naturecommunications