Establishing the connection between large-scale tectonics and volcanic activity requires a conceptualisation of magmatic systems and the identification of potential factors that influence their temporal evolution. At convergent margins, the dehydration of subducting oceanic lithosphere induces partial melting of the mantle wedge, together with a contribution from decompression melting (Grove et al., 2009; England and Katz, 2010; Karlstrom et al., 2014). The thermal structure of the mantle wedge, together with the amount of water supplied by the dehydrating slab are important factors controlling where and how much melting is produced (Syracuse and Abers, 2006; Cagnioncle et al., 2007; van Keken et al., 2008; England and Katz, 2010; Karlstrom et al., 2014; Turner and Langmuir, 2015a, b). Modeling indicates that convergence rate and dip of the subducting slab exert a dominant control on the thermal structure of the mantle wedge (England et al., 2004; van Keken et al., 2008; Grove et al., 2009). Additionally, the age of the slab influences the amount of melting occurring within the mantle wedge. Older oceanic lithosphere is more likely to be affected by hydrothermal alteration, which results in larger amounts of H₂O being released during subduction leading to higher degree of partial melting of the mantle wedge (Poli and Schmidt, 2002). The effect of large-scale subduction architecture can be appreciated using the global geochemistry of magmatic rocks, in which incompatible trace elements correlate negatively with the thermal parameter of the slab (a product of the vertical descent rate and age of the subducting slab; Kirby et al., 1991), and positively with the thickness of the arc crust (Turner and Langmuir, 2015a, b). Additionally, a positive correlation has been observed between the arc-normal rate of plate convergence and the density of calderas, which was suggested as an expression of the association between high magma productivity in the mantle (assumed proportional to the convergence rate) and large volcanic eruptions (Hughes and Mahood, 2008).

Large volcanic eruptions are commonly fed by upper crustal reservoirs in which magma may reside for up to several hundreds to thousands of years (Cooper, 2019). The volumes and timescales over which eruptible magma is present depend on the thermal, chemical and physical evolution of magma reservoirs (Marsh, 1981; Jellinek and DePaolo, 2003; Huber et al., 2009; Caricchi and Blundy, 2015), which is intrinsically linked to the flux of magma intruded and rate of heat loss to the surrounding crustal rocks (Annen, 2009; Degruyter and Huber, 2014; Karakas et al., 2017). The ability to accumulate large volumes of eruptible magma in upper crustal reservoirs will depend upon the thermomechanical environment of the surrounding crust, and it’s ability to store rather than erupt magma (Jellinek and DePaolo, 2003; Karlstrom et al., 2017). The thermal maturation of the crust will promote magma storage (e.g., drop of crust viscosity; Jellinek and DePaolo, 2003) and the accumulation of the large volume of magma required to feed the largest eruptions (Magnitude ≥ 8; de Silva and Gosnold, 2007; Paterson and Ducea, 2015). For smaller eruptions, the timescales for eruptible magma accumulation are shorter (Cooper, 2019) and therefore the thermomechanical behavior of the upper crust does not change over time (de Silva, 2008). For these systems the accumulation of eruptible magma depends upon the depth(s) of magma storage and the geometry of magmatic bodies, which can be influenced by the crustal stress regime (Glazner, 1991; Vigneresse, 1995; Cruden, 1998; Menand et al., 2010; Gudmundsson, 2012; Chaussard and Amelung, 2014; Rivalta et al., 2015). The regional stress regime reflects the regional kinematics of subduction (Acocella and Funiciello, 2010), and alongside the thermal structure of the crust will influence the transport and storage of magma. Hence, the interaction between regional tectonics and volcanism extends beyond processes only associated with magma production in the mantle wedge, which is observed when analyzing the time-averaged mass eruption rate in different arcs (Acocella and Funiciello, 2010).

As summarized above, much research has been performed to understand how the dynamics of a subduction zone affect the generation, transport, storage and eruption of magma. However, the influence and importance of each of these processes on the size of eruptions in the volcanic record remains unclear. With the recent development of global databases summarizing both the volcanic eruption record and the kinematics of subduction zones, there is an opportunity to study the link between tectonic setting and the variability in the size of volcanic eruptions. Here we critically examine this link between tectonics and volcanism by using a wide suite of tectonic parameters and the Holocene eruption record of 162 arc volcanoes in convergent settings (Figure 1 and Supplementary Table A1). To perform our analysis we characterize the eruption record by quantifying the probability of recording eruptions of different magnitude. Because of the complex nature of magmatic plumbing systems and the multiple links between tectonics and volcanic activity, no clear relationships can be identified between distribution of eruption magnitudes and any single tectonic parameter (Acocella and Funiciello, 2010). Consequently, we employ a structured learning method that enables us to assess the strength of the probabilistic relationships between all the tectonic and volcanic parameters under consideration.

Our analysis focuses on convergent margin volcanism and utilizes two databases: a record of large-magnitude explosive eruptions (LaMEVE v.3.1; Brown et al., 2014), and a synthesis of subduction-related parameters (Lallemand et al., 2005; Laske et al., 2013; Figure 1). We make an assumption of exchangeability in which each volcano represents a possible state of the common “subduction zone magmatic system” and hence each volcanic record represents one possible sequence of eruption events that could occur at any volcano at any time during the evolution of a magmatic system (Sheldrake, 2014). Such a hierarchical structure lends itself to the Bayesian analysis that is fundamental to our methodology. Our exchangeability assumption is based on the fact that each volcano in our analysis is located in an arc possible of producing large-magnitude volcanic eruptions between Magnitudes 4 and 7. Consequently, implicit in our assumption of exchangeability is that the injection of magma does not alter the thermomechanical behavior of the upper crust and influence magma accumulation, which occurs for the largest magmatic reservoirs feeding super-eruptions (de Silva and Gregg, 2014). We refine the assumption by removing volcanic arcs undergoing significant extension (e.g., New Zealand; Stern, 1985). The assumption of exchangeability then allows us to investigate tectonic parameters that may influence variability in the behavior of each volcano, but not specifically control the size or style of eruptive activity of a single volcano.

As we are interested in the effect of regional-scale tectonic processes it is important to average the behavior of many volcanoes (e.g., Turner and Langmuir, 2015a). We are not attempting to explain why an individual volcano had a large or a small eruption. Rather, we are attempting to understand whether groups of volcanoes have larger or smaller eruptions on average, and if this different average behavior relates to their subduction setting. Grouping multiple volcanoes reduces the variability arising from local processes that may influence an individual volcanic record over a short time period. Selecting volcanoes based on their geographic location is a natural approach to grouping (e.g., Hughes and Mahood, 2008) as it can minimize within-group variance. For each volcanic arc we compare the mean of each tectonic parameter with the average eruptive activity recorded (Supplementary Table A2). The mean represents an average of all possible eruptive states for a subduction zone magmatic system.

In total, our analysis consists of 162 individual volcanic records constrained by the availability of tectonic parameters (Figure 1) and effects of under-recording, which is described in the following section. Due to the multi-dimensional character of subduction zones, in which parameters are likely related to each other, we first use a structure-learning algorithm to define a graphical model, which quantifies the probabilistic relationships between different variables, and guides our comparison of the eruptive behavior in different geographical regions.

In total, the eruption records of 197 volcanoes satisfy the under-recording criterion presented in Sheldrake and Caricchi (2017). Of these, 162 volcanoes have corresponding tectonic parameters that we analyze in this study. Variability in the magnitude of volcanic eruptions is calculated using the Holocene eruption record during the previous 11,700 years before 2000 C.E. This includes 278 Magnitude 4 events, 202 Magnitude 5 events, 63 Magnitude 6 events, and 15 Magnitude 7 events (Figure 3 and Supplementary Table A1). We also record the latitude and longitude of each volcano, which allows us to calculate the modeled crustal thickness at that location (Laske et al., 2013), and to attribute tectonic parameters from the closest corresponding arc segment (Lallemand et al., 2005; Table 1 and Figure 3). We do not consider the impact of climate on volcanism (e.g., Rawson et al., 2016; Sternai et al., 2017) or the history of volcanism at a given volcano.

The results of the structural learning algorithm are presented in Figure 6, in which the darkness of the edge is proportional to its probabilistic strength. In total only 22 possible networks were identified by the learning algorithm (Figure 6A). The strength of each edge was calculated by averaging its frequency across all 1000 networks (Figure 6B). Solid lines indicate edges have passed the data-driven significance threshold, which is 0.34. The variable that shows the greatest probabilistic relationship with eruption size is the obliquity of convergence. Other parameters that are related to the distribution of eruptions sizes (α) and are greater than the significance threshold are the thickness of the crust, the rate of subduction normal to the trench, and the age of the subducting slab. The subduction-parallel component is not significant.

The role of convergence obliquity is observed when we compare the eruptive behavior of different volcanic arcs. Volcanic arcs with the lowest values of alpha (α) are characterized by subduction that is normal to the arc (Figure 7). There is no clear linear relationship, however, between decreasing α and increasing convergence obliquity, which shows why it is complicated to identify relationships between single tectonic parameters and eruptive behavior of volcanoes located in different arcs. This is due to the roles of other parameters on eruption size, such as the age of the subducting slab and the rate at which the slab is subducting. Both of these parameters do not exhibit any clear trend with α, and volcanic arcs with the highest values for the rate of subduction or age of the slab do not produce proportionally larger eruptions on average (Figure 7). This corroborates the results of the graphical model that the obliquity of convergence has the greatest probabilistic strength with eruption size (Figure 6). Finally, we can see the influence of crustal thickness for some groups of volcanoes, such as those with moderate convergence obliquity (Figure 7A) or the oldest slab ages (Figure 7B), where the value of alpha decreases as crustal thickness increases.

To explore the combined effect of these tectonic parameters on volcanic activity in more detail we distinguish volcanic arcs into two groups based on the following parameter,

where A is the age of the slab and Qs = Sn / (Sn + Sp) is the subduction motion, which represents the motion of the plate boundary, and is generally reflected in the motion of the arc (Nakamura et al., 1977; Acocella and Funiciello, 2010). We utilize this combined parameter (H) to represent the average size of magma reservoir in the crust. It allows us to distinguish the effects of high mantle productivity and low convergence obliquity, which combined favor the development of larger crustal reservoirs. Mantle productivity increases with the age of the slab (A in Eq. 2; Poli and Schmidt, 2002). Neglecting the cases of arcs experiencing significant extension, increasing normality of convergence (Qs in Eq. 2) promotes the formation of deeper magma reservoirs (Chaussard and Amelung, 2012, 2014) by modifying the trajectory of dikes to favor the formation of sills (Menand et al., 2010; Menand, 2011; Daniels and Menand, 2015; Rivalta et al., 2015; Kavanagh, 2018). Therefore, with increasing H the likelihood of formation of larger reservoirs crust also increases. Using H, we identify two groups of arcs characterized by; (1) High mantle productivity and normal subduction (H > 3.0); and (2) Low mantle productivity combined with oblique convergence (H < 3.0). The value of 3 separating the two regimes was identified using a k-means clustering algorithm.

Volcanic arcs with higher values of H have moderate crustal thicknesses with an average of ∼30 km, except for Izu-Bonin with an average value of ∼19 km (Figure 8). The low-H group is characterized by the complete range of crustal thicknesses from ∼20 to ∼55 km (Figure 8).

We estimated graphical models for volcanoes in each of these two groups, which enables us to identify within-group distinct dominant behaviors (Figure 9). For the high-H group the rate of subduction has no significant correlation with the magnitude distribution of volcanic eruptions (Figure 9A). Instead, the age of the slab becomes the dominant parameter followed by the obliquity of convergence and then the thickness of the crust (Figure 9A). The smallest values of α are found in volcanic arcs that have moderate values for the convergence obliquity and age of the slab, whereas older ages and/or lower obliquity are associated with lower values of α (Figure 9B). For the low-H group the value of α decreases with increasing crustal thickness (Figure 9D). This is also visible in the results of the graphical model, which indicates that the strongest probabilistic link with eruption size is with crustal thickness (Figure 9C). Here the only exception is the Central Aleutians, which has exceptionally large values of S_p for this group (Figure 9D).

We have compared the tectonic setting to the eruptive records of individual volcanoes, using a Bayesian network. We have further compared the tectonic setting of different volcanic arcs to their average eruptive behavior. These two approaches have provided strikingly similar results and identified the obliquity of convergence, combined with the age of the slab and crustal thickness, as the most important parameters controlling the size distribution of volcanic eruptions (Figures 6–9). Using these results we discuss the possible processes that influence the size of volcanic eruptions up to magnitude 7 in different volcanic arcs.

When comparing different eruption records it can be difficult to isolate the influence of a single tectonic or crustal parameter. Subduction zones are characterized by multiple parameters and the bootstrapped approach to develop a probabilistic graphical model has allowed us to identify significant correlations between parameters. We have been able to validate this approach by comparing the results of a graphical model for 162 volcanoes, with the average behavior of 17 volcanic arcs. Our results shed light on how tectonic parameters influence the competition between magma accumulation and eruption in volcanic arcs. For large-magnitude explosive volcanic eruptions up to Magnitude 7, we identify two broad regimes: (1) high-H regimes, characterized by normal convergence that favors the accumulation of larger volumes of magma, which are thermally sustained by high mantle productivity associated with older slabs; and (2) low-H regimes, characterized by oblique convergence that favors transport and eruption of smaller volumes of magma. Within both regimes, we suggest that more frequent perturbation and eruption of upper crustal reservoirs will favor smaller eruptions on average. Within the high-H regime, sufficiently high mantle productivity is required to favor the accumulation of large volumes of eruptible magma. However, the highest mantle productivity will favor perturbation and eruption of magma and so the oldest slabs are not associated with the highest proportion of larger eruptions. Within the low-H regime mantle productivity is a less important control on the accumulation of magma and instead crustal transport associated with thinner crust will lead to more frequent and thus smaller eruptions. For the largest eruptions on Earth (e.g., Magnitude 8), however, we suggest the effects of slab age and crustal thickness on triggering eruptions becomes negligible, and instead the accumulation of eruptible magma is dominantly controlled by thermomechanical modification of the crust.

Within the field of Earth Sciences many complex highly multidimensional datasets exist, providing an opportunity to use graphical models and structured learning as an exploratory tool to disentangle probabilistic relationships within heterogeneous data sets. Where data are plentiful the directionality of edges in graphical models can be estimated, which can aid the distinction of correlation and causality between different variables (Scutari and Denis, 2014). For example, structural learning algorithms could be employed to understand what parameters control the dynamics of subduction, regional and local controls on the geochemistry of volcanic eruptions (e.g., Till et al., 2019) or the location (e.g., Andikagumi et al., 2020), geometry (e.g., Geyer and Marti, 2008), or morphology of volcanoes (e.g., Grosse et al., 2012). Such a structured approach to data analysis reduces the impact of pre-conceptual biases and provides a framework to develop data-driven hypotheses for physical relationships by incorporating additional data into the graphical models, to test hypotheses using numerical models, and ultimately help to shed light on geological processes. Where data are plentiful we envisage that our approach can be focused to understand more subtle variations in eruption sizes between volcanic regions, by incorporating parameters such as the history of volcanism or geometry of the mantle wedge. We hope our results and the methodology that we have used can stimulate further research into understanding the connection between tectonic kinematics and volcanism.

All datasets and analytical codes generated for this study are included in the Supplementary Material.

TS designed the study and performed the analysis. All authors helped with the manuscript preparation and discussion. MS helped with the statistical analysis.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.