Large explosive eruptions may be dominated by pyroclastic flows instead of buoyant plumes: insights from a global data compilation

Abstract

The dynamics of Earth’s largest, most voluminous eruptions (≥ 100 km³ ejecta; VEI 7) are poorly understood. Here, we explore the question of whether these very large eruptions can be treated as scaled-up versions of moderate volume historical eruptions (0.1–10 km³ ejecta; VEI 4–6), or whether they should be treated as fundamentally different Earth system phenomena. To examine this, we compile fall deposit and ignimbrite volume data for 74 explosive eruptions worldwide that are magnitude 4 or greater, and use this dataset to assess how material is partitioned into buoyant plumes versus pyroclastic density currents as a function of eruption magnitude. Importantly, we filter our results by overall distance from seas/oceans, so that we can focus in on the eruptions for which preservation of the deposits is as reliable as possible. After filtering, we find that the largest eruptions are dominated by ignimbrites and not fall deposits, implying that, co-ignimbrite plumes notwithstanding, there may be little or no buoyant plume component to eruptions of the highest known magnitudes. This result is consistent with model simulations showing that the pyroclastic materials produced during larger events can be emplaced in density currents alone, and highlights important considerations for contemporary eruption simulations, the fate of volcanic gases relative to solid mass, and subsequent appraisals of the climatological and environmental impacts of explosive volcanism on Earth.

Introduction

The largest volcanic events in Earth history are complex, underrepresented in the geological record, and thus poorly understood (Mason et al. 2004; Miller and Wark 2008; Wilson et al. 2021). Explosive volcanic eruptions emit prodigious mass of solid material in the form of tephra and volatile gases such as SO₂ and H₂O (Carn et al. 2016; Jones 2015; Oppenheimer et al. 2011). Both SO₂ and H₂O can cause a sizeable disturbance to the Earth’s atmosphere and environment (Mather, 2015), and the potential for ocean-cryosphere-climate feedbacks resulting from volcanic aerosol forcing could even push these effects into millennia (Baldini et al. 2018, 2015). The very largest events (e.g., ≥ 100 km³ ejecta dense rock equivalent erupted volume; DRE) occur roughly once every ~ 20,000 years (Rougier et al. 2018), and no such event has been witnessed or studied directly. Furthermore, the most recent magnitude (M) ≥ 7 event of Mount Tambora in 1815 CE was not the subject of direct scientific observation. Although that eruption dispersed ~ 42 (± 4) km³ of material across the Indonesian archipelago (Kandlbauer and Sparks 2014), its dynamics and immediate impacts have been reconstructed entirely from historical accounts, geological evidence, and more recent analyses of its deposits and global climatic effects (e.g., Oppenheimer et al. 2003; Sigurdsson & Carey, 1992). Taken together, it is clear that the dynamics of exceptionally large eruptions are less constrained than for lower relative eruption magnitudes that are regularly observed.

Jets of gas and solid particles are emitted during an explosive eruption. These jets entrain and heat air so that they may become buoyant relative to the surrounding atmosphere, forming plumes (Sparks and Wilson 1982), which can generate two primary phenomena. First, if the jet becomes buoyant before losing its initial kinetic energy, the plume may continue to rise to form an umbrella cloud, transporting large amounts of ash (diameter < 2 mm) into the atmosphere, dispersing it and then depositing it as fall deposits up to hundreds of kilometres away from the source (Bonadonna and Costa 2013; Constantinescu et al. 2021). Second, if buoyancy is not attained prior to losing its initial kinetic energy the plume may collapse, generating pyroclastic density currents (PDCs), capable of blanketing large areas of the landscape surrounding the volcano with ignimbrite (Branney and Kolekaar 2002; Giordano and Cas 2021; Roche et al. 2022). The extent to which ignimbrite and fall deposits account for an eruption’s total volume will vary between eruptions owing to a combination of atmospheric and eruption factors. For example, plume height and stability will be influenced by atmospheric stratification and moisture content at the time of eruption (Bursik 2001; Sparks et al. 1997). Plume evolution may also be dictated by source magma composition, changes in vent size and shape, and mass eruption rate (e.g., Jessop et al., 2016; Woods and Wohletz 1991), among other factors.

Many ignimbrites occur in conjunction with one (or several) basal Plinian fall deposits (Baines and Sparks 2005; Jessop et al. 2016). However, the geological record appears to contain evidence for eruptions that are associated with only very small Plinian fall volumes or even no fall deposits at all; these are most commonly the eruptions responsible for producing some of the largest ignimbrites preserved on Earth (Cas et al. 2011; De Silva et al. 2006). For example, the ~ 2 Ma Cerro Galán Ignimbrite eruption (Argentina; Folkes et al. 2011), ~ 27 Ma Fish Canyon Tuff eruption of La Garita caldera (USA; Gleadow et al. 2015), and the ~ 164 ka Xáltipan eruption of Los Humeros caldera (Mexico; Cavazos-Álvarez and Carrasco-Núñez 2020). Where no fall deposits are recorded, this may be because the eruption did not form a buoyant plume at all; but this may also be because the deposits were eroded, buried, or reworked into subsequent density currents. Clearly then, detailed field work across proximal and distal distances from the eruptive vent are needed. Fluid-dynamical models have demonstrated that eruption plumes can become critically unstable under increasingly high mass eruption rates (MER; e.g., 10⁸ to 10¹¹ kg/s; Costa et al., 2018; Costa and Martí, 2016; Roche et al. 2021, 2016) and/or if erupted from especially large vents (Costa and Martí 2016; Jessop et al. 2016). This instability would favour ground-hugging pyroclastic currents as the dominant mode by which erupted material is distributed across the landscape and, therefore, show that increasingly large eruptions would produce increasingly large currents in response to increasing MER, and thus increasingly large ignimbrite volumes (Woods and Wohletz 1991).

Field, laboratory, and computer-based studies provide clear indications that fundamental differences exist between the eruptive behaviour of the largest explosive eruptions in the geological record, and the comparatively smaller events observed in historic times. For example, volcanoes that produce exceptionally large eruptions exhibit much greater diversity in their sub-surface architecture which can include deep, vertically extensive systems composed of multiple (compositionally distinct) magma bodies (e.g., Toba caldera; Pearce et al. 2020), or a single source of magma that is layered and/or mixed with different compositions (e.g., Aso caldera; Keller et al. 2021): both of which can be influenced by processes such as pre/syn-eruptive mixing, and/or rejuvenation and transport of crystals and melt (e.g., Taupo caldera; Allan et al. 2017). Hence, magmas that feed these eruptions can come from one of several storage regions, be zoned or convectively mixed, and have been assembled over timescales ranging from centuries, to many millions of years (Wilson et al. 2021).

Models also suggest that the mass flow rates associated with exceptionally large eruptions could range from 10⁸ to 10¹¹ kg/s – rates requiring repeat eruptions from multiple vents and/or continuous emission along dykes (Costa et al. 2018; Costa and Martí 2016; Roche et al. 2021, 2016). Many ignimbrites overlay several basal fall deposits in their host stratigraphic sequences, pointing to collapse of a vertical plume during the source eruption and/or an increase in eruption intensity (waxing); both indicative of a shift in eruptive dynamics (Baines and Sparks 2005; Jessop et al. 2016). However, some eruption stratigraphies exhibit a total, or near-total, absence of both proximal and distal fall deposits. On one hand, the absence of these fall deposits from the geological record may reflect post-depositional erosion and resuspension processes (e.g., Mastin et al. 2014), and/or suggest that fall was re-worked into the larger currents (Branney and Kolekaar 2002). However, it may also signal that these eruptions did not form stable and buoyant eruption columns capable of distal ash dispersal (Cas et al. 2011; De Silva et al. 2006).

Recent work suggests that the largest explosive events on Earth cannot be fully understood as "scaled-up" versions of smaller, historical events. By analysing the links between eruption size—defined by magnitude, ignimbrite, and fall deposit volumes—and plume behaviour, we can more closely examine how eruption dynamics vary based on different geophysical factors, and so provide critical information for assessment and modelling of the hazards posed by various volcanological phenomena (Barclay et al. 2019; Roche and Carazzo 2019). Here, we compile data from 74 eruptions (M ≥ 4 at > 0.1 km³ up to M ≤ 9.1 at > 10,000 km³) representing a wide range of eruption types, locations, and magma compositions, alongside information on ignimbrite and fall deposit volumes, eruption magnitude, and proximity to the nearest sea or ocean. Using this dataset, we investigate how eruption magnitude could have influenced volume partitioning across various scenarios.

Methods

Data acquisition

The database presented here was constructed with the objective of compiling data for the volume of ignimbrite, fall deposit, and total erupted volumes for a range of explosive eruptions worldwide. We chose ≥ 4 as the magnitude (M) threshold for defining ‘explosive eruptions’, which encapsulates eruptions with bulk ejecta volumes exceeding 0.1 km³. Eruptions of this size cause sizable environmental impacts (Bonadonna et al. 2016) and this threshold also aligns with established definitions of explosive volcanism (Crosweller et al. 2012; Papale 2018; Rougier et al. 2016).

Both direct (instrumental) and indirect (geological) evidence exists for eruptions that produced both effusive and explosive features; termed ‘hybrid’ eruptions (Bonadonna et al. 2022; Castro and Walter 2021; Wadsworth et al. 2020). There are two reasons we do not include eruptions of this type in our compilation. First is because our focus here is solely on explosive eruption processes, yet, it appears that hybrid events are characterized predominantly by processes at the effusive end of the spectrum, with explosive dynamics constituting only a small proportion of the eruption itself. For example, the 1730–1736 CE Timanfaya eruption of Lanzarote emplaced 0.51–0.54 km³ of material throughout the course of the eruption. Only 14(± 2)% of this volume is fall deposit and there is no evidence for ignimbrite emplacement within the eruption stratigraphy, suggesting that the remaining 86% likely comprises material ejected by basaltic lava effusion, and cone-building Strombolian activity (Muller and Longpré 2024). Fall deposits also account for 12(± 5)% of the estimated total bulk volume of the 2022 CE Tajogaite eruption of La Palma (Bonadonna et al., 2022), ~ 5% of the 1783–1785 CE Lakí eruption of Grímsvötn in Iceland (Thordarson & Self 1993), and < 50% of total mass emplaced during the 2011–2012 CE eruption of Puyehue-Cordón Caulle (Delgado et al. 2019; Pistolesi et al. 2015), with field-based evidence suggesting the remaining material was not emplaced solely by pyroclastic density currents. The second reason is because hybrid events are associated with highly complex eruptive dynamics; such that eruptive behaviour is (typically) associated with different magma ascent and mass eruption rates can occur simultaneously, with highly variable physical manifestations at the Earth’s surface (Wadsworth et al. 2022). Thus, their magnitude is more difficult to classify by traditional means (e.g., the VEI scale; (Newhall and Self 1982)).

Our database was constructed in two stages. Stage (1) involved compilation of a list of eruptions corresponding to the chosen magnitude range, and hence with the potential for inclusion. Our first source of this information was the Large Magnitude Explosive Eruption database (LaMEVE), which contains data for 1,883 explosive eruptions between M ≥ 4 and M ≤ 8.8, produced by 471 discrete volcanoes during the Quaternary (2.58 Ma to present; Crosweller et al. 2012). There are almost certainly a large number of eruptions missing from LaMEVE, given that a lack of geological investigations, failure to publish data, or difficulties accessing data and/or materials can all hinder the identification, characterization, and thus inclusion of discrete eruption events (Brown et al. 2014). Nonetheless, our objective here was not to acquire the most ‘complete’ eruption record, and was instead to explore as great a range of eruptions as possible for which volumetric data was available. Thus, LaMEVE provided a suitable starting point from which we could identify eruptions for inclusion. Following this initial compilation, eruptions were removed if they did not have (1) a robust estimate for total bulk deposit volume $V$ (in dense rock equivalent; discussed later), and (2) robust estimates of fall deposit volumes ${V}_{FALL}$ and/or ignimbrite volumes ${V}_{FLOW}$.

To ensure our database suitably accounted for eruptions at the upper end of the magnitude range, we also included data for M ≥ 8 events. However, only four events of this magnitude occurred during the Quaternary period, meaning most that are known fall outside the scope of the LaMEVE database. To acquire data for these events, we used the Collapse Caldera Database (CCDB), and the compilation of large-volume silicic explosive eruptions presented in Mason et al. (2004). The CCDB database contains information for ~ 473 caldera-forming eruptions that occurred prior to the Quaternary (> 2.58 Ma) including dimensions, morphology, age, magma composition, and geodynamic setting, and the compilation presented by Mason et al. (2004) lists information for 47 deposits formed by known eruptions of M ≥ 8. Of the 47 deposits listed in Mason et al. (2004), 40 have associated calderas and so are also listed in the CCDB.

Stage (2) involved quality-checking the remaining LaMEVE and CCDB eruptions through literature review. As part of the review process, the quality of volumetric estimates provided by the relevant literature was assessed based on whether these studies provided a thorough description of volume calculation methods, alongside clear consideration of error and uncertainties in volume calculations. Eruptions for whom the associated literature did not satisfy these criteria were subsequently removed. To independently verify the quality of the volume estimates provided in the reviewed literature, we conducted a thorough intra-study comparison by cross-referencing with instrumental observations (e.g., satellite measurements), and data acquired by a range of techniques (e.g., computer modelling, field cartography, petrography, stratigraphy). Dataset quality was also assessed by cross-comparison of geochronological and rheological information given for each eruption; given their importance for assessment of eruption dynamics and accurate calculation of eruption mass partitioning, respectively.

Following stage (2), a total of 1546 (95%) eruptions from the LaMEVE (2580–0 ka) dataset, and 354 (98%) from the CCDB dataset (> 2580 ka), had been removed. From this, 74 eruptions of magnitude $4\le M\le 9.1$ and spanning a diverse range of eruption types, locations, and magma compositions were deemed suitable for analysis. For each database entry that remained following this validation process, sources of the listed volumetric information are listed and fully referenced.

Magnitude (M) is used as the primary quantitative measure of eruption size in the dataset compiled herein. Our database does not include values for Volcanic Explosivity Index (VEI) or eruption intensity. We acknowledge that both VEI and intensity can provide effective scales with which eruptions of different sizes and expressions can be compared. However, both become less effective when considering the very largest eruptions (M > 7), given that there are no modern observations of comparably substantial events. For example, the VEI calculation requires both plume height and bulk volume estimates, which become subject to increasingly large uncertainties with increasing eruption age and/or size (Houghton et al. 2013; Mason et al. 2004; Pyle 2015). The ordinality of the VEI scale also means it cannot simultaneously reflect intra-eruptive variation in magnitude and intensity, complicating its application to continuously erupting or long-term eruptions (Neal et al. 2019), and/or unusually intense eruptions (Wilson 2001).

Data consolidation

Volume estimates for each eruption featured in our database are provided in both bulk cubic kilometres and as the dense-rock-equivalent (DRE) volume values. DRE value is the volume of magma or erupted material if it were non-porous and is therefore more reliable for comparison across systems. DRE volume of a given deposit ${V}_{i}$ (where $i$ refers to either fall ( ${V}_{FALL}$) or flow (${V}_{FLOW}$) for the fall deposit or ignimbrite case, respectively) is computed from the bulk porous volume ${V}_{b,i}$ by

$${V}_{i}={V}_{b,i}\frac{{\rho }_{T}}{{\rho }_{g}}$$

(1)

where ${\rho }_{T}$ is the tephra density, ${\rho }_{g}$ is the dense rock matrix density. For eruptions whereby ${\rho }_{T}$ is not available, we use ${\rho }_{T}\approx 1000 \text{kg }{\text{m}}^{-3}$. Similarly, where ${\rho }_{g}$ is not given explicitly, we use ${\rho }_{g}$ of $2200$, $2500,$ or $2700$ $\text{kg}.{\text{m}}^{-3}$ for rhyolitic, andesitic, and basaltic cases, respectively (Table S1). The bulk composition of the eruptive products from a given eruption is taken from the Global Volcanism Programme (GVP) database (v. 5.2.0; Global Volcanism Program, 2024), and instances in which we have applied this conversion are indicated the dataset. Our database can be used to define the ratio of total DRE volume to fall deposit DRE volume $V_{FALL}/V$.

With each eruption entry, we provide values for distance to the nearest ocean $D$, and distance to nearest ocean if travelling in the direction of the primary dispersal axis elucidated by field observations ${D}_{PDA}$ (Fig. 1). $D$ and ${D}_{PDA}$ values were obtained using the ‘ruler’ function in GoogleEarth™. This is discussed later.

Fig. 1

Schematic diagram illustrating how distance-to-ocean $D$ and distance-to-ocean along the primary dispersal axis ${D}_{PDA}$ were obtained for each eruption in our database. For caldera systems with no clear peak, eruptive centres were determined by taking the ejecta point source defined in the associated literature. Note that the geography depicted here is for illustrative purposes only and these two cases do not relate to a specific eruption

Uncertainty

Volumetric estimates for eruption deposits are usually subject to large uncertainties. Several factors can contribute to this, including the landscape in which materials are preserved, intense erosion of deposits following emplacement, logistical constraints on sample access and acquisition, and discrepancies in data collection and analysis methods (Engwell et al. 2015). All are effectively amplified when considering older eruptions (as we do here) and are variably documented or reported in the literature. Where uncertainty estimates on total eruption, ignimbrite, and/or fall deposit volumes were provided in the source literature, these were automatically added to the dataset, and listed alongside all relevant parameters. However, from the 74 eruptions included in our dataset, this information was provided for only 29 (39%). Hence, for the remaining 45 eruptions where uncertainty estimates on total eruption, ignimbrite, and/or fall deposit volumes were not provided, we follow Giordano & Cas (2021) and apply a conservative constant uncertainty of 45% of the total volume estimated.

Dataset overview

Our database compiles volume estimates and/or reconstructions for 74 eruptions of magnitude $4\le M\le 9.1$, spanning a diverse range of eruption types, locations, geological settings, and magma compositions (Table 1, Fig. 2). Approximately 20% of these eruptions took place in the Southern Hemisphere (SH), with Chile and Indonesia constituting the best-represented SH regions in the dataset. From the 74 subaerial eruptions listed in the dataset, 79% are situated in the Northern Hemisphere (NH). The NH contains the majority of the global landmass, meaning that eruptions taking place in this domain are likely to deposit large volumes of material on land to render them more easily accessible, and their source eruptions more likely to be documented (Buckland et al. 2020). The NH also possesses longer historical records, greater availability of resources for geological investigations, and is home to economies from which scientific funding is much more readily available (Rakotonarivo and Andriamihaja 2023). In the NH, Italy is best represented by number of eruptions. However, ~ 60% of these eruptions are produced by a single volcano (Mount Vesuvius), meaning that the USA leads if considering the number of eruptions produced by discrete volcanic systems (n = 7).

Table 1 Names, units, and definitions of data fields included in the database

Fig. 2

(a) Physical locations of the volcanoes attributed to the 74 eruptions included in our compilation, with individual volcanoes marked by triangles. Also shown is the number of $(M\ge 4)$ eruptions in our compilation corresponding to each 10° latitude zone. (b) Relationship between the ratio of total fall deposit volume (${V}_{FALL}$) to total volume erupted $(V)$ and eruption magnitude ($M$). (c) Relationship between the ratio of total ignimbrite volume ${(V}_{FLOW})$ to total volume erupted $(V)$ and eruption magnitude $(M)$. Error bars are marked in grey, and calculated according to the uncertainty estimates given in the source literature or assigned where these values were unavailable

Over 90% of the eruptions in the database are < 2.58 Ma old, with ~ 10% of these less than 50 years old. Three eruptions have listed ages > 2.58 Ma, and all took place in North America: the Fish Canyon Tuff (La Garita; 27,800 ± 0.85 ka), the Peach Spring Tuff (Silver Creek; 18,780 ± 0.2 ka), and the Rattlesnake Tuff (High Lava Plains; 7050 ± 1 ka). The most recent eruption listed in our database is that produced by Hunga Tonga Hunga Ha’apai volcano (Tonga Islands; −0.07 ka).

Data overview

Our database contains information for eruptions spanning several orders of magnitude, and thus a broad range of volumetric data (Fig. 3). The highest-magnitude (and most voluminous) eruption listed in the database is the Fish Canyon Tuff (FCT) eruption of La Garita caldera (Colorado, USA), with a listed magnitude of 9.1 and DRE volume of ~ 4,500 (± 2,025) km³. The crystal-rich FCT is also the oldest unit currently listed in the dataset (27,800 ± 850.5 ka) (Lipman et al. 1997). The largest Quaternary-age eruption in the database is the M8.8 Youngest Toba Tuff (YTT) eruption of Toba caldera (Indonesia; 73.8 ± 0.32 ka), with a DRE volume of ~ 3,783 (± 1,702) km³ (Costa et al. 2014). The smallest eruption listed in the dataset is the 2010 CE eruption of Eyjafjallajökull volcano (Iceland; M4.0), which erupted ~ 0.31 km³ (DRE) material over the course of 39 days, primarily in the form of fine-ash (Gudmundsson et al. 2012).

Fig. 3

Eruptions listed within the database presented as a function of latitude. Eruptions are colour-coded by their proximity to the nearest ocean using distances (a) $D$ or (b) ${D}_{PDA}$. White triangles (grey outline) mark eruptions for which the primary dispersal axis is not known. For both panels (a) and (b), the eruptions exhibiting the largest and smallest values are labelled. Age is presented on a logarithmic axis, given that the majority of eruptions listed in the database are those that occurred since < 1 Ma

The only major continent for which our database doesn’t contain at least one eruption is Antarctica (Fig. 2a). Volcanism is known to have played a critical role in the construction of the continent since ~ 200 Ma (Smellie 2020), and explosive eruptions throughout its history have left distinct signatures in ice cores extracted from the Antarctic Ice Sheet (e.g., Cole-Dai et al. 2021; Dunbar et al., 2017; Narcisi et al. 2017) and marine sediment successions extracted from the surrounding Southern Ocean (e.g., Di Roberto et al. 2021). However, the primary deposits produced by these eruptions are often near-totally covered by thick ice and snow (Burton-Johnson et al. 2016), and harsh conditions mean that mass burial, erosion, and degradation of deposits following an eruption is highly likely (Smellie 2021). This problem worsens over time, with older eruption deposits more likely to be buried under thick layers of ice, snow, and tephra from successive eruptions, eroded by tectonic activity, reworked by changes in sedimentation, and/or geochemically altered (Geyer et al. 2023). These issues are all compounded by the logistical challenges presented by Antarctica’s location, which makes it difficult to access and study sites critical for accurate volumetric estimates to be made. Nonetheless, advances in analytical methods, geochronology, and marine tephrochronology all offer the promise of data that can directly improve characterization of these events (Di Roberto et al. 2021), and so permit their future inclusion in this database.

Filtering data by geographical location

The 74 eruptions hosted within our database show a broad range of ${V}_{FALL}/V$ and ${V}_{FLOW}/V$ (Fig. 2b, c). There is no statistically significant pattern in the relationship between those ratios ${V}_{FALL}/V$ and ${V}_{FLOW}/V$ and eruption magnitude (Fig. 2), and so we suggest that the data may require filtering in a way that addresses the well-known issue that many explosive eruption deposits are poorly preserved globally. The main reason for this is deposition in remote and/or aquatic environments. Cores extracted from oceans (e.g., Mahony et al. 2020), lakes (e.g., Leicher et al. 2021), and ice sheets (e.g., Cole-Dai et al. 2021) can provide exceptional volcanic records. However, deposit preservation may be compromised by aqueous processes and/or limited by logistical challenges: both reducing overall preservation quality in these settings (e.g., Cassidy et al. 2014). To address this, we first consider the relationship between these ratios and the total land area per latitudinal domain, which shows that there are a large number of eruptions in latitudes of high land cover. Hence, it may be possible to filter the dataset to only examine eruptions that are at high values of distances $D$ and ${D}_{PDA}$ from oceans, increasing the chance that we are observing the full volume of fall deposits and ignimbrites.

Considering distance to nearest ocean $D$, 20% of the eruptions listed in the database are within $D$ <10 km from the ocean. These eruptions correspond to a range of magnitudes (4.4 to 7.7) and almost all are located in narrow, rifted arc or back-arc settings (LaFemina 2015). Five volcanoes within the dataset are (near) entirely surrounded by water: Krakatau (Indonesia), Kikai (Japan), Santorini (Greece), Kos (Greece), and Hunga Tonga Hunga Ha’apai (Kingdom of Tonga). Values for $D$ and ${D}_{PDA}$ remain the same for all listed eruptions produced by these three volcanoes, with the bulk of erupted products deposited into the ocean. The most recent of these five eruptions was the 2022 CE Hunga Tonga Hunga Ha’apai event, for which detailed characterization of emplaced deposits has been made possible by highly sophisticated sea-floor mapping techniques, satellite measurements, and tephrostratigraphic analyses (e.g., Clare et al. 2023; Seabrook et al. 2023). Similarly extensive sea-floor mapping and marine seismology has been combined with deep drilling projects, and detailed sedimentological analyses, to reconstruct the sub-aerial extent, volume, and distribution of the older (pre- Common Era) volcanic deposits produced by eruptions of the South Aegean (Santorini and Kos) and Pacific (Kikai) arcs (Geshi et al. 2017; Karstens et al. 2023; Kutterolf et al. 2021; Shimizu et al. 2024). Krakatau volcano lacks similarly detailed investigations, and difficulties in ascertaining the true volume of the 1883 eruption can be ascribed to a lack of direct instrumental measurements and/or observations, coupled with large-scale erosion of both marine and on-land deposits (Madden-Nadeau et al. 2021; Yokoyama 2015). Non-ocean island volcanoes $D$<10 km from their nearest ocean include Vesuvius (6 eruptions), and São Miguel (1 eruption).

Fifteen eruptions were produced by volcanoes located > 200 km away from the ocean by direct distance $D$. This includes the largest eruption listed (Fish Canyon Tuff; 27,800 ± 850 ka; M 9.1), and over 80% of these eruptions are M > 6. The dominance of high-magnitude eruptions in this bracket may be related to their continental tectonic settings, with extension of active continental margins in broader extensional belts (e.g., Yellowstone caldera, Valles caldera), and/or intraplate to rifted continental environments (e.g., Corbetti caldera) known to be associated with production of particularly high-volume (> 1000 km³) silicic eruptions (Bryan et al. 2010). If distance by the primary ash dispersal direction ${D}_{PDA}$ is considered, 30 eruptions occurred > 200 km from the nearest ocean. Similar to the upper $D$ bracket the majority of these are M > 6 and produced by volcanoes located in continental tectonic settings. However, deposits from seven (23%) of these eruptions have been detected in the marine realm. All are M > 6.5, and all except one were produced by volcanoes located in Central America; the only eruption in this group not located in this region is the Kidnappers’ Tuff eruption of Mangakino volcano, New Zealand (1010 ± 10 ka; (Wilson et al. 2009). Both New Zealand and Central America are characterized by a relatively narrow landmass, meaning that secondary/radial dispersal trajectories could cause simultaneous ash deposition in the ocean, even if (a) primary ash dispersal occurred on land, and/or (b) these eruptions at the upper end of the magnitude spectrum. Three eruptions of Mount Vesuvius, Italy (Avellino, Pollena, Pomici di Base) also fall within this > 200 km ${D}_{PDA}$ group, predominantly due to the north-eastern primary ash dispersion recorded for these three eruptions – sending ash toward the interior of the Italian landmass rather than into the Gulf of Naples (Rolandi et al. 2007).

Eight fewer eruptions exhibit distances < 10 km if we filter these events by their ocean proximity, where distance is measured in the direction of the primary dispersal axis ${D}_{PDA}$ rather than direct distance $D$. The difference in number can be ascribed to the removal of two eruptions from Campi Flegrei (Campanian Ignimbrite, Neapolitan Yellow Tuff), and 6 from Vesuvius (Mercato, 1631CE, 79CE, Pomici di Base, Pollena, Avellino). Similar to the example of the 1883 CE Krakatau eruption, this filtering further underscores geographical discrepancies in marine-based volcanological research for specific regions, and the lack of deep drilling in waters directly corresponding to certain dispersal patterns. For example, both the Fogo A eruption of Sao Miguel volcano (4.95 ± 0.36 ka) and the 1963 CE eruption of Agung volcano are < 10 km from the coast ${D}{\prime}$, and land-based investigations have identified ash dispersal trajectories from these two eruptions that would suggest substantial material deposition in the ocean. Yet, ash produced by these eruptions is yet to be identified in the marine realm.

Three eruptions show a > 3000 km difference between $D$ and ${D}_{PDA}$ values, and all are located on the west coast of the USA. These are the Mazama eruption of Mount Mazama (7.63 ± 0.15 ka; Egan et al. 2015), the 1980CE eruption of Mount St Helens, and the Bishop Tuff eruption of Long Valley caldera (764 ± 0.6 ka; Andersen et al. 2017). Field-based investigations have shown all to have dispersed ash to the east of the source vent, and thus into the heart of the North American continent.

Ocean distance, magnitude, and mass partitioning

Determination of mass-partitioning between ignimbrite and fall deposit volume is complicated by the fact that the majority of volcanic eruptions at the surface will deposit some (if not all) of its material in ocean basins. Identifying and sampling underwater deposits requires specialist equipment and vessels and often involves expensive and complex approaches (Mahony et al. 2020): both issues can severely limit the number of studies available to researchers for certain volcanoes, and may be further compounded by data reporting being incomplete and/or non-transparent (Aubry et al. 2021). Marine processes can also cause irreversible alteration of volcanic deposits in the ocean in the days, years, and millennia following an eruption, and disproportionately affects distal fall deposits (Cassidy et al. 2014). Hence, degrees of uncertainty in volume estimates will be broadly larger for eruptions whose material deposition was predominantly over water (e.g., Ballard et al. 2023; Vidal et al., 2015).

Figure 4 shows that eruptions occurring inland exhibit the strongest volume partitioning relationship between magnitude and ${V}_{FALL}/V$ when both $D$ and ${D}_{PDA}$ are considered. We have no reason to expect a strictly linear relationship between ${V}_{FALL}/V$ and $M$, but we nevertheless use a linear regression to characterise the extent to which there is a relationship of any kind between these parameters. This is testing the hypothesis that the largest eruptions may preferentially produce ignimbrites with only limited buoyant plume generation. Therefore, the Pearson’s coefficient of determination ${R}^{2}$ is given for each regression in Fig. 4. It is clear that the correlation between these parameters is strongest when $D>200 \text{km}$ is considered and that when the dataset is not filtered by $D$ or ${D}_{PDA}$, the correlation becomes substantially weaker.

Fig. 4

Relationship between the ratio of total erupted ash volume (V_FALL) to total volume erupted (V) and eruption magnitude for the 74 eruptions considered, filtered by their distance to nearest ocean ($D$; a-d). Plots shaded by grey (${D}_{PDA}$; e–h) filter eruptions based on their distance to nearest ocean if travelling in the direction of the primary ash dispersal axis listed in the literature. The significance of relationships is determined by Pearson’s correlation analyses, with R² values displayed here to indicate the ‘goodness’ of fit for the linear model applied to the two variables (full results in Table S3). Error bars are marked in grey, and calculated according to the uncertainty estimates given in the source literature or assigned where these values were unavailable (see methods)

Two effects could explain the relationship between $M$ and ${V}_{FALL}/V$ observed in Fig. 4, once the effects of both $D$ and ${D}_{PDA}$ have been considered. First is the tendency for intra-continental eruptions to be larger, due to the unique conditions that can result in the accumulation of more magma in the sub-surface reservoir. For example, upwelling of magma beneath the continental crust can deliver a consistent and long-term supply to shallow (eruptible) reservoirs over extended periods, and interactions between magma and the continental crust can lead to the development of highly evolved, volatile-enriched silicic magmas (e.g., Bryan et al. 2010; Van Zalinge et al. 2022). The sheer size of the ignimbrites emplaced during these events (typically > 1000 km³; Giordano & Cas, 2021) mean they are more likely to be preserved, and the intra-continental setting also offers more stable tectonic conditions compared to plate boundaries, where frequent tectonic activity can disturb or erode deposits. For example, the ocean surrounding Quaternary-active calderas in the Sumatran subduction zone has limited characterization of their fall deposits and ignimbrites (e.g., Bouvet de Maisonneuve & Bergal-Kuvikas., 2020; Forni et al. 2024; Kutterolf et al. 2023); especially when compared to the volcanic fields of the Western USA (e.g., Roche et al. 2022, 2016).

The second effect could be a greater propensity for ignimbrite formation over fall deposits as the magnitude of an eruption increases. This reconciles the fact that tectonic location alone cannot fully account for the shift in deposit type with increasing eruption magnitude, and would suggest that eruption dynamics also play a critical role in determining deposit proportions. This propensity could be explained by processes that occur exclusively in large vent structures. An example of such a process could be “eruption gargling” (Wilson and Hildreth 1997). Computer simulations show that interactions between the erupting (juvenile) mixture of gas and volcanic particles and dense pre-existing (fill) debris would cause uplift of a dense sheath of fill material along the margins of the erupting jet, causing pulsing behaviour that would prompt rapid collapse of the buoyant column, and production of localised pyroclastic currents (Valentine and Cole 2021). Ignimbrites produced by eruptions of Novarupta volcano (1912 CE; Hildreth and Fierstein 2012), Okmok volcano (~ 2 ka; Burgisser et al. 2024), and Mount Vesuvius (79 CE; (Scarpati et al. 2025)) provide further evidence for this effect and, given all were M > 6, could imply that it is one whose significance scales with magnitude (Valentine and Cole 2021). Therefore, size-mediated changes in source processes could reduce the proportion of fall material relative to the total erupted volume, and hence explain the decrease in ${V}_{FALL}/V$ with increasing magnitude (Fig. 4).

Given that none of the > M7 eruptions in our dataset were studied in real-time, we cannot rule out two alternative reasons for the absence of fall deposits associated with several of these events. First are effects produced by the complex, multi-faceted nature of debris entrainment. Computer model simulations have revealed the true scale of variability that can arise during eruption of polydisperse debris mixtures of varying ages, densities, temperatures, and speeds (e.g., Koyaguchi et al. 2018; Valentine and Sweeney 2018). For example, more stochastic jet behaviour can occur if entrainment is dominated by fine debris particles (Valentine 2023), which could result in the erosion, burial, and/or entrainment of any coeval fall deposits that are deposited within the vicinity of PDC runout; rendering them unidentifiable in the geological record. These deposits may also not be entrained instantaneously, and instead they may be entrained by successive eruptions characterized by (a) larger PDCs, (b) greater eruption intensity, and/or (c) a unique entrainment regime (e.g., Dowey & Williams, 2022). Phreatomagmatism poses an additional layer of complexity, where evidence has shown the presence of groundwater in an erupting vent, or introduction of water at various stages of the eruption, can influence PDC dynamics through generation of steam, increases in flow and plume density, and hydrothermal deposit alteration (e.g., dome collapse) (e.g., Wilson 2001).

Another reason for the absence of fall deposits associated with the largest eruptions in our dataset could be the risk of post-depositional degradation in distal locations, which increases as a function of time. Both the fall and PDCs produced by large eruptions can have an extremely wide spatial coverage, depositing material in both terrestrial and marine environments (e.g., Cisneros de León et al. 2021; Vidal et al. 2016). Even if predominantly over land, deposits can be altered by fluvial, aeolian, and gravitational processes ranging from small scale bioturbation, local topographic alteration, large scale debris flows, and resuspension by surface winds (Buckland et al. 2020). The effectiveness of preservation can also vary relative to the physical properties of the deposit itself, such as its composition, grain size, and density (Cutler et al. 2018; Dugmore et al. 2018), and ultimately reflect competing influences of plume behaviour, particle characteristics, atmospheric dynamics, and can all bias how deposits are retained in the geological record (Major 2023). Therefore, for the largest (and often oldest) eruptions in our dataset, we cannot rule out the possibility that fine ash particles were carried far beyond the maximum extent of the ignimbrites, but ultimately lost to post-depositional alteration, reworking, and degradation.

Discussion

Analysis of the database shows that eruptions of low relative magnitude are more likely to be dominated by buoyant plumes capable of producing both proximal and distal fall deposits, whereas larger relatively high magnitude eruptions (e.g., caldera-forming events) are dominated by ignimbrites (Fig. 4). It could be true that taller vertical plumes (high V_FALL/V, and therefore lower magnitude) may inject volatile magmatic gases higher into the stratosphere, and so exert a more significant (and longer-lasting) radiative effect (Timmreck 2012). Together, this would imply that a combination of physical, thermodynamic, and magmatic processes may all determine the extent to which an explosive eruption influences the surface environment, and the broader climate system, and so insights gleaned from our compilation already underscore a range of important, multi-disciplinary implications. Here, we outline a selection of these implications, and discuss their relevance to modern scientific research.

Explosive volcanic eruptions emit prodigious quantities of gas and solid material into the atmosphere. It is tempting to assume that the dose of gas/ash to the atmosphere (and particularly the stratosphere) scales with magnitude – larger eruptions mean more gas and ash. However, Fig. 5 suggests this may not be the case. Rather, that extremely large (M > 7) eruptions may interact with the atmosphere in fundamentally different ways. This is ultimately down to the ‘carrying capacity’ of buoyant plumes. It is clear that for moderately sized explosive eruptions, the entrainment and heating of ambient air is sufficient to impart buoyancy to an ash-gas mixture, lofting it to the atmosphere efficiently (Bursik and Woods 1996; Engwell et al. 2016; Woods and Wohletz 1991). However, as eruption magnitude is increased, there is a limit above which fountaining of material into pyroclastic density currents is favoured over buoyant plume production (implied by Fig. 4): an effect that is captured in both plume models (e.g., Baines and Sparks 2005; Jessop et al., 2016; Costa et al. 2018), and the geological record (e.g., Cas et al., 2011; De Silva et al. 2006). While such fountaining does not preclude fall deposits because large pyroclastic density currents produce co-ignimbrite plumes, such secondary plume formation will not be associated with stratospheric injection of volcanic gas (Marti et al. 2016). For these reasons, ultimately the fate of magmatic gases in those larger eruptions where ash fountaining is favoured is less well known and it is likely that the gas follows the ash flow to some extent, but it could equally be true that the gas is liberated and lofted independently (e.g., in vog plumes). It is at least clear that if fountaining dominates, then it is far less likely that volatile magmatic gases will reach the stratosphere.

Fig. 5

Conceptual diagram illustrating key differences between (a) high V_FALL/V ratio eruptions, and (b) low V_FALL/V ratio eruptions. For (a), the erupted mass is concentrated into a vertically-ascending, fine-ash enriched plume. This plume facilitates the injection of both fine-ash and volatile gases (e.g., SO₂) into the upper troposphere/stratosphere; the latter driving notable atmospheric disturbance through formation of aerosols, and ozone destruction. The mass erupted during (b) is also enriched in fine-ash particles but concentrated primarily in the thrust region, where materials are partitioned into local pyroclastic density currents via pyroclastic fountaining and plume collapse. Not only would this mitigate the net flux of SO₂ into the stratosphere, but it would also favor exceptionally high (~ 700 °C) temperatures, enabling more effective ‘scrubbing’ of SO₂ prior to atmospheric injection. Arrow tips indicate the direction of transport, and colour denotes the emission type. ‘Volatile gases’ (e.g., SO₂, CO₂, methane) are released from magma during explosive eruptions, with effects on atmospheric chemistry that are known to evoke significant climate disturbance (von Glasow et al. 2009). ‘Tephra’ is used here as an all-encompassing term for the solid products produced during an explosive eruption, including all grain sizes, compositions, and emplacement processes (Lowe 2011). Vertical dashed lines mark proximal, medial, and distal zones of deposition. Note: this figure is schematic and so not drawn to true scale. However, the stratosphere and troposphere are labelled as points of reference, and equate to altitudes of ∼10–50 km, and ∼0–10 km, respectively. No co-ignimbrite plumes are pictured

Explosive eruptions wax and wane in their eruptive intensity (e.g., Williams et al. 2014). During either waxing or waning, pyroclastic density currents can be produced, and will form forward-stepping or back-stepping onlap deposit structures in ignimbrites (Branney and Kolekaar 2002). Another feature of very large magnitude eruptions that is unknown is the extent to which they are pseudo-steady (inferred from the production of massive beds in ignimbrites; Branney and Kolekaar 1997) or time-varying. Where they are time-varying (e.g. see radial ‘unzipping’ of caldera walls; Williams et al. 2014), it is unclear if the waxing or waning phases of a given eruptive pulse exhibit a transition from no-plume to plume onset; more geological evidence is needed. However, it is useful to acknowledge here that waxing and waning during exceptionally large magnitude eruptions could result in shifts in plume dynamics, and thereby acknowledging that large magnitude eruption dynamics could be highly unsteady and time-varying.

An important application of the dataset produced here for (paleo)climatology would be to explore how the partitioning of erupted mass could influence the source-to-atmosphere pathway of SO₂, and hence the ensuing climate response. There is already a known mismatch between SO₂ budgets calculated using erupted volcanic products and glass inclusions, and SO₂ budgets directly recorded by satellites, which suggests there is ‘excess' sulphur (Carn et al. 2017; Shinohara 2008). Mass partitioning may further complicate the relationship between the SO₂ erupted at the vent, and the SO₂ dose provided to the atmosphere where aerosols can form.

There remain sizable gaps in our understanding of where plume dynamics fit into our present understanding of how volcanoes interact with the climate system. Our data suggest that eruptions exhibiting a high V_FALL/V ratio may inject volatile gases into the atmosphere more effectively (Fig. 5), and so be expected to produce more severe climate perturbations. If the size-mediated differences in eruption mass partitioning highlighted in our study do influence the dominant sink for erupted gases (e.g., SO₂), this would also influence the availability of SO₂ for climatologically significant processes, such as formation of sulphate aerosols (e.g., Diallo et al., 2017; Mather et al. 2003; Timmreck et al. 2010), and the destruction of ozone (e.g., Brenna et al. 2019; Klobas et al. 2017). For example, if the most voluminous eruptions do partition their erupted materials into relatively local pyroclastic density currents via pyroclastic fountaining and plume collapse (Fig. 6), this would limit the net flux of SO₂ into the stratosphere, and hence limit their capacity for climate perturbation. Supporting this proposition is a lack of large climate perturbations corresponding to, or following, many of the largest known volcanic eruptions of the past ~ 100-kyr (Baldini et al. 2015; Lohmann and Svensson 2022). Our dataset, therefore, provides a compelling case for further testing of the extent to which eruption-specific factors (e.g., plume morphology) may influence the climate system response to explosive volcanism relative to, and in conjunction with, changes in latitude (e.g., Fuglestvedt et al., 2024; Toohey et al., 2019), and magma composition (e.g., Kutterolf et al. 2023; Metzner et al., 2014).

Volcanic eruptions are events produced by a single, localised source, but also carry local, regional, and global consequences (Ayris and Delmelle 2012; Jones 2015; Marshall et al. 2022). Globalisation has added another layer of complexity to how we quantify the impacts of explosive volcanism, and with that a requirement for new ways to effectively manage their associated risks. However, the global-scale risks posed by explosive volcanism have historically focussed (almost entirely) on so-called ‘super-eruptions’ (M > 7), and largely assumed these events are just scaled-up versions of smaller, more recent events (Marzocchi and Papale 2019; Newhall et al. 2018). However, our study demonstrates that that high-magnitude eruptions may not produce buoyant plumes in the same way smaller eruptions do, and so the impacts of large (M ≥ 7) eruptions cannot be predicted by simply ‘scaling up’ the impacts observed following smaller, more recent eruptions. A combination of geographical, thermodynamic, and magmatic processes will all determine the type of phenomena produced by an explosive eruption, although the type and severity of the hazards they pose reflect interactions between these physical factors, and a vast combination of socio- and geo-political factors (Barclay et al. 2019). Finding new ways to explore and interrogate the geological record will undoubtedly yield a better understanding of these processes, and subsequently aid formulation of robust hazard assessments, early warnings, short-term forecasts, awareness campaigns, and communication strategies across a broad spectrum of potential eruption scenarios (Fearnley et al. 2018).

Conclusions

Are the largest explosive eruptions on Earth are fundamentally different from the smaller, more moderate events witnessed in recent history? To investigate this question in terms of eruption dynamics, particularly the partitioning of erupted mass, we compiled data on ignimbrite volumes, fall deposit volumes, eruption magnitude, and geographical location for 74 eruptions ranging from M ≥ 4 to M ≤ 9.1. Our analysis reveals a statistically significant, negative relationship between eruption magnitude and ${V}_{FALL}/V$, indicating that as eruption magnitude increases, the proportion of fall deposits decreases relative to the total erupted volume. This effect is most pronounced when eruptions are near the ocean, where the full volume of both fall deposits and ignimbrites is more likely to be preserved. Taken together, our findings suggest that smaller eruptions tend to produce more fall deposits and/or fewer pyroclastic currents, while larger eruptions are typically dominated by ignimbrites. Although Earth’s largest eruptions (M > 7) have not been studied directly, use of the geological record to gain a quantitative understanding of how their dynamics change in response to different physical, thermodynamic, and magmatic processes is an important, and valuable, next step toward reconstructing the environmental and climatic impacts of volcanism in the geological past. Moreover, the insights gleaned from these investigations may also assist in assessment and prediction of the hazards that future events may pose in vulnerable parts of the world.