Gold Open Access: This paper is published under the terms of the CC-BY license.

Abstract

Statistical analysis of temporal relationships between large earthquakes (Mw ≥ 7) and volcanic eruptions suggests that seismic waves may trigger eruptions over great (>1000 km) distances from the epicenter, but a robust relationship between volcanic and teleseismic activity remains elusive. Here we investigate the relationship between dynamic stresses propagated by surface waves and a volcanic response, manifested by changes in sulfur dioxide (SO2) emissions measured by the spaceborne Ozone Monitoring Instrument (OMI). Surface wave amplitudes for a catalog of 69 earthquakes in A.D. 2004–2010 are modeled at 12 persistently degassing volcanoes detected by the OMI. The volcanic response is assessed by examining daily OMI SO2 measurements in 28 day windows centered on earthquakes meeting a variable peak dynamic stress threshold. A positive volcanic response is identified if the average post-earthquake SO2 mass was at least 20% larger than the pre-earthquake SO2 mass. We find two distinct volcanic responses, correlating strongly with eruption style. Open-vent, basaltic volcanoes exhibit a positive response to earthquake-generated dynamic stress (i.e., the earthquake triggers increased SO2 discharge), and andesitic volcanoes exhibit a negative response. We suggest that the former is consistent with disruption or mobilization of bubbles, or magma sloshing, in low-viscosity magmas, whereas the latter observation may reflect more dominant controls on degassing in viscous magmas or a post-earthquake reduction in permeability. Overall this analysis suggests that the potential effects of large earthquakes should be taken into account when interpreting trends in volcanic gas emissions.

INTRODUCTION

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 (SO2) 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. SO2 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 SO2 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 SO2 emissions permits such an analysis on a large scale for the first time.

DATA AND METHODS

To quantify daily volcanic SO2 emissions, we use the operational OMI SO2 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 SO2 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 Repository1). 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 SO2 loading measured by OMI between 1 October 2004 and 30 September 2010 (e.g., Fig. 1). Volcanic plume altitude affects the OMI SO2 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 SO2 product was used, assuming that passively degassed SO2 would rise no higher than 3–4 km at the target volcanoes. Daily SO2 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 SO2 degassing were considered as a single source. However, recent work (Carn et al., 2017) indicates that the dominant SO2 sources in these regions are Ambrym, Rabaul, Fuego, and Semeru-Bromo (emissions from the latter are essentially indistinguishable in satellite data), respectively. The daily SO2 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 SO2 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 (Mw) ≥ 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 Mw ≥ 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 Mw 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.

Synthetic seismograms were computed for each earthquake in our catalog, using the target volcano locations as receiver positions. We used the synthetic seismograms to estimate the PDS with the simplified shear-wave method of Velasco et al. (2004): 
graphic
where v is ground velocity, µ is rigidity, and β is shear wave velocity. We used average crustal values of 3.3 × 1010 Pa for µ and 3.5 km/s for β (Velasco et al., 2004). Rather than calculate separate PDS for Rayleigh and Love waves, the three components of each synthetic seismogram were combined into a single velocity vector, from which peak velocity was derived for the time period spanning the arrivals of 40 s to 15 s period waves. Clearly our method does not properly account for Rayleigh wave PDS, but serves as a quantitative proxy for total surface wave PDS.

The OMI-derived SO2 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 SO2 mass to pre-earthquake average SO2 mass was calculated. The relationship between earthquake occurrence and SO2 mass ratio was evaluated with variable thresholds of 1.2–2.0 (Tables DR3–DR9). The inverse of each SO2 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 SO2 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 SO2 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 SO2 response of two volcano clusters (Ambrym-Gaua and Merapi-Semeru-Bromo) by stacking the SO2 emissions in a 31 day window centered on each qualifying earthquake generating PDS > 5 kPa. This provides a visualization of SO2 degassing patterns before and after each earthquake and reveals how degassing re-equilibrates after an earthquake.

RESULTS AND DISCUSSION

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 Mw 7.3 event generated a higher PDS than all other earthquakes up to Mw 8.0 (Fig. DR1). Full results of the 14 day SO2 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 SO2 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 SO2 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 SO2 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 SO2 emissions following 63% of earthquakes generating PDS of at least 5 kPa (Fig. 2), whereas higher SO2 would be expected following only 41% of days selected at random. A stacked SO2 emission plot (Fig. 3) indicates that SO2 degassing at Ambrym (the dominant SO2 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 SO2 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 SO2 emissions (Fig. 3A) indicate that peak SO2 degassing occurs after a few days, suggesting that SO2 release at the surface lags behind the triggering event (however, we have no constraints on the behavior of other volatiles, e.g., H2O).

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 SO2 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 SO2 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 SO2 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 SO2 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).

1GSA Data Repository item 2017230, additional results, Tables DR1–DR9, and Figures DR1and DR2, is available online at http://www.geosociety.org/datarepository/2017/ or on request from editing@geosociety.org.