Regional variability in the frequency and magnitude of large explosive volcanic eruptions

Quantifying the frequency at which volcanic eruptions of different size occurs is important for hazard assessment. Volcanic records can be used to estimate the recurrence rate of large-magnitude eruptions (magnitude ≥4), but recording biases that impact data completeness complicate analysis. To overcome these biases, we conceptualize the volcanic record as a series of individual and unique time series associated by a common behavior. Thus, we approach issues of completeness on a volcano-by-volcano basis and use a hierarchical Bayesian approach to characterize the common frequency-magnitude (f-M) behavior for different groups of volcanoes. We identify variations in the f-M relationship between different volcano types and between different volcanic arcs. By accounting for systematic under-recording in the volcanic record, we also calculate the global recurrence rates for large-magnitude eruptions during the Holocene, which are similar to previous estimates. However, higher recurrence rates for smaller-magnitude events are observed, which is a result of our adjustments for data completeness. Quantifying how the f-M relationship varies between different groups of volcanoes provides an opportunity to understand how the tectonic setting influences f-M behavior, which is important to quantify long-term regional volcanic hazard.


INTRODUCTION
Various physical processes, occurring across a range of spatial and temporal scales, control the frequency and size of volcanic eruptions (Caricchi et al., 2014).Quantifying the frequency-magnitude (f-M) relationship for different groups of volcanoes, or for different regions on Earth, provides an opportunity to link these physical processes to the tempo and size of volcanism.Calculation of the f-M relationship is based on analyzing the eruption record (e.g., Mason et al., 2004;Deligne et al., 2010).The eruption record, however, suffers from issues of data completeness related to natural (aleatory) sources of bias associated with preservation of deposits (Deligne et al., 2010;Brown et al., 2014) and human-induced (epistemic) sources of bias resulting from the short duration of historical records (Coles and Sparks, 2006;Furlan, 2010;Rougier et al., 2016).This systematic variation in under-recording limits the ability to analyze records from different time periods or even different regions on Earth.We develop a method to account for under-recording biases, which permits the analysis and comparison of the eruption record for different groups of volcanoes and different volcanic regions on Earth.

DATA
The size of an eruption can be represented using the magnitude scale (Pyle, 2000), which is a continuous measure of mass: M = log 10 [erupted mass (kg)] -7. (1) We use the LaMEVE database (version 3.1), which provides a record of large-magnitude explosive eruptions (M ≥ 4) throughout the Quaternary (Crosweller et al., 2012).The data suffer from a bias associated with measurement error, meaning a disproportionally large number of events are rounded to the nearest integer (Brown et al., 2014).To overcome this bias we group eruptions into different magnitude classes (e.g., magnitude 4.0-4.9= M4), rather than try to fit a parametric distribution to the data.
In our analysis, we make an assumption of exchangeability (Bebbington, 2014;Sheldrake, 2014), which implies that fundamental processes leading to large-magnitude explosive eruptions, associated with the accumulation of eruptible magma and triggering of an eruption, are shared by different volcanoes.The frequency of eruptions of M < 4 is potentially influenced by additional processes such as second boiling (Tait et al., 1989), and hence smaller eruptions may not satisfy the exchangeability assumption.Similarly, we restrict the upper bound of the data set to M7, as the processes triggering larger eruptions could differ from those responsible for eruptions of lower magnitude (Deligne et al., 2010;Gregg et al., 2012;Caricchi et al., 2014).To further satisfy the assumption of exchangeability and common geo-tectonic processes, we restrict the analysis to arc volcanoes.

UNDER-RECORDING ANALYSIS
Traditionally, completeness in the eruption record is analyzed at a group level (e.g., globally or regionally).Change points in data completeness are identified based on analyzing the rate of volcanic events (Coles and Sparks, 2006;Furlan, 2010;Mead and Magill, 2014;Rougier et al., 2016), which allows return periods for different eruption magnitudes to be estimated.In the Bayesian approach we adopt, the f-M behavior is quantified based on the proportion of each eruption magnitude from individual, but exchangeable, volcanic records.Thus, we attempt to statistically account for under-recording based on the proportion of each eruption magnitude in a period of time, rather than the number of events for each magnitude.
We search for an observation window (t start t present , where t is time; Fig. 1) in which the processes that control the preservation of volcanic deposits are consistent through time, so that the proportion of each eruption magnitude is constant.Essential to our approach is the notion of unique records, where volcanoes are only included in the analysis if their oldest recorded eruption occurred before a particular date (t unique ; Figs. 1A and 1B).This is based on the assumption that for the whole tephrostratigraphy at one volcano, any deposit that has not been eroded and lies above (i.e., is younger than) the oldest observed eruption is recorded.For an observation window, as the value of t unique changes, the proportion of each eruption magnitude must remain constant.Importantly, changing the value of t unique only determines whether the record of each volcano is included in the analysis (Fig. 1B).If a volcano is included (i.e., it has at least one eruption older than t unique ), all data within the observation window are considered and are used to calculate the proportion of each eruption magnitude for each value of t unqiue (Fig. 1C).
Using a value of 2000 CE for t present , we perform a chi-square test of homogeneity for the proportion of eruption magnitudes for two values of t start : 50 ka and 11.7 ka (i.e., the Holocene) (see the GSA Data Repository1 ).The hypothesis of homogeneous under-recording is rejected for t start = 50 ka (p ≈ 0; Fig. 2A) but accepted for t start = 11.7 ka (p = 1; Fig. 2B).The lack of homogeneity in the longer observation window may be due to a combination of decreasing number of volcanoes satisfying the t unique criterion and preferential recording of larger eruptions (Kiyosugi et al., 2015).Alternatively, these changes may also represent a natural process related to changes in the rate of preservation associated with glacial erosion, or perhaps to changes in the dynamics of magmatic systems as a result of deglaciation (Watt et al., 2013;Rawson et al., 2016).
For the Holocene observation window, when the value of t unique is younger than ca.0.5 ka, the proportion of M4 events is not constant (Fig. 2B), which is due to a change point in the recording rate of M4 events (Fig. 2C).This suggests that for the Holocene, once volcanoes that have only a historical record (i.e., dominated by smaller M4 events) are removed from the analysis, the rate of under-recording is constant for all eruption magnitudes (Deligne et al., 2010).Wellstudied volcanoes with records older than 0.5 ka are not subject to a historical bias.Indeed, the majority of volcanoes meeting the 0.5 ka criterion are in regions where extensive research has been undertaken (e.g., Japan, n = 44), and very few are in regions where the geological record is scarce (e.g., Indonesia, n = 6).

FREQUENCY-MAGNITUDE RELATIONSHIP
To statistically characterize the f-M relationship for volcanic eruptions, we use the hierarchical Bayesian method set out by Sheldrake (2014) (see the Data Repository).In this approach, the proportion of eruptions of different magnitude is quantified considering groups of volcanoes associated by a common behavior (e.g., arc volcanoes), while maintaining that each volcano is characterized by a unique behavior.The statistical model is set up as uninformative, so that prior to observing the data, each eruption size is equally likely and there is no subjective judgement.Once updated by the data, the Bayesian model is used to calculate (1) a prior distribution that characterizes the common global f-M behavior, and (2) a series of posterior distributions that characterize the f-M behavior at individual volcanoes.The prior distribution is modeled using a Dirichlet distribution, as we characterize eruption magnitude as a discrete multivariate data set (i.e., mutually exclusive events) where the probability of the different events adds to unity.The Dirichlet distribution is advantageous, as it does not put any restrictions on the shape of the distribution, allowing different behaviors to be identified for different groups of volcanoes.
Based on the analysis of under-recording, we include only volcanoes that have a record older than 0.5 ka and count events that occurred only during the Holocene.In total, 197 volcanoes and 664 eruptions (our global data set) meet these requirements and are included in the statistical analysis.The posterior probabilities (Pr) for different magnitude events (m) occurring in different groups of volcanoes are characterized by powerlaw behaviors (Fig. 3; see the Data Repository): where g is a decay parameter and is the inverse of the power-law exponent, which enables the comparison of f-M behavior between different groups.Not all posterior distributions for individual volcanoes exhibit a power law (Fig. 3) because all states of the common magmatic processes are not necessarily observed at a single volcano, especially over the time scale of the Holocene.As an example, a volcano might be characterized by frequent magma replenishment leading to closely spaced smaller-magnitude eruptions, such that the inter-eruptive time is not sufficient to accumulate magma to feed larger-magnitude eruptions.A series of intervals in observation window is chosen (t unique = {t start , t 1 , t 2 }), and if individual volcanic record has no events older than t unique , then volcano is not included in analysis for that t unique (e.g., volcano A for all values of t unique ).C: For each magnitude, number of events from all volcanoes that satisfy t unique criterion are summed (Σ M ), and relative proportions (P M ) calculated for each value of t unique .Importantly, an aleatory level of under-recording still persists in the eruption record, as not all volcanic deposits are well preserved (e.g., 1991 CE eruption of Mount Pinatubo, Philippines; Gran et al., 2011).Consequently, the slopes of the calculated power laws are unlikely to represent the true f-M behavior for volcanic eruptions.Nevertheless, the level of aleatory under-recording is stationary in time (Fig. 2B), allowing different groups of volcanoes to be compared.
To investigate the robustness of the method, we distinguish volcanoes based on their morphology, as presented in the LaMEVE database.The two most populous volcano types in our analysis are calderas (n = 34) and stratovolcanoes (n = 121).The power-law exponent describing the common behavior in these two groups is lower for caldera volcanoes (a = 2.27; Fig. 3A) than for stratovolcanoes (a = 2.63; Fig. 3B), which is expected given that calderas are associated with larger volcanic eruptions (Brown et al., 2014;Cashman and Giordano, 2014;Whelley et al., 2015).This gives us confidence that the methodology correctly characterizes the volcanic record, and allows us to investigate different volcanic arcs (Figs.3C-3E).
Variation in the f-M relationship is observed between different volcanic regions, with powerlaw behavior that can be either steeper (i.e., proportionally more smaller-magnitude events) or shallower (i.e., proportionally more largermagnitude events) than the global f-M relationship (Fig. 3).This ability to compare different volcanic regions offers a unique opportunity to understand how tectonic setting influences the frequency of explosive eruptions between M4 and M7.As an example, here we look at the role of crustal thickness.At a global scale, there is no clear relationship between crustal thickness and the steepness of the power-law behavior (Figs.3C-3E), although our data restrain us from exploring this relationship fully as different parameters influence the rate and explosivity of arc volcanoes (e.g., Hughes and Mahood, 2008;Acocella and Funiciello, 2010).At a regional scale, however, when comparing different portions of the Japanese arc where magma output is similar (e.g., Izu-Bonin versus Honshu; Table 1), the steeper power-law behavior is observed for sub-regions with the thinnest crust, which also have the lowest caldera density (Table 1).Nevertheless, it is clear that at a global scale, other parameters such as convergence rate or obliquity of convergence (Table 1) will be more important in controlling the f-M relationship.
Quantifying the f-M relationship for Holocene eruptions is an important step to perform long-term volcanic hazard assessment at regional and global scale, but requires accounting for aleatory sources of under-recording.We calculate the global f-M relationship (see the Data Repository) by assuming that the level of underrecording (l) varies linearly with the power-law behavior, meaning that natural sources of underrecording (e.g., erosion) remove from the record a proportionally greater number of smaller events with respect to larger ones: (4)  Additionally, on the basis of previous studies, we assume that (1) the global recording rate of M4 events since 1961 CE is complete, based on a calculated change point in underrecording (Furlan, 2010), which allows us to compute the expected number of M4 events during the Holocene; and (2) the completeness in the global record of M7 events in the Holocene is 70% (Brown et al., 2014).Importantly, there is no time dependence in Equation 4, as our calculations are based on the proportion of events, rather than a count of events, and so we assume that any time dependence is removed in the under-recording analysis.
Fitting the global relationship from Figure 3F, we obtain a value of l = 0.319, which we use to estimate the number of eruptions per unit time.At a global scale for arc volcanoes (~88% of the LaMEVE database), the expected number of eruptions of each magnitude in the Holocene is similar in order of magnitude to that obtained from the analysis of a single global data set of 2000 yr (Mason et al., 2004;Fig. 4).However, the ability of our method to adjust for incompleteness is reflected in the estimates for M4, M5, and M6 eruptions that are, respectively, 50%, 33%, and 14% larger than those of Mason et al. (2004), supporting the consensus that smaller eruption magnitudes suffer greater levels of under-recording in the Holocene (Brown et al., 2014).These observations provide support to the validity of our method, suggesting the Bayesian approach can be used to identify variations in the f-M behavior for different groups of volcanoes, and potentially for other geological processes such as earthquakes.

CONCLUSION
We present a new method to approach underrecording, based on viewing the eruption record as a group of unique time series rather than a single global data set.With this approach, the rates of under-recording are systematically constrained.Using a hierarchical Bayesian approach, these unique records have been analyzed to understand common f-M behavior for groups of volcanoes.The results indicate that f-M behavior varies for different volcanic arcs, which is essential to developing appropriate volcanic hazard models for regions with a record of volcanic eruptions.Furthermore, quantifying variability in f-M behavior potentially provides insights into how tectonic setting controls volcanic activity.
Tackling the analysis of unique volcanoes to understand common processes is a major frontier in volcanology (Cashman and Biggs, 2014).Bayesian methods provide an opportunity to understand common physical and geological processes from a series of unique data sets, which here are individual eruption records.With increases in the size of global data sets, we envisage that the application of Bayesian methods to other multivariate data sets, such as geochemistry, may provide an opportunity to link magmatic processes to the size and tempo of volcanism.

Figure 1 .
Figure 1.Diagram presenting methodology to analyze group completeness for unique volcanic records (volcanoes A-D).A,B: For a specific observation window (t start -t present , where t is time), number of events of different magnitude (indicated by bars of different height and color) is counted.A series of intervals in observation window is chosen (t unique = {t start , t 1 , t 2 }), and if individual volcanic record has no events older than t unique , then volcano is not included in analysis for that t unique (e.g., volcano A for all values of t unique ).C: For each magnitude, number of events from all volcanoes that satisfy t unique criterion are summed (Σ M ), and relative proportions (P M ) calculated for each value of t unique .

Figure 2 .
Figure 2. Analysis of under-recording for magnitude M4-M7 eruptions (using the same color scheme as Fig. 1) from unique records at arc volcanoes for two observation windows (see Fig. 1 for explanation of observation window and parameters).A: For t start = 50 ka (t is time), value of relative proportion P M for a particular magnitude is not stationary (each diamond represents different value of t unique ).B: For t start = 11.7 ka, value of P M is approximately stationary and level of under-recording constant.C: Number of M4 events in our Holocene data for different values of t unique , with change point in recording rate calculated at 0.47 ka (see Data Repository [see footnote 1]).

Figure 3 .
Figure 3. Frequency-magnitude (f-M) behavior for different volcano morphologies (A,B), different volcanic arcs (C-E), and all volcanoes (F) in the analysis.Plotted on each chart are posterior results for individual volcanoes (dashed lines), fitted power-law behavior for respective group (crosses), number of volcanoes in each group (n), and value of power-law exponent (a; 2 standard deviations in parentheses).Inset on each chart is probability (prob.)density plot for crustal thickness at volcanic centers for each group using data from the global crustal model Crust 1 (Laske et al., 2013).