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

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.

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.

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 Repository¹). 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). Well-studied 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).

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.

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.

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 (α = 2.27; Fig. 3A) than for stratovolcanoes (α = 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 power-law behavior that can be either steeper (i.e., proportionally more smaller-magnitude events) or shallower (i.e., proportionally more larger-magnitude 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.

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 under-recording (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 λ = 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.

We present a new method to approach under-recording, 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.

Many thanks to Mark Bebbington, Valerio Acocella, and an anonymous reviewer for their informative comments, which greatly improved the manuscript. Sheldrake and Caricchi are supported by ERC grant 677493-FEVER.