Triggering of volcanic degassing by large earthquakes

The potential for triggering of volcanic activity by large earthquakes has been pondered for centuries: Pliny the Elder addressed the topic in A.D. 77–79 (Holland, 1601). Volcanic responses to earthquakes are difficult to quantify in a robust manner, however, largely because only a small proportion of global volcanoes are adequately monitored. Despite this, there is mounting evidence for earthquake-triggered volcanic activity on various time scales (e.g., Linde and Sacks, 1998; Brodsky et al., 1998; Manga and Brodsky, 2006; Cigolini et al., 2007; Harris and Ripepe, 2007; Walter and Amelung, 2007; Walter et al., 2007; Watt et al., 2009; De la Cruz-Reyna et al., 2010; Delle Donne et al., 2010). Some of these studies have been catalyzed by decade-scale satellite remote sensing observations of volcanic activity (e.g., Delle Donne et al., 2010). Quantitative measurements of changes in both the earthquake-generated state of stress and volcanic activity (including ground deformation, heat flux, and degassing) may provide greater insight into the triggering mechanism as well as the types of volcanoes that may be susceptible to triggering. However, to date there has been no analysis of the volcanic response to earthquakes based on direct observations of volcanic degassing.

Here we assess the links between volcanic degassing and dynamic (not static) stresses propagated by large earthquakes by modeling surface wave propagation from M7 and greater earthquakes in 2004–2010 to a selection of target volcanoes with sulfur dioxide (SO₂) emissions measurable by satellite remote sensing. Surface wave modeling and the associated peak dynamic stress (PDS) calculations provide a robust method for identifying the earthquakes most likely to trigger a change in volcanic activity. SO₂ is the optimal indicator of volcanic degassing because it is easily distinguishable against the atmospheric background; for this reason it is the most commonly measured volcanic gas, and thus it is also crucial to understand processes that modulate its emission. We use daily measurements of volcanic SO₂ emissions from the Ozone Monitoring Instrument (OMI) on the NASA Aura satellite (Carn et al., 2013) that provide a consistent measure of volcanic degassing to compare with calculated PDS. The availability of daily satellite measurements of volcanic SO₂ emissions permits such an analysis on a large scale for the first time.

To quantify daily volcanic SO₂ emissions, we use the operational OMI SO₂ product derived using the linear fit algorithm (Level 2 OMSO2 data product; Yang et al., 2007), obtained from the NASA Goddard Earth Sciences Data and Information Services Center (DISC; http://disc.sci.gsfc.nasa.gov/). We selected the following target volcanoes, covering a variety of eruption styles in diverse arc settings, from those with SO₂ emissions detectable by OMI (e.g., Carn et al., 2017): Ambrym and Gaua (Vanuatu); Bagana, Rabaul, and Ulawun (Papua New Guinea); Fuego and Pacaya (Guatemala); Merapi, Semeru, and Bromo (Indonesia); Turrialba (Costa Rica); and Villarrica (Chile) (Table DR1 in the GSA Data Repository¹). These volcanoes were all in an open-vent state in 2004–2010 (i.e., with a magma-atmosphere interface at or close to the surface), and therefore may be more sensitive to small earthquake-induced pressure changes. The eruptive styles represented range from active lava lakes (e.g., Ambrym, Villarrica) to open-vent degassing (e.g., Turrialba, Ulawun) to lava dome growth (e.g., Merapi, Semeru, Bagana).

For each target, we created a time series of daily SO₂ loading measured by OMI between 1 October 2004 and 30 September 2010 (e.g., Fig. 1). Volcanic plume altitude affects the OMI SO₂ retrievals (Yang et al., 2007), and so OMSO2 products are provided for several assumed altitudes (Carn et al., 2013), including lower troposphere (∼3 km altitude) and mid-troposphere (∼8 km). For plumes with a designated altitude in the Smithsonian Institution volcanic activity reports (Global Volcanism Program; Venzke, 2013), the OMSO2 product closest to the reported altitude was used. In all other cases the lower troposphere SO₂ product was used, assuming that passively degassed SO₂ would rise no higher than 3–4 km at the target volcanoes. Daily SO₂ loading was calculated for each target volcano in a 2° × 2° region centered on the vent. Where the 2° × 2° regions for 2 or more volcanoes overlapped, a larger box was established with boundaries 1° from the edges of the volcano cluster. This resulted in the clustering of Ambrym and Gaua (Fig. 1), Rabaul and Ulawun, Fuego and Pacaya; and Merapi, Semeru, and Bromo; in these cases multiple sources of SO₂ degassing were considered as a single source. However, recent work (Carn et al., 2017) indicates that the dominant SO₂ sources in these regions are Ambrym, Rabaul, Fuego, and Semeru-Bromo (emissions from the latter are essentially indistinguishable in satellite data), respectively. The daily SO₂ mass was smoothed using a 1 week moving average to mitigate the impact of daily variations in OMI viewing geometry (e.g., Flower et al., 2016), and OMI ultraviolet reflectivity data were also used to screen excessively cloudy SO₂ observations.

Our earthquake catalog (69 events) was derived from the U.S. Geological Survey global catalog (https://earthquake.usgs.gov/earthquakes/search/), and includes all moment magnitude (M_w) ≥ 7 earthquakes at depths of ≤35 km between 1 January 2004 and 30 September 2010 (Table DR2 in the GSA Data Repository). In addition, all M_w ≥ 8 events during the period are included, regardless of depth. We use Computer Programs in Seismology (CPS; Herrmann, 2013) to model surface wave amplitudes. CPS produces synthetic seismograms based on a layered Earth velocity model (ak135-F, http://ds.iris.edu/ds/products/emc-ak135-f/), the focal mechanism, depth, and M_w of the earthquake, and incorporates the distance and azimuth of the target from the epicenter. To validate our synthetic waveform data set, we compared a subset of synthetic seismograms to real data at a variety of azimuths and distances. We found that surface wave amplitudes in the synthetic seismograms were typically within a factor of 2–3 of the real data, consistent with other studies (e.g., van der Lee, 1998). Although our modeling was relatively simple and did not account for finite fault complexity (i.e., directivity), we do not feel that more sophisticated modeling would significantly improve the discrepancy between the synthetic and real seismograms.

The OMI-derived SO₂ mass was averaged in two 14 day windows, one preceding and one following each earthquake, including the day of the event (e.g., Fig. 1). For each earthquake, the ratio of post-earthquake average SO₂ mass to pre-earthquake average SO₂ mass was calculated. The relationship between earthquake occurrence and SO₂ mass ratio was evaluated with variable thresholds of 1.2–2.0 (Tables DR3–DR9). The inverse of each SO₂ mass ratio was also calculated, to investigate if the PDS changes caused a decrease in volcanic degassing (i.e., a negative response). These calculations were repeated for PDS from 5 to 100 kPa, to investigate the effect of different levels of dynamic stress. To minimize the influence of meteorological clouds on the SO₂ data, we rejected days with a scene average reflectivity value >35%; in addition, if fewer than 8 days in either window met the reflectivity criterion, the earthquake was rejected.

In order to compare the occurrence of positive response to an earthquake to coincidental occurrence of elevated SO₂ emissions, we repeated the same windowing analysis for each day of the study period. This revealed any notable deviation in occurrence of response to an earthquake from that expected from random coincidence. We also further examined the variation in SO₂ response of two volcano clusters (Ambrym-Gaua and Merapi-Semeru-Bromo) by stacking the SO₂ emissions in a 31 day window centered on each qualifying earthquake generating PDS > 5 kPa. This provides a visualization of SO₂ degassing patterns before and after each earthquake and reveals how degassing re-equilibrates after an earthquake.

PDS values calculated using modeled seismograms contain more useful information than moment magnitude because the effects of focal mechanism and distance are included. For example, there are instances where the PDS generated by a smaller earthquake is higher for specific volcanoes than the PDS generated by a larger earthquake (e.g., Fig. DR1). Notably, an M_w 7.3 event generated a higher PDS than all other earthquakes up to M_w 8.0 (Fig. DR1). Full results of the 14 day SO₂ data window analysis for each target volcano or volcano cluster are reported in Tables DR3–DR9, and summarized in Figure 2.

There are at least two challenges to assessing the volcanic degassing response to earthquakes with our method: the frequency of major earthquakes at shallow depths is relatively low, and we expect any SO₂ degassing response to be modest and so not always detectable from space. However, considering our target volcanoes as an ensemble, two robust results emerge (Fig. 2): open-vent basaltic volcanoes show a clear occurrence of elevated SO₂ output following earthquakes that generated high PDS (positive response; e.g., Ambrym, Villarrica), and volcanoes undergoing cyclic activity that may include conduit plugging by more viscous magmas seem to demonstrate reduced SO₂ degassing following such earthquakes (negative response; e.g., Merapi, Fuego, Bagana). In this study, only Rabaul-Ulawun did not show a measurable degassing response to earthquake-induced dynamic stress; this could provide insight into the state of their magmatic plumbing systems.

Here we focus on our most significant results, i.e., the positive response of Ambrym-Gaua and the negative response of Merapi-Semeru-Bromo (Fig. 2). Ambrym-Gaua showed increased SO₂ emissions following 63% of earthquakes generating PDS of at least 5 kPa (Fig. 2), whereas higher SO₂ would be expected following only 41% of days selected at random. A stacked SO₂ emission plot (Fig. 3) indicates that SO₂ degassing at Ambrym (the dominant SO₂ source) consistently increased in response to earthquake-generated PDS. Ambrym hosted at least one active lava lake during the study period, likely maintained by the convection of low-viscosity basalt within the lake and conduit, driven by degassing-induced density variations (e.g., Stevenson and Blake, 1998). This suggests that a steady-state volcanic system that can maintain active lava lakes is more easily perturbed by small dynamic stresses (5 kPa; Fig. 3), and that the system responds by degassing SO₂ as it re-equilibrates.

Given that our calculated PDS values are likely far below the supersaturation pressures required for bubble nucleation (a few MPa in silicic magmas; e.g., Hurwitz and Navon, 1994), possible mechanisms to explain earthquake-triggered degassing include rectified diffusion (pumping of gas into preexisting bubbles; e.g., Brodsky et al., 1998) or sloshing (or another mass redistribution process) in a bubbly magma body (e.g., Carbone et al., 2009; Namiki et al., 2016). Disregarding rectified diffusion as inefficient (e.g., Ichihara and Brodsky, 2006), and given that the estimated PDS is commensurate with the stresses expected to weaken crystal-rich magma (e.g., Sumita and Manga, 2008), we favor a scenario whereby pockets of accumulated bubbles in the plumbing system are disrupted by dynamic strains, either mobilizing bubbles or increasing permeability via bubble wall collapse. Stacked plots of SO₂ emissions (Fig. 3A) indicate that peak SO₂ degassing occurs after a few days, suggesting that SO₂ release at the surface lags behind the triggering event (however, we have no constraints on the behavior of other volatiles, e.g., H₂O).

Merapi-Semeru-Bromo exhibited a negative response at all PDS levels (e.g., Fig. 2). This suggests minimal susceptibility to the effect of dynamic stress, at least with respect to degassing activity detectable from space. Figure 3B indicates a consistent decrease in degassing following each earthquake, with small post-event spikes that are less than pre-earthquake SO₂ emissions. We attribute this to the influence of higher viscosity basaltic-andesite or andesite magmas at these volcanoes, which typically erupt lava domes. In such systems, bubble rise is inhibited and degassing occurs primarily along permeable fracture networks in the conduit (e.g., Edmonds et al., 2003). Earthquake-induced stresses could act to modify or close existing fractures, reducing SO₂ emissions until the degassing pathway is reestablished. The cyclic nature of many dome-forming eruptions (e.g., Melnik and Sparks, 2005) may also be a dominant control on degassing, rendering them less susceptible to earthquake-induced effects. It is notable that Bebbington and Marzocchi (2011) found no statistical evidence for significant earthquake triggering of eruptions at Merapi, Semeru, or Bromo.

The earthquake sample size was lower for our other target volcanoes (Tables DR3–DR9), but the results concur with our proposed interpretation. Villarrica hosts an active lava lake and shows a positive degassing response (Fig. 2). Persistent, open-vent SO₂ emissions at Turrialba (e.g., Campion et al., 2012), coupled with its positive response, point to the involvement of low-viscosity magma supplying the emitted gases. This demonstrates the insight into volcanic systems that analysis of their response to earthquakes can yield.

Negative responses similar to Merapi-Semeru-Bromo were also observed at Fuego-Pacaya and Bagana (Fig. 2). The near-surface viscosity of conduit plugs and lavas likely exerts a strong control on activity at both Fuego and Bagana (e.g., Lyons et al., 2009; McCormick et al., 2012), with permeability governed by transient pore networks and cracks that could be susceptible to dynamic stresses (e.g., Heap and Kennedy, 2016). Thus in these cases permeability controls on degassing are likely dominant. The only neutral result in our analysis was found at Rabaul-Ulawun, which displayed no clear degassing response. There are many possible reasons for this, including opposing responses at each volcano (canceling out any net degassing signal), or the effects of static stress changes or fault rupture directivity (e.g., Delle Donne et al., 2010).

Our results indicate a striking dependence of the volcanic degassing response to earthquake-induced stresses on magma viscosity and eruptive style. We therefore recommend that interpretation of SO₂ emission time series at actively degassing volcanoes takes into account the potential effects of large earthquakes on observed trends. More significant insights into the nature of volcanic responses to earthquake-induced stresses are sure to arise from the coming years of observations of volcanic activity from space.

We acknowledge R. Hermann for his help with the CPS (Computer Programs in Seismology) code. Ben Kennedy and Michael Manga provided constructive reviews that greatly improved the paper. This work was supported by the Michigan Space Grant Consortium (Avouris) and NASA (Aura Science Team grant NNX11AF42G and MEaSUREs grant NNX13AF50G to Carn).