Geochemical variability as an indicator for large volume eruptions in volcanic arcs

Caldera-forming eruptions have the potential to impact global climate and induce drastic socioeconomic change. However, the criteria to identify volcanoes capable of producing large magnitude eruptions in the future are not well constrained. Here we compile and analyse data, revealing that volcanoes which have produced catastrophic caldera-forming eruptions in the past, typically show larger ranges of long-term erupted bulk-rock geochemistry compared to those that have not. This observation suggests that geochemical variability is a measure of a magmatic systems size. Using a 2D thermal model that simulates the growth and evolution of crustal-scale magmatic systems by stochastic injection of dikes and sills, we show that such behaviour is consistent with differences in crustal magma fluxes. Higher injection rates accumulate greater melt volumes in more extensive crustal plumbing systems, leading to more variable distributions of temperatures and thus melt composition. We conclude that compositional variability should be included in the catalogue of criteria to identify volcanic systems with greater probabilities of producing future large eruptions. Importantly, this allows to identify stratovolcanoes with caldera-like geochemical signatures, which have not yet been recognized as systems with greater probabilities of producing large magnitude eruptions.


INTRODUCTION
Caldera-forming eruptions, classified with Volcanic Explosivity Index (VEI) of 7, occur on average once or twice every thousand years (Deligne et al. 2010;Sheldrake and Caricchi, 2017;Rougier et al. 2018). Analysis of past large magnitude eruptions provides evidence of considerable environmental impacts, including widespread cooling (Fischer et al., 2007;Sigl et al., 2015), which resulted in cultural disruption including forced migration and increased mortality (Nooren et al. 2017). In today's technologically interconnected world, an eruption of this type would likely have cascading impacts that go beyond historical accounts, stressing the need to address challenges in forecasting capabilities. At present, available geophysical imaging studies have not identified a particular volcano with large, interconnected melt volumes in the shallow crust, capable of producing a massive caldera-forming eruption. Although crustal magma bodies have accumulation histories that span hundreds of thousands to millions of years (Cooper, 2019;Frey et al., 2018;Weber et al., 2020a;Liu et al., 2021), diffusion studies indicate that the final melt mobilisation process may operate over timescales of years to centuries for large and smaller eruptions alike (Costa et al., 2020;Cooper et al., 2017;Druitt et al., 2012). This disparity between accumulation and mobilisation timescales may impede the long-term identification of candidate volcanoes by geophysical means. To anticipate volcanoes that could produce large eruptions in the future Newhall et al. (2018) developed and applied six criteria, including high rates of magma supply from the lower to upper crust. Such rates of magma transfer are inherently difficult to assess but recent progress in thermochemical modelling of magmatic systems indicates that the chemical diversity of volcanic rocks provides a record of crustal magma fluxes (Weber et al. 2020a). Here we further explore this hypothesis, starting Preprint from a data perspective and show that larger volcanic systems produce a greater variety of erupted bulk-rock compositions. Finally, we use a 2D thermal model to simulate the trans-crustal evolution of magma plumbing systems and show that such behaviour can be explained by differences in accumulation rates.

Data compilation
Accurately quantifying the bulk-rock chemical spread of volcanic systems requires that: 1) the long-term eruptive history is sufficiently well resolved; and 2) biases are identified, evaluated and, if problematic, eliminated. Given the large range of sampling strategies in geochemical studies and differences in the extent to which eruptive histories of volcanoes can and have been reconstructed, we argue that no single set of filtering criteria for existing large geochemical databases should be utilized to quantify the spread in bulk-rock geochemistry of volcanoes. Therefore, we carefully compiled data, using the GEOROC database as a starting point, which was complemented by literature search. Volcanoes were only included in the compilation with at least one published study, specifically targeting the long-term geochemical evolution of the system. If that condition applied, all available references presenting geochemical data for the volcano were examined to remove detailed studies on a particular eruption and oversampling of prominent eruptions, which could potentially bias the compositional spread towards lower values. We note, however, that both potential biases are rare based on available data. Further, studies with unclear spatial or temporal relation to the volcanic centre under investigation were removed. To be conservative in data inclusion, no filters based on analytical totals or loss on ignition were applied, as the risk of introducing biases for different rock groups this way Preprint outweighs slight differences in overall compositional spread resulting from analytical totals.
Clearly spurious data were, however, removed by only including analyses within a SiO2 contents between 48-78 wt.%. We focus only on arc volcanoes to avoid problems of intercomparison between different tectonic settings. The primary volcano type was assigned for each system based on the Smithsonian Holocene volcano list (https://volcano.si.edu/list_volcano_holocene.cfm) and calderas dataset of Hughes and Mahood (2008). For the sake of simplicity, we only distinguish stratovolcano(es), caldera(s) with diameter >10 km, and complex volcanoes, which includes larger dome complexes and composite systems for which a unique classification is not meaningful. We further assigned volcano lifespans based on geochronological data (Data Repository, Table S2) and the magnitude of the largest eruption taken from either the LAMEVE database (Crosweller et al., 2012) or a literature search for systems where data is available. The full data compilation used in this study is available in the Data Repository.

Thermal modelling
We use heat conduction theory to simulate the temporal evolution of temperatures in the Earth crust experiencing injection of magmatic sills and dikes. The heat conduction equation in two dimensions is given by: , (1) where T is the temperature, t is time, rho is the density of 2700 kg m -3 , c is the specific heat of 1 kJ kg -1 K -1 , k is the thermal conductivity of 2.7 W m -1 K -1 and Lc is the latent heat of fusion of 350 kJ kg -1 . The equation was discretized on a numerical grid using a fully explicit second order finite-difference scheme. Timestep and spatial resolution were chosen by convergence test. In all Preprint simulations we choose an initial linear geotherm of 30°C/km. Zero flux boundary conditions were applied at all sides except the surface which was set to 0°C. To include latent heat release, the non-linear melt fraction (Mf) temperature relation of Nandedkhar et al. (2014) was implemented. Latent heat release was implemented using the heat capacity method and a nonlinear iteration loop. Magma batches were injected as 80 m thick vertical dikes, randomly varying in length between 0.5 and 2 km and horizontal sills of 2 to 3 km length. Injection sites were randomized within a subdomain of 4 to 16 km depth and 10 km horizontal extent of the total 20x20 km model. Injected magma batches were accommodated by advection of crustal rocks to the sides in case of dikes and downwards for sills.

LONG-TERM COMPOSITIONAL DIVERSITY OF ARC VOLCANOES
We identified 54 volcanoes, distributed across 11 arcs, with sufficiently well characterised histories to quantify their bulk-rock compositional spread (Fig. 1). Within this compilation, stratovolcanoes comprise the largest group with 69%, followed by calderas (20%) and complex volcanoes (11%). To quantify the bulk-rock diversity for each volcano, we calculate the 95th percentile of SiO2 contents in weight percent, hereafter termed 'range', which includes 95% of the data by cutting off 2.5 % of each side of each volcano's distribution to reduce sensitivity to outliers. We note that the usage of a different index of differentiation (e. g. K2O) or different measures of compositional variability, like standard deviation or inter quartile range, do not substantially alter the structure of the data (Figs. S1, S2; Data Repository).
Comparing the ranges of SiO2 content for each volcano to ranges calculated for data from the unfiltered GEOROC database, most volcanic systems fall close to a 1:1 line ( Fig. 2A). Notable differences between the compilation and the GEOROC database can be explained by the Preprint overrepresentation of the prominent eruptions, and data of unclear provenance that has been removed while compiling data. It should be noted that the number of analysed samples does not correlate with the range of erupted bulk-rock compositions ( Fig. 2A).
Erupted whole-rock SiO2 distributions show large differences in the range and median of erupted compositions and reveal systematic differences between volcano types (Fig. 2B). SiO2 ranges in our compilation vary between 4 wt.% for Uturuncu volcano (Bolivia) and 28 wt.% for the Los Humeros caldera in Mexico, with many systems recording intermediary ranges of ~15 wt. %. The most important observation from Fig. 2B is that caldera volcanoes clearly show a systematically higher spread in erupted bulk magma chemistry compared to stratovolcanoes.
Using a Brown-Forsythe test for the equality of group variances, we show that erupted bulk-rock SiO2 ranges differs significantly among stratovolcanoes (mean range: 13.3 wt.%) and calderas (mean range: 21.7 wt.%) with p-value of 0.6646. To assess the impact of the variable sample size, calculated SiO2 ranges were resampled 1000 times with replacement, confirming our result that caldera volcanoes show systematically higher ranges in erupted chemistry. (Fig. S3, Data repository). Interestingly, several stratovolcanoes like Three Sisters (Cascades) or Ceboruco (Mexico) have erupted compositional distributions, which are more similar to calderas rather than other stratovolcano or complex systems in the compilation.
Before discussing underlying mechanisms and implications of this finding, we turn to the evaluation of how the observed differences may be impacted by spatial sampling effects.
Calderas typically comprise a larger spatial area compared to most stratovolcanoes, which could potentially bias the results due to differences in the extent of sampling. It is therefore pertinent to compare the spread of bulk-rock geochemistry for stratovolcanoes and calderas based on samples within equal areas. We use the Mexican Volcanic Belt as a test case to visualize and compare the Preprint spatial distribution of erupted compositions, given that it contains a comparably large number of volcanoes in our compilation (Fig. 3). Compositional distributions for eight volcanoes were calculated by including all available data of the GEOROC database within equal 'buffer areas' of 30 km diameter surrounding the main vent. Clearly, the shape of the resulting distributions based on spatial sampling confirm our results obtained by analysis based on a volcanic system approach (Fig. 2, 3). All calderas in the Mexican Volcanic Belt show broader SiO2 distributions compared to four out of five stratovolcanoes. Again, Ceboruco volcano shows more caldera-like geochemical patterns, which may relate to fundamental differences in the architecture or dynamics of its magmatic systems compared to other Mexican stratovolcanoes.

CONTROLS ON VOLCANIC BULK-ROCK DIVERSITY
Caldera-forming eruptions require the accumulation of large magma volumes, suggesting that the conditions favouring this also result in wide ranges of erupted bulk-rock compositions ( Fig S4, Data Repository). Due to the coupling between melt temperature and chemistry (e.g., Sisson et al. 2005;Ulmer et al., 2018), we suggest that geochemical variability is related to thermal heterogeneity in spatially extensive magmatic plumbing systems (Weber et al., 2020a).
We suggest the primary factor favouring variability is higher magma injection rates, which is shown using 2D thermal modelling of stochastic magma injection and evolution in the crust (Fig.   4). Whilst magmatic systems may experience geochemical variability with time (e.g., Bacon and Lanphere, 2006;Hora et al., 2017;Singer et al., 2008), we suggest at the global scale volcano age is less significant than magma flux, as the lifespan of volcanoes with large geochemical variability varies considerably (Fig. 2b). Furthermore, there are many examples of long-lived volcanoes (e.g., Uturuncu, Aucanquilcha, Xianantecatl) that have produced monotonous Preprint andesitic-dacitic magma chemistry throughout their >1.5 Ma histories. Nevertheless, some of the longest-lived caldera volcanoes exhibit high geochemical variability, suggesting that the maturity of a magmatic system can still be an important factor to accumulate large volumes of eruptible magma.
The development of large caldera systems results from the complex interplay of tectonic and crustal processes that influence the ability to accumulate large magma volumes (Hughes and Mahood, 2008;Sheldrake et al., 2020). The regional impact of these processes is observed in the Trans Mexican Volcanic Belt. Here the active volcanic front is dominated by similar lowvariability magmas, whereas further from the active front volcanic centres produce a diverse array of magma geochemistry (Fig. 3). An identical pattern is observed in the type of volcanoes, with the active front characterised by andesitic stratovolcanoes in comparison to large calderas (Ferriz & Mahood, 1984). One exception is Ceboruco, which has not produced a large caldera, but lies behind the active front. Nevertheless, Ceboruco exhibits large geochemical variability (Fig. 3). We suggest this variability results from the same processes that favour the accumulation of magma feeding large caldera-forming eruptions at other systems and is possibly related to slab depth (Fig. 3). Such close spatial association of compositionally monotonous and highly variable volcanic systems emphasises the important role of local differences in crustal magma fluxes (Kent, 2013;Wörner et al. 2017;Till et al., 2019;Weber et al., 2020a).

IMPLICATIONS FOR VOLCANIC HAZARD MITIGATION
The relationship between magma diversity and igneous plumbing system size, as indicated by our compilation and modelling, can be used as an indicator of hazardous volcanoes with greater likelihood of producing caldera-forming eruptions. This is important in the pre-caldera phase (i.e., de Maissonneuve et al., 2021), especially for stratovolcanoes that can transition to largemagnitude caldera forming eruptions but have not experienced such events previously. Several stratovolcanoes in our compilation (e.g., Three Sisters, Ceboruco, St Helens) show distributions akin to those of calderas, indicative of larger plumbing systems size, but these systems have previously not been recognized as capable of producing VEI7 eruptions based on the criteria laid out in Newhall et al. (2018). We therefore suggest that the compositional spread of volcanoes is included in the list of criteria to evaluate the long-term risk stemming from different volcanoes worldwide.