Continuous seismic data have been the primary data stream used to forecast eruptive activity during volcanic crises since modern volcano monitoring began. Examples of forecasts from seismic data are numerous (e.g., Endo et al., 1981; Malone et al., 1983; Klein, 1984; Swanson et al., 1985; Power and Lalla, 2010; Chouet et al., 1994; Harlow et al., 1996; Soosalu et al., 2005; Power and Lala, 2010; Buurman and West, 2010; and; Ruppert et al., 2011). The most established seismic indicators of incipient eruption include increasing number and amplitude of seismic events and shifts in types of seismic events from VTs (volcano-tectonic) to LFs (low frequency) and tremor as magma nears the surface (e.g., Minakami, 1974; Voight, 1988; Chouet, 1996; McNutt, 1996; Boué et al., 2015; McCausland et al., 2019; White and McCausland, 2019). Many of these methods are most reliable for explosive eruptions, whereas small and slowly evolving eruptions can be more difficult to forecast reliably (e.g., Cameron et al., 2018) and require a detailed understanding of a volcano’s structure, tectonic framework, hydrologic system, and eruptive history.

Seismicity escalated dramatically in the months prior to the 2017–2018 eruption at Mount Agung, Bali, Indonesia, providing ample warning of an impending eruption (Syahbana et al., 2019) and raising fears of a violent eruption, similar to the VEI 5 in 1963 (Kusmandinata, 1964). Although the eruption was expected in the long term, forecasting the size and onset time of the eruption in the short term was difficult. Ultimately, the eruption was only a VEI 2, characterized by relatively small explosions and the extrusion of 24 million m³ of lava at the summit crater (Syahbana et al., 2019). Ash fall affected local farmland and small lahars traveled down drainages on the N and S flanks, and there were no fatalities.

Several challenges presented themselves throughout the crisis:

(1) Scientists with CVGHM (Indonesia’s Center for Volcanology and Geologic Hazard Mitigation), who were responsible for monitoring, had very few automatic processing tools to assist them in their work. Aside from RSAM (Realtime Seismic Amplitude Measurement: Endo and Murray, 1991), staff manually counted and classified earthquakes to interpret volcanic activity and provide information to decision makers and the public.

Small eruptions can be difficult to forecast precisely because precursory signals can be subtle, but in densely populated areas small eruptions still pose a great threat. Thus, new approaches and analysis tools are needed. Studies of repeating earthquakes, or earthquake families, have been successful, though often in retrospect, of showing notable changes in eruptive activity at a variety of volcanoes (e.g., Green and Neuberg, 2006; Umakoshi et al., 2008; Arambula-Mendoza et al., 2011; Thelen et al., 2011; Buurman et al., 2013; Rodgers et al., 2013). Repeating earthquake families can include either broad-band brittle failure or low-frequency earthquakes and are observed when seismic sources are stationary, non-destructive and of the same mechanism (Geller and Mueller, 1980). Processes that cause them vary widely, including brittle failure associated with deep pressurization and conduit building under the summit (Buurman and West, 2010; Deshon et al., 2010; Budi-Santoso and Lesage, 2016), shallow fluid movement and lava dome growth at the surface (Stephens and Chouet, 2001; Rowe et al., 2004; Green and Neuberg, 2006; Umakoshi et al., 2008; Matoza and Chouet, 2010; Thelen et al., 2011), magma-water interaction in the shallow system (Rodgers et al., 2013; Rodgers et al., 2015), and degassing (Caplan-Auerbach and Petersen, 2005; Waite et al., 2008; Matoza et al., 2015).

In this retrospective study, we examine the time history of earthquake families at Mount Agung in search of additional insight into the temporal changes in the shallow crust (upper 20 km) prior to eruption. Specifically, we analyze the period of declining seismic activity ∼5 weeks prior to the eruption. We show that the evolution of earthquake families illustrates that Mount Agung was progressing toward eruption even though overall earthquake rates and seismic-energy-release declined.

After a 50 + year period of repose, unrest at Mount Agung began in mid-July with anomalous levels of seismicity, thermal anomalies, and an increase in steaming at the volcano’s summit (Syahbana et al., 2019). Earthquake rates and magnitudes intensified in September, and reports of felt seismicity and locations for the largest earthquakes (>M2.3)—computed by Indonesia’s Meteorological, Climatological, and Geophysical Agency (BMKG), which runs the national seismic network—suggested the earthquakes were located between the edifices of Mount Agung and Batur (NW of Mount Agung) at depths between 10 and 20 km below sea level (Figure 1). The growing swarm of brittle failure VTs was interpreted as the seismic response of the crust as a deep intrusion under the summit of Mount Agung (White and McCausland, 2016). Later work helped clarify the path magma took to reach the surface. A series of InSAR images from 21 September 2017 UTC to 8 November 2017 UTC show a pattern of inflation that is consistent with a dike intrusion between 7 and 13 km depth on the NW flank of Mount Agung (Albino et al., 2018). The location of the magma reservoir below the depth of the dike is not known, but Syahbana et al. (2019) favor a model where a the dike propagated laterally from a deeper magma reservoir (12 to >15 km) beneath the summit of Mount Agung before magma continued to the summit.

One day before the peak seismic earthquake rates on 23 September 2017, CVGHM raised the alert level to IV (the highest of 4 levels). Earthquake magnitudes reached a new maximum (M4.2) a day later on 23 September 2017 (Syahbana et al., 2019), but subsequently earthquake rates and RSAM values decreased significantly (Figure 1). After a month of reduced earthquake rates and declining RSAM values, the alert level was lowered to III on 29 October 2017 (Syahbana et al., 2019).

During this time of uncertainty, CVGHM staff were vigilant for any signs of shallowing magma, such as the emergence of tremor and low-frequency (LF) seismicity or a shift of seismicity toward the summit (Chouet, 1996; White and McCausland, 2019). On 9 November 2017, a M4.9 earthquake, the largest of the crisis, and a series of aftershocks occurred in a prominent new location NE of the mountain at ∼10 km depth, according to BMKG. According to global solutions by the U.S. Geological Survey (USGS) National Earthquake Information Center ² and the Global Centroid Moment Tensor ³ (Dziewonski et al., 1981; Ekström et al., 2012), the earthquake occurred on an approximately E-W striking, south-dipping thrust fault—consistent with seismicity in the back-arc of the Sunda convergent margin (Supendi et al., 2020). On 12 November, CVGHM staff noted the first significant tremor as well as the first appearance of LF earthquakes. These tremor bursts were ∼40–120 s each, were broadband (1–10 Hz), and were barely strong enough to rise above the amplitudes of daytime cultural noise on the analog stations TMKS and PSAG (Syahbana et al., 2019). Together, the M4.9 earthquake and appearance of tremor and LFs were interpreted by CVGHM as a possible shift toward a more phreatic phase of precursory unrest as heat interacted with the shallow hydrologic system and fluids and gases moving toward the surface began opening the conduit. Later work by (Sahara et al., 2021) (submitted) lend credence to this interpretation by showing a gradual shift in earthquake locations toward the summit.

The fact that tremor did not continue and that shallow LF earthquakes were few in number, however, suggested that if processes in the volcanic conduit were progressing toward eruption, they were doing so slowly. The lack of convincing evidence for shallow conduit seismicity and the increased lag time since the energetic VT swarm increased CVGHM’s uncertainty that an eruption would occur in the near future (Syahbana et al., 2019), and CVGHM opted not to raise the alert level. In hindsight, another InSAR image from 20 November showed significant depletion from the inflation source under the NW flank, which Albino et al. (2018) interpret as a sign that magma had begun to migrate toward the summit. The first gas measurements, which were taken by unoccupied aircraft systems early on 21 November, also showed high levels of CO₂, which was interpreted as a significant sign of unrest (Syahbana et al., 2019).

Earthquake rates remained at low levels and fluctuations in RSAM stayed within normal bounds leading up to the first phreatomagmatic eruption on 21 November. LF earthquakes and tremor increased after the first phreatomagmatic eruption, and on 25 November, CVGHM staff documented 21 larger LF signals with a dominant frequency of ∼2 Hz over a 90-min period (Syahbana et al., 2019). These earthquakes were closely associated in time with the onset of lava at the summit, and within just a few days, 24 million m³ of lava had filled one-third of the summit crater, and regular explosions occurred until mid-January 2018 (Syahbana et al., 2019).

We used REDPy (Hotovec-Eills and Jeffries, 2016) to build an earthquake catalog at Mount Agung prior to, during, and after the late-November 2017 eruptions, and we group those earthquakes into families based on waveform similarity to assess changes in seismicity over time. We focus on the time period of greatest forecasting uncertainty—i.e., when seismicity decreased following peak rates in late September, and then in the period after the eruption started.

REDPy is designed to comb through real-time or archived continuous data to produce a catalog of triggered signals using a short-time-average/long-time-average (STA/LTA) algorithm. It then cross-correlates all triggered signals against each other, and earthquakes with cross-correlation coefficients above a certain threshold are grouped into families of repeating earthquakes. Earthquakes from a single family are presumed to have originated from a similar location and a similar source mechanism, thus providing some information about earthquake source properties even if the network is insufficient to compute hypocenters or focal mechanisms (Geller and Mueller, 1980). Earthquakes whose cross-correlation coefficients with prior events are below a defined threshold are labeled “orphans” and remain in the catalog to be compared to future earthquakes. We also use REDPy to compute the Frequency Index (FI) of each earthquake by comparing energy in a high frequency band (5–10 Hz) to energy in a low frequency band (1–2.5 Hz) (Buurman and West, 2010). Frequency index is an arithmetical method to describe earthquake frequency content, and in our case, we empirically define the boundary between VT and LF earthquakes at 0 (Figure 2).

Prior to the crisis, the seismic monitoring network at Mount Agung consisted of two short-period, vertical seismometers that were transmitted via analog telemetry back to the observatory post at Pos Rendang. The stations were ∼4 and 5 km away from the summit. Between mid-October and early November, CVGHM and VDAP (USGS-USAID Volcano Disaster Assistance Program) installed seven new seismometers — 6 broadband and 1 short-period—with digital telemetry back to the observatory post (Figure 1). Four of those stations were located <10 km from the volcano while the other three were located farther away. For this analysis, we used 5 stations (the two original short-period analog stations (TMKS, PSAG), the short-period digital station (ABNG) and two of the closest broadbands (CEGI, YHKR) (Figure 1). We did not include the other proximal broadband, DUKU, because it had significant data outages during our study.

We filtered all data between 1 and 10 Hz before applying the recursive STA/LTA algorithm (short-time widow: 3 s; long-time window: 8 s; trigger on threshold: 1.8; trigger off threshold: 1.3). We required a detection at 3 stations in order to produce a trigger, and we required a cross-correlation coefficient of 0.8 at 3 stations in order to be considered a matching waveform. Finally, we required families comprise five earthquakes or more for the events to be grouped as a family and counted as a repeater for the purposes of Figures 1–3 and Table 1.

We reviewed all earthquake waveforms in Swarm ⁴ to verify that regional earthquakes and noise were not included in our analysis. We also repeated our analysis with cross-correlation coefficients of 0.6 and 0.7 to verify that the interpretation of our results is not sensitive to this value (Supplementary Table S1). Finally, we extended our analysis to the beginning of the crisis (early September) by using just the 2 original analog stations and a cross-correlation coefficient of 0.8 to verify that observations made in the study period are representative of the beginning of the crisis.

Our analysis resulted in a catalog of 6,508 earthquakes from 21 October 2017 to 15 February 2018 (Table 1). The calculated daily counts correspond well with manual counts conducted by CVGHM during the crisis (Figure 2). Over the course of the entire study period, 43% of earthquakes in our catalog—including 9 of the 24 BMKG earthquakes recorded on three local stations or more (Supplementary Table S2) — group into one of 92 families of 5 earthquakes or more (Table 1). Several characteristics of repeating earthquake families—including the percentage of repeating earthquakes as a portion of total seismicity and the longevity of each family—change throughout the course of the unrest and eruption.

From the beginning of the study (21 October) through the onset of tremor on 12 November 2017, the pattern of repeating earthquakes was defined by a large number of simultaneous families, each with a large number of events (Figure 2; Table 1). In total, this time period included 4,065 earthquakes, 2,395 (∼59%) of which grouped into 65 different families (Table 1). These earthquakes ranged in magnitude from <M1 to M3, and all were brittle-failure or volcano-tectonic (VT) events.

Coinciding with the beginning of our study period, the BMKG catalog shows that earthquake locations began to occur N and NE of the volcano as well as on the W flank near Abang where a large majority of earthquakes had previously occurred (Syahbana et al., 2019). None of the BMKG-located earthquakes were part of an earthquake family until some of the 50 + located aftershocks that followed the M4.9 earthquake on 9 November 2017 05:54 SGT grouped into various earthquake families (Supplementary Table S2).

On 12 November, the behavior of earthquake families changed for a brief period. Most (54 of 65, or ∼83%) pre-existing families ceased on or shortly after 12 November, and new families appeared. The first new family (Family 66) appeared during the onset of tremor. Over the course of the next three days, a total of five new families appeared, comprising 35 earthquakes. The new families ranged in longevity from ∼1 h to ∼25 h. Waveforms were emergent and low-frequency (<5 Hz). None were large enough to be located with the BMKG network. Overall, from 12 November to 15 November, there were 183 total earthquakes, 69 (∼38%) of which were in earthquake families—including those that belonged to pre-existing families (Table 1).

The nature of earthquake families changed again several days prior to the onset of the eruptions. Starting on 15 November, earthquake families were defined by a small number that comprised a relatively low percentage of the overall seismicity but were long-lived. A majority of earthquakes during this time period were small (<M3) and were a mix of VT events and LF events. Overall, only 14% of earthquakes (311 out of 2260) were part of earthquake families (Table 1). Several families remained active through the end of the study. After 15 November, the BMKG catalog includes 13 earthquakes with sufficient data for our study, and 5 of them had repeating waveforms. All belong to families that appeared after 15 November (Supplementary Table S2).

Seismic rates remained low in the days prior to the initial phreatomagmatic eruption on 21 November, and there were no unique patterns among the few repeating earthquakes that occurred. The swarm of larger, low frequency earthquakes that roughly coincided with the onset of lava effusion on 25 November (see Eruption Timeline and Response), however, is clearly highlighted in our results. Starting at ∼0530h SGT, a new family of 21 highly repetitive (cross-correlation coefficient, CCC ∼0.95) earthquakes appeared (Supplementary Figure S1). These earthquakes had a strong peak energy at ∼2 Hz and are clearly identifiable by their frequency index in Figure 2B.

The remainder of the study period, which included continuous to semi-continuous ash explosions until 19 January 2018, was characterized by low rates of seismicity, low numbers of repeating earthquakes, and long-lasting families (Figure 2).

Extending the analysis to the beginning of the crisis using just the 2 analog, short-period stations resulted in a larger number of detected earthquakes. Many of the additional earthquakes were either lower amplitude with lower signal-to-noize ratio or were larger magnitude, clipped signals from the most intense part of the unrest. This resulted in a smaller proportion of repeating earthquakes overall (∼30%), but certain key observations remained consistent between this longer term analysis with lower quality data and our primary results—a large number of simultaneous families at the beginning of the crisis, the cessation of a most families (>90%) on or shortly after 12 November, the overall decrease in repeating earthquake percentage over time, a prominent swarm of 20 + LFs on 25 November, and notably longer-lived families after the onset of lava effusion.

Earthquake families—based on cross-correlation analyses—have long-been studied at a large number of volcanoes. However, very few examples document earthquake family evolution associated with the lateral emplacement of a dike, followed by magma migration and eruption (White et al., 2011). The evolution of earthquake families at Mount Agung illustrates a changing volcanic system. We interpret the transitions in seismicity as the manifestation of a three-phase physical model. The three phases inferred by our analysis in context of other observations are: 1) the Intrusion Phase, 2) the Transitional Phase, and 3) the Eruptive Phase (Figure 2; Figure 3).

Although our study was completed retrospectively, we can speculate on the role of our analysis in real-time eruption forecasting. In general, the operational value of repeating earthquakes depends on the ability to recognize patterns and interpret the processes that drive those patterns. This depends, in part, on the characteristics of the crust, the conduit, and the magma itself. It also depends on the availability of seismic records from past eruptions because volcanoes can display similar precursors to multiple eruptions.

Traditionally, seismologists rely on metrics such as earthquake rates, amplitudes, and frequency-based classifications to interpret unrest at a particular volcano. Repeating earthquake analyses provide additional metrics—such as the repeater percentage of total seismicity, the number of concurrent families, and the longevity of each family—that can aid the interpretation of unrest.

We are generally successful at detecting and alarming on the precursors we expect to see, while we are far less successful if we do not know what to expect from a volcano. Automatic detection of changes in repeating earthquakes could be an exceptionally valuable tool for eruption forecasting at well-understood volcanoes. For example, some studies suggest that certain sequences may be alarmable. Explosions at Volcán de Colima in 2004–2005 (Arámbula-Mendoza et al., 2011) and at Augustine in 2006 (Buurman and West, 2010) were preceded by diagnostic occurrences of certain earthquake families. In addition, cyclic patterns of dome deflation at Soufrière Hills in 1997 were preceded by the same pattern of repeating earthquakes (Green and Neuberg, 2006).

Repeating earthquakes have been most well-studied during episodes of dome growth, thus providing several models and common observations to explain their occurrence. For example, Thelen et al. (2011) note that earthquake families tended to be longer-lived during stable phases of effusion across multiple eruptions at Mount St. Helens and Bezymianny. This is consistent with our observation of relatively long-lived families after lava extruded at Agung.

Although repeating earthquakes have been most well-studied at domes, earthquake sources repeat throughout the crust beneath volcanoes of all types. Therefore, all phases of unrest that produce seismicity are likely to produce repeating seismicity. When we lack a clear paradigm for eruption precursors at a volcano, we interpret repeating earthquakes by applying generic models that incorporate knowledge of the magma properties, the stability of the conduit and degassing pathways, and the overall geomechanics of the crust. For example, increasing numbers of VT earthquake families at a closed-system volcano, particularly in combination with increasing seismic amplitudes, may be a concerning sign of pressurization of the edifice as more faults are activated. In contrast, increasing numbers or amplitudes of LP earthquakes within stable earthquake families at a volcano open to degassing, may simply mean that the degassing rate has increased without a significant increase in the overall hazard. Or during an intrusive phase of unrest, for example, rapidly changing sets of earthquake families may indicate migrating magma, while stable earthquake family may indicate rapid reloading of shear-stress on the same faults.

At Agung, little was known about crustal and conduit conditions or about seismicity associated with prior eruptions. However, clear changes in repeating earthquakes—such as those on 12 November and 15 November—could have given observatory scientists an indication that the volcanic system was evolving when other monitoring parameters were stable. This could have increased vigilance and encouraged more cross disciplinary discussion of the volcano’s current state. During a protracted unrest sequence, such as this one, scientist fatigue can be a significant problem, and any clear indication of meaningful changes in monitoring parameters is an advantage.

Seismicity is often the most readily available data stream for tracking unrest during volcanic crizes. Interpreting changes in earthquake rates, earthquake magnitudes, earthquake classifications, and RSAM has led to many successful forecasts. When unrest does not progress in a traditional pattern, however, relying on these metrics leads to uncertainty. In these scenarios, other tools add value. We analyzed the evolution of earthquake families with time during the Mount Agung volcanic crisis. We find that when considered in context of other observations, the time history of earthquake families provides insight into the evolution of the stress distribution in the volcanic edifice, the development of the volcanic conduit, and seismogenesis of magma effusion.

CVGHM staff noted the new appearance of volcanic tremor—nominally a harbinger for an evolving magmatic system—on 12 November 2017 after months of elevated seismicity at Mount Agung, but volcanic tremor and LF seismicity were subtle and earthquake rates were not high compared to prior months. Thus, the time frame over which to expect eruption remained uncertain.

Our retrospective study shows that earthquake families that dominated seismicity during the early stages of unrest ceased near the onset of tremor, highlighting that stress in the crust had been re-distributed as magma migrated toward the summit. Furthermore, earthquake families at Mount Agung took on characteristics commonly observed during effusive phases of eruptions on 15 November—a full six days before the first phreatomagmatic eruption on 21 November 2017 and a full ten days prior to the first magmatic eruption on 25 November 2017.

It is feasible to conduct advanced analyses like this in near-real time thanks to the availability of high-quality open-source codes written by the seismology community in recent years. This study demonstrates how analyzing earthquake families can be used to improve future eruption-response efforts.