Dynamic earthquake triggering response tracks evolving unrest at Sierra Negra volcano, Galápagos Islands.

The propensity for dynamic earthquake triggering is thought to depend on the local stress state and amplitude of the stress perturbation. However, the nature of this dependency has not been confirmed within a single crustal volume. Here, we show that at Sierra Negra volcano, Galápagos Islands, the intensity of dynamically triggered earthquakes increased as inflation of a magma reservoir elevated the stress state. The perturbation of short-term seismicity within teleseismic surface waves also increased with peak dynamic strain. Following rapid coeruptive subsidence and reduction in stress and background seismicity rates, equivalent dynamic strains no longer triggered detectable seismicity. These findings offer direct constraints on the primary controls on dynamic triggering and suggest that the response to dynamic stresses may help constrain the evolution of volcanic unrest.

INTRODUCTION

Dynamic stress perturbations result from the radiation of seismic waves from a source rupture (, ). Earthquakes triggered by dynamic stress perturbations (–) have been observed thousands of kilometers away from the source (–), well beyond the reach of static stress perturbations (, ). Although the mechanisms underpinning dynamic triggering remain enigmatic (), earthquakes are more likely to be triggered on faults already close to failure (). Therefore, as well as playing an important role in controlling earthquake sequences (), dynamic triggering also offers a means to probe the stress state of the subsurface (, ). To do this, however, it requires knowledge of the triggering response as a function of stress state and stress perturbation. Efforts to quantify this behavior have relied on statistical studies of the response to many triggering earthquakes in many different crustal volumes (, –). As it is difficult to measure directly, seismicity rates are commonly used as a proxy for the state of stress, with high seismicity rates corresponding to a large number of faults close to failure (). However, generally, the background stress state (shear stress or pore fluid pressure) does not change substantially during the observation window. Consequently, it is challenging to undertake a longitudinal study of the same region over an extended period at different stress states.

Dynamic triggering of seismicity at volcanoes is well known (–, –). As some volcanoes experience rapid, high-amplitude changes in stress due to magma accumulation and withdrawal, they potentially offer good conditions to study the factors that influence dynamic triggering. However, the non-Poissonian, swarm-like, nature of much volcanic seismicity means that triggering responses can be difficult to detect and quantify. So far, no evidence has been reported for a systematic dependence of the triggering response at volcanoes on either the background earthquake rate or amplitude of the dynamic stress perturbations ().

Sierra Negra is a large basaltic shield volcano in the Galápagos Islands, located 1000 km west of the coast of mainland Ecuador () (Fig. 1). It has a prominent summit caldera and an intracaldera “trapdoor” fault system (). Recent eruptions are associated with a deformation cycle controlled by inflation and deflation of a sill located 2 km below the caldera floor (). Eruptions are preceded by years of surface uplift () and accompanied by rapid subsidence (). The 26 June to 25 August 2018 eruption () was preceded by >6.5 m of uplift in the 13 years from the end of the 2005 eruption and was accompanied by a total of 8.5 m of subsidence (Fig. 2A). Uplift rates in the final 6 months before the eruption exceeded 1 m/year. A seismometer at site VCH1 on the northern caldera rim (Fig. 1C) recorded large numbers of volcano-tectonic earthquakes through this deformation cycle (Fig. 2B). High-precision source locations obtained using a temporary deployment of seismometers from late April 2018 () indicate that almost all earthquakes occurred on the trapdoor fault system, at depths of <2 km (Fig. 1C). Earthquake rates increased with increases in uplift rate and total uplift, increasing notably in the year before the eruption (Fig. 2, B and C). Coeruptive earthquake rates were initially very high and decreased as deflation slowed. The polarity of earthquake focal mechanism flipped as deflation began, indicating a reversal of the stress field (). At the end of the eruption, background earthquake rates had fallen to well below those at the start of monitoring in 2010.

Fig. 1.

Location of Sierra Negra volcano, Galápagos Islands.

(A) Volcanoes of the Western Galápagos islands of Isabela and Fernandina and permanent monitoring network of Instituto Geofísico de la Escuela Politécnica Nacional (IGEPN). (B) Location map showing position of the Galápagos archipelago 1000 km west of the coast of mainland Ecuador, and major tectonic features. (C) Epicenter locations of earthquakes recorded by the IGUANA network, 22 April to 26 June 2018. Location of station VCH1 indicated by red triangle.

Fig. 2.

Geophysical unrest at Sierra Negra, 2010–2020.

(A) Total inflation at the center of the caldera measured by cGPS. Vertical green dashed line indicates time of eruption onset. (B) Magnitudes (blue circles) and total number of M > 0.8 earthquakes (red line) detected at VCH1. (C) Thirty-day “background” rate of earthquakes. Scale is clipped at 1200 to exclude very high rates of earthquakes accompanying the eruption.

Interferometric Synthetic Aperture Radar (InSAR) and continuous Global Positioning System (cGPS) data suggest that pre-eruptive deformation at Sierra Negra is predominantly elastic and driven by pressure changes in the shallow sill, which is effectively in permanent hydrostatic connection with a deeper reservoir (). Accompanying seismicity reflects the evolution of local stress conditions on the weak trapdoor fault system. This process results in a smoothly changing seismicity rate, with statistics more closely approximating a Poisson process than observed in many other volcanic settings (fig. S1).

The proximity of the Galápagos Islands to eastern Pacific subduction zones means that they experience a high rate of large dynamic stress perturbations (fig. S3). Therefore, as seismicity at Sierra Negra is characterized by large systematic fluctuations in rate and relatively simple statistics, this makes it an ideal setting to study processes of dynamic earthquake triggering at volcanoes. In this study, we identify episodes of dynamic triggering at Sierra Negra and then seek to quantify the triggering response to different amplitude stress perturbations at different volcanic unrest states.

RESULTS

For the largest dynamic stresses experienced at Sierra Negra, there is clear evidence of earthquake triggering. For each potential trigger event, we quantify the amplitude of stress perturbation using peak dynamic strain, “PDS” (, , ), based on measured ground velocities during teleseismic arrivals at VCH1 in a 30- to 2-s band. Figure 3 illustrates dynamic triggering of local seismicity due to the moment magnitude (Mw) of 8.2 Chiapas earthquake, Mexico. The unfiltered data (Fig. 3, A to C) show the high-amplitude, lower-frequency arrivals associated with the teleseismic earthquakes, with a peak ground velocity in the Rayleigh waves of 2 mm/s. The dynamic strain (, , ) in a 30- to 2-s band (Fig. 3D) increases during the body wave arrivals and again during the surface wave arrivals reaching a peak value of 10–6.4. When high pass–filtered at 3 Hz, the lower-amplitude, higher-frequency local seismicity is evident (Fig. 3, E and F). The highest seismicity rate coincides with highest velocity surface wave arrivals, although, in this instance, locally triggered seismicity is also evident at the peak velocity of the body wave arrivals. The 13 local earthquakes of Ml > 0.8 in the 30 min following the surface wave arrivals reflect a ~30-fold increase over the background rate at this time.

Fig. 3.

Seismic data recorded at station VCH1 during arrivals of Mw 8.2 Chiapas earthquake, Mexico, 8 September 2017.

(A) Vertical (HHZ), (B) radial (HHR), and (C) tangential (HHT) components. P-wave and surface wave arrivals indicated by green and white dotted lines, respectively. (D) Dynamic strain. (E) Data high pass–filtered at 3 Hz. Dashed vertical light blue lines indicate times of local earthquakes in our template-matched and manually picked catalog (Methods and Materials). (F) Spectrogram of data.

Figure 4 shows spectrograms for teleseismic arrivals recorded at Sierra Negra with 10 of the largest dynamic stress perturbations between 2010 and 2020 in order of time: 7 before the eruption and 3 after the eruption. Figure 5 shows the same data in the form of dynamic strain, including the times of local Ml > 0.8 earthquakes. Triggered seismicity is most intense during the highest amplitude surface wave arrivals and decays quickly after this peak. By inspection, the number of triggered earthquakes apparently increased toward the eruption (Fig. 4, A to H), and then is absent afterward (Fig. 4, I and J).

Fig. 4.

Spectrograms of the vertical component of 10 of the largest PDSs associated with surface wave arrivals recorded at station VCH1 during 1 January 2010 to 30 July 2020.

(A to J) P-wave and surface wave arrivals indicated by green and white dashed lines, respectively. All spectrograms use the same color scale. Note the prominent T-phases in (B), (C), (E), (I), and (J).

Fig. 5.

Dynamic strain of 10 of the largest PDSs associated with surface wave arrivals recorded at station VCH1 during 1 January 2010 to 30 July 2020.

(A to J) P-wave and surface wave arrivals indicated by green and white dotted lines, respectively. Catalog of template-matched and manually picked local earthquakes indicated by vertical dotted black lines.

We investigate the triggering response for any potential triggering event with a PDS > 10–7.5 but exclude the few moderate-size earthquakes within the Galápagos Islands that might be associated with static stress changes at Sierra Negra, leaving a catalog of 60 events. This selection makes the manual catalog picking viable and includes all events for which we expect to see a triggering response. The highest PDS experienced was 10–5.8 during the 2016 Mw 7.8 Pedernales earthquake (1100 km from Sierra Negra). We build a catalog of local high-frequency earthquakes for Sierra Negra based on detections at VCH1 (Materials and Methods). We use template matching and manual picking to identify closely spaced earthquakes that occurred during the triggering window and define a triggered earthquake rate using the number of events that occur in a 30-min window after the surface wave arrivals. The difference between the triggered seismicity rate and the 30-day background rate is then quantified in terms of a number of SDs using the β-statistic (, ). The β-statistic assumes a Poisson distribution of earthquake rates, which is more appropriate at Sierra Negra than other systems (). We determine 95% confidence limits for the β-statistic as a function of time, using a bootstrap method based on 500,000 random samples from the earthquake catalog ().

Figure 6A shows the β-statistic between 2010 and 2020. Thirteen events are associated with values of β greater than the 95% confidence limit, all but one before the onset of the 2018 eruption, although some with low absolute numbers of earthquakes in the triggering window. In Figure 6B, the rate of triggered earthquakes is generally higher for events with larger PDS. In particular, events in later 2017 and early 2018 with PDS ~10–6.5 are associated with higher rates than larger events in 2010 and 2016. Figure 6C illustrates the response in terms of the β-statistic, thus accounting for changing background seismicity rate. Above a PDS of ~10–6.8, the β-statistic apparently increases linearly with log10PDS. We can see that the higher triggered seismicity rates during the 2017 and 2018 earthquakes are due to the combination of high PDS and the elevated background seismicity rates as the eruption approached. At lower PDS, some events have values of the β-statistic above 0 but less than 10. These are associated with single earthquakes during the 30-min triggering window.

Fig. 6.

Statistics of dynamic triggering at Sierra Negra from 2010 to 2020.

(A) β-statistic of seismicity associated with large teleseismic earthquakes. Size of circle indicates PDS during surface waves arrivals. Horizontal black lines represent the time-dependent bootstrapped 95% confidence limit for β-statistic. Light gray shading indicates seismic data gaps. Color indicates time through the sequences, with red points after the end of the eruption on 25 August 2018. Vertical red dashed line indicates onset of eruption. (B) Earthquake rate and (C) β-statistic during surface wave window as a function of PDS. Colors correspond to those in top panel. (D) Predicted (circles) and observed (blue line) 30-day background rate.

At the time of writing, there have been three earthquakes associated with PDS exceeding 10–6.8 since the end of the eruption on 25 August 2018. No triggering is evident for any of these events by inspection (e.g., Figs. 4 and 5, H to J).The largest PDS after the eruption was 10–6.5 associated with the 26 May 2019 Mw 8.0 Lagunas earthquake in Peru (Figs. 4 and 5H). The density of data from available observations is not sufficient to confidently constrain any particular relation between PDS and β. However, if we assume a linear relation between log10PDS and the β-statistic consistent with the available data, then we would predict a β-statistic of ~10 for this event. With the low rates of seismicity at that time, this value of β would correspond to a predicted 30-min earthquake rate of 1.4 events, and the probability of observing no earthquakes is 25%. Consequently, our observation of no triggered earthquakes for this event would be consistent with a constant relation between β and PDS before and after the eruption.

If we assume this speculative linear relation between log10PDS and the β-statistic, then using the observed number of triggered earthquakes for a known PDS, we can make predictions as to what would be a corresponding background rate of earthquakes. This relation needs more data points to be sufficiently constrained, but an example prediction of background rates, and the observed values, are shown in Fig. 6D.

DISCUSSION

Our observations suggest that dynamic stress perturbations with larger peak dynamic strains resulted in a larger triggering response in terms of the rate change of volcano-tectonic seismicity, measured by the β-statistic. As the β-statistic is a measure of rate change relative to the background rate of seismicity, this observation also implies that, for a given dynamic stress perturbation, a greater absolute number of earthquakes are triggered if the background rate of seismicity is high (e.g., immediately before the eruption) than if it is low (e.g., after the eruption). We find no evidence to suggest that the physics of the triggering process changed after the deflation and relaxation accompanying the eruption; the absence of triggering after the eruption is entirely consistent with a drop in the background rate of seismicity, as fewer faults are close to failure. To our knowledge, this is the first time that such a relation has been observed at a volcano. In this context, we argue that the triggering response probably tracks the stress state of Sierra Negra. Hence, observations of dynamic triggering at volcanoes could provide additional evidence for elevated volcanic unrest, alongside more dependable measures such as the rate of seismicity and deformation at a volcano.

We only observe convincing evidence for dynamic earthquake triggering at Sierra Negra for events with PDS > 10–6.8 and, even then, only when the rates of seismicity are elevated before the eruption. These are quite large strains. For example, dynamic triggering has been reported at more than two orders of magnitude smaller strains albeit using a statistically aggregated interevent time ratio within a large regional catalog (). We have no evidence whether 10–6.8 is a real physical threshold below which triggering does not occur at Sierra Negra, or, more likely, a detection limit dependent on the background rate of seismicity. Such evidence will statistically emerge as many more moderate-sized triggering events are documented, ideally at high seismicity rates during the approach to eruption. Although we suggest a potential logarithmic relation between β and PDS, more data are needed to test this model, and alternatives (such as a step function) cannot be excluded from these available data. We await more large earthquakes to test and build on our observations here.

All the evidence for earthquake triggering at Sierra Negra is associated with the highest surface wave (and exceptionally, body wave) amplitudes and their immediate aftermath. We see no evidence for delayed triggering in these data. The most intense triggering coincides with the highest velocities recorded during the Rayleigh wave arrivals, but there is evidence for triggered seismicity in the Love waves and body waves. Future studies of the details of the timing of triggering may help to better constrain its physical mechanisms. Cross-correlation of the waveforms of triggered earthquakes with those of nontriggered seismicity as part of the template matching process shows that many of the triggered earthquakes belong to the same waveform families as earthquakes that happened many months before or after the triggered event (e.g., fig. S4). These repeating earthquakes represent recurrent failure of the same or similar patch of the trapdoor fault. Consequently, we interpret the triggered earthquakes as occurring on critically stressed faults that would probably have slipped eventually without the additional dynamic stresses. However, the detailed physics of the triggering process remain unclear and are likely to be more complex than simply a shear stress increase on a fault close to failure.

The prolonged accumulation of magma in the shallow edifice before the 2018 eruption drove predominantly elastic deformation, probably increasing the shear stress on faults within the weak trapdoor fault system (). These elevated shear stresses meant that more faults were close to failure, and seismicity rates were high. Shear stresses relaxed during the eruption, and the posteruption seismicity rate was low (and remains low at the time or writing). Although pore pressure is clearly important in controlling seismicity, there is no indication that the hydrothermal system underwent substantial systematic changes at Sierra Negra before or after the eruption. Therefore, to first order, we believe that changes in shear stress on the trapdoor fault system are controlling the evolving background seismicity rates at Sierra Negra.

It is not clear whether Sierra Negra is particularly sensitive to dynamic triggering, i.e., whether the values of the β-statistic are particularly high for a given PDS compared to other volcanoes or tectonic areas. During the same period, no triggering responses were observed at any of the other Galápagos Island volcanoes. However, at the most likely candidate volcanoes, seismicity rates were generally lower than at Sierra Negra, and the nearest seismometers are located >10 km away from the edge of the caldera (fig. S2), so any triggered earthquakes might not be detectable. We hope that this work motivates future systematic studies of dynamic earthquake triggering in volcanic settings, including the deployment of dense monitoring networks to detect and locate low-magnitude seismicity.

MATERIALS AND METHODS

Seismic data

Because of the remote location of Sierra Negra, the extent of seismic monitoring at the volcano has been limited. A permanent regional network operated by the Instituto Geofísico de la Escuela Politécnica Nacional (IGEPN) was installed in late 2012 (, ) and used to construct a seismic catalog for the Galapagos. The IGEPN broadband station “VCH1” is located at Volcan Chico on the northeast rim of the caldera and provides an excellent record of seismicity at Sierra Negra (Fig. 1A). The station GS09 of the SIGNET campaign deployment () was located at the same site as VCH1, extending the record of seismicity back to early 2010. In late April 2018, a dense campaign network (IGUANA) was installed at Sierra Negra, recording the seismicity associated with final approach to eruption, the eruption and deflation, and renewed inflation. Gaps exist in the data between 2011 and 2012, after the SIGNET campaign deployment and before the IGEPN permanent network was installed. Substantial gaps due to instrument or network failure are present in 2013, 2015, and 2016. The 2018 eruption began on 26 June and finished on 25 August. A continuous GPS network recorded the surface deformation across the caldera through this time (Fig. 2A).

The large distances between stations in the IGEPN seismic network means that earthquake locations have substantial uncertainty and the detection threshold is high. Most small earthquakes are only recorded on the nearest station. Since 2010, seismicity in the Galapagos was numerically dominated by events located at Sierra Negra, apart from a small number of short-lived episodes associated with eruptions at La Cumbre volcano, Fernandina Island (), and an intrusion at Cerro Azul volcano in 2017 (fig. S2). Outside these times, high-frequency earthquakes detected at VCH1 are almost all associated with the unrest at Sierra Negra. Therefore, VCH1 provides an excellent record of seismicity at Sierra Negra, even though it is not possible to determine precise locations. Almost all seismicity recorded by the IGUANA deployment is located on the trapdoor fault system within the caldera (Fig. 1C). The characteristics of seismicity recorded at VCH1 from 2010 are similar to those recorded during the IGUANA deployment (frequency content, s-p times, etc.), supporting our contention that most of the seismicity detected at VCH1 from 2010 is also located on the trapdoor fault system.

We initially identify the 500 largest candidate trigger earthquakes at Sierra Negra in the global U.S. Geological Survey catalog using a surface wave amplitude-distance relation (fig. S3). We exclude events within 200 km of the volcano to avoid any potential static stress triggering effects. For each candidate event, we then measure the actual peak ground velocities for the surface waves at VCH1 in a 30- to 2-s band. We use these velocities to estimate peak dynamic stress and strain () associated with each potential trigger event. During method development, we compare strain estimates based on single components of the data (radial, transverse, and vertical), root mean square as opposed to peak values, and measures of displacement, velocity and acceleration. In all cases, we find similar results to PDS and, so for consistency with previous studies (, , ), use this as our measure of the triggering amplitude. We select the 60 events with the largest PDSs (PDS > 10–7.5).

We use an Short-term Average/Long-term Average (STA/LTA) algorithm, implemented in ObsPy (), to detect local high-frequency earthquakes over 10 years of data at GS09 and VCH1. We first bandpass the data between 5 and 10 Hz. To define an event, we require a trigger on all three components. We then remove any events with anomalous frequency content and very short trigger durations. We also remove any events with very short interevent times (<10 s), within the typical waveform duration. We have visually verified the high reliability of the algorithm on hundreds of events.

To supplement the STA/LTA detections during high rates of seismicity coincident with high-amplitude surface wave arrivals, we use a template-matching methodology (, ) to identify masked or overlapping events. We use ~15,000 10-s duration templates from our detected event catalog (not considering the coeruptive events where the stress field reversed). We cross-correlate each of these events with a 1-hour window encompassing each stressing event (30 min before and 30 min after the surface wave arrivals). We specify a cross-correlation value of >0.6 and an interearthquake separation of >5 s to identify an event. We then manually pick these data, inspecting each STA/LTA and template-matched pick and adding new picks where appropriate. The template-matching method does detect some additional events during the highest rates of seismicity. Manual picking adds more and removes a small number of erroneous picks. Most of the detected events are small and below the magnitude threshold of 0.8.

For each event, we determine the root mean square amplitude of the vertical component of the velocity in a 10-s window around the peak amplitude, in a 1- to 20-Hz frequency band. We then estimate a magnitude, based on a linear calibration of log10 (amplitude) against the magnitude of located events in the IGUANA catalog, incorporating an empirical relation between amplitude attenuation and epicentral distance.

We estimate a completeness magnitude of 0.8 for our catalog of high-frequency earthquake detections at Sierra Negra (fig. S1, C and D). We apply a range of methods for estimating the completeness magnitude, including the b value stability test (, ). During the highest rates of seismicity before the eruption, the magnitude-frequency distributions follow a Gutenberg-Richter relation, allowing reliable estimation of an instrumental completeness magnitude. Earlier data are not consistent with a Gutenberg-Richter relation, and so standard methodologies for estimating Mc are not reliable [we note that the apparent reduction in b value before the eruption is consistent with an apparent increase in stress; ()]. However, as instrumentation is the same, and rates of seismicity are lower, we conclude that 0.8 is a reliable estimate for the duration of our data (excluding the coeruptive period). The value of 0.8 is also consistent with day/night event ratios. We note that we measure the triggering effect using the β-statistic—measure of relative rate changes. Consequently, small changes in catalog completeness unlikely result in large changes in β.

We quantify the triggering response using the absolute rate of triggered earthquakes and the β-statistic (). We use a 30-day window to determine the background rate, although we test shorter and longer time intervals and again find that they have little influence on our results. The rate of seismicity varies slowly at Sierra Negra during inflation.

The β-statistic quantifies the change in rate in terms of the number of SDs and is useful for comparing rates in windows of different durations ()Here, N is the number of earthquakes recorded in the period of duration T after the arrival of the surface waves. We define the surface wave arrival time using a group velocity of 4 km s−1 and confirm this by direct inspection of the data, although results are generally insensitive to small changes in this velocity. N is the number of earthquakes recorded in the time period T, comprising a window before the surface wave arrivals (used to estimate the background rate) plus the time period after the surface wave arrivals. The period T therefore includes the period T. In this study, we prefer T = 30 min and T = 30 days. These values are shorter than used previously () but consistent with the immediate triggering effect observed in our data. We test our results using different values of T from 5 min to 1 day and fixed numbers of events rather than time windows. For different combinations of Mc, T, and T, we find slightly different numbers of triggered and background earthquakes, but overall patterns of our results are robust and support the qualitative observations in Figs. 4 and 5.

We use a postulated logarithmic relation between β and PDS to “predict” a background rate of seismicity from the observed rate of triggered earthquakes and the PDS. We find the value of N that minimizes the misfit in the relation: β = 16 log10 PDS + 120, where β is given by Eq. 1.