Introduction

Kīlauea volcano is the youngest and most active of the five volcanoes that form the Island of Hawai‘i. Currently, the volcano consists of a summit caldera overlying a magma reservoir complex which includes at least two storage zones, a shallow zone at 1–2 km and a deep zone between 3- and 5-km depth (Poland et al. 2014). Extending from the summit are two rift zones—the East Rift Zone (ERZ) and the Southwest Rift Zone (SWRZ) (Fig. 1). Prior to May 2018, Kīlauea had been erupting semi-continuously for the past 35 years, with activity focused at two vents: a middle East Rift Zone (MERZ) vent, Pu‘u ‘Ō‘ō, located 20 km down-rift from the summit and active since 1983 (Heliker and Brantley 2002); and a summit lava lake within the Halemaʻumaʻu pit crater, active since 2008 (Patrick et al. 2018).

Fig. 1
figure 1

Seismic ambient noise interferometry network map of Kīlauea, Hawai‘i. (Main) Three-component (triangle) and single-component (circle) seismometers used in the Southwest Rift Zone (yellow), summit (blues/red), and East Rift Zone (green) ambient noise interferometry (ANI) groups. The approximate extent of the two rift zones are given by the dashed lines. An InSAR interferogram spanning the final month of the eruption (July–August 2018) is overlain in the Southwest Rift Zone. Also shown are the areas of erupted lava in the Lower East Rift Zone (red) and the location and focal mechanisms for two earthquakes (focal spheres), the pre-eruptive 2017 Ml5.3, and the syn-eruptive 2018 Mw6.9. (Inset right) The summit seismic network and location of the summit lava lake (dark red) within the Halema‘uma‘u pit crater. Summit stations are colored by subgroupings (north/south); see Fig. 6 for details. Sublabels under summit stations UWE and RIMD are GPS (blue) and tiltmeter (green) station codes

Full size image

On April 30, 2018, the Pu‘u ‘Ō‘ō eruptive phase dramatically ended when the vent collapsed and magma propagated down-rift into the lower East Rift Zone (LERZ). The collapse of Pu‘u ‘Ō‘ō was followed by the opening of 24 eruptive fissures in the LERZ, eventually extending ~ 6.8 km (Neal et al. 2018). The 3-month-long LERZ eruption was one of the volcano’s most voluminous eruptions in 500 years, with lava eruption rates exceeding 100 m3/s and lava covering ~ 36 km2 (Neal et al. 2018). At the summit, the draining of summit-stored magma led to near-daily collapses of a portion of the caldera and ultimately up to 500 m of summit subsidence. While deformation data indicated that the summit and middle East Rift Zone were inflating for the previous several years (Fig. 3), why Pu‘u ‘Ō‘ō collapsed and what initiated down-rift dike propagation remains unclear.

Seismic ambient noise interferometry (ANI) is a developing method for time-lapse monitoring of changes in relative seismic velocity in the shallow crust (Lecocq et al. 2014). The method expands on the well-proven technique to reconstruct empirical Green’s functions by the cross-correlation of ambient noise waveforms from seismic station pairs (Campillo and Paul 2003). These empirical Green’s functions, also referred to as noise cross-correlation functions (NCFs), are a proxy for the Earth’s response between the stations (Weaver and Lobkis 2002; Campillo and Paul 2003). ANI relies on time-progressive comparisons of NCFs made to a reference function to quantify changes in relative seismic velocity, and thereby changes in the Earth’s response, likely from changes in material properties, intensive parameters (pressure, temperature, etc.), and/or crustal structure. These changes are predominantly observed on NCFs calculated from the vertical components of seismic instrument station pairs (ZZ), indicative of changes in Rayleigh surface-wave velocity. However, changes can also be observed on horizontal components, e.g., radial-radial (RR), the subordinate Rayleigh-wave polarization, and transverse-transverse (TT), from changes in Love-wave velocity.

While the surface-wave phases of time-progressive NCFs can be compared directly to reference NCFs, ANI analyses are usually constrained to the coda-wave portions of the NCFs, where ideally multi-scattered Rayleigh and Love waves dominate the signal. These multi-scattered waves have path lengths longer than the direct arrival phases and are thereby more sensitive to changes in the medium and less affected by variations in the noise source (Grêt et al. 2005; Obermann et al. 2013). Under ideal multi-scattering conditions, without surface-to-body wave conversions (Obermann et al. 2013), the later-arriving portion of the coda-wave NCFs should not preserve along station-station pair directionality. Observed changes in seismic velocity should represent isotropic averages. However, in the early-arriving portion of the coda wave direct Rayleigh and Love wave phases may still dominate TT and RR NCFs, respectively. While these direct surface wave contributions decrease with increasing time into the well-scattered (later) coda, studies of where this transition occurs are limited (Obermann et al. 2013; Takano et al. 2019). With these caveats, anisotropic responses may be observed and interpreted for TT and RR NCFs; however, caution should be exercised when extending these interpretations to quantitative analyses.

At volcanic settings, changes in relative seismic velocity include those from seasonal hydrologic variations (Sens-Schönfelder and Wegler 2006; Haney et al. 2014; Hotovec-Ellis et al. 2014; Wang et al. 2017); magma pressurization changes, which may be accompanied with associated crack dilation or compression (Brenguier et al. 2008a; Duputel et al. 2009; Nagaoka et al. 2010; Anggono et al. 2012; Obermann et al. 2013; Donaldson et al. 2017; Donaldson et al. 2019); structural changes in the volcanic edifice due to dome or caldera collapse (Baptie 2010; Hotovec-Ellis et al. 2015); changes in the hydrothermal system (Baptie 2010; Mordret et al. 2010; Caudron et al. 2015; Taira and Brenguier 2016); and changes in poroelastic contraction from geothermal production and fluid withdrawal (Taira et al. 2018). Emplacement or removal of magma into pre-existing void space and changes in the physical properties of magma as it cools, e.g., density, bulk, and shear modulus, may also affect seismic velocities.

Here, we use results from both a real-time ANI seismic-velocity monitoring network—implemented by the US Geological Survey within days following the onset of the first LERZ fissure eruption—and two independent larger retrospective analyses, to investigate precursory conditions that may have contributed to down-rift dike propagation and describe relative seismic velocity changes observed throughout the eruption (Fig. 1).

Summary of eruption precursors and eruption onset

Starting in March 2018, three swarms of lower-crustal earthquakes were observed beneath the summit of Kīlauea (Fig. 2 and Movie S1). The first swarm began on March 7, lasted for approximately 1 day, and was focused between 6- and 12-km depth (~ 55 earthquakes, M1-2.3). The second and largest swarm began on April 11 with a similar depth range and duration, but included numerous shallow, < 6-km deep, earthquakes (~ 170 earthquakes, M1-2.3). A third, smaller swarm was observed on April 16 (~ 32 earthquakes, M1-2.2; Fig. 2). While similar swarms were common prior to 2005, the March–April 2018 swarms were three of only seven similar swarms to occur in the previous 10 years.

Fig. 2
figure 2

Kīlauea Caldera and upper East Rift Zone pre-eruptive seismicity for 2 months before the May 3rd, 2018 eruption. (Main) Pre-eruptive seismicity from April 1 to May 2, 2018. Only earthquakes shallower than 15-km depth with magnitudes greater than M1.0 are shown. Earthquake sizes are scaled by magnitude and colored by depth below sea level. The April shallow-crustal earthquake swarm in the connection between the summit and upper East Rift Zone, synchronous with the Halema‘uma‘u pit-crater lava-lake overflows, is highlighted in turquoise. The approximate extent of the two rift zones is given by the dashed lines. (Inset top right) The same pre-eruptive seismicity shown in the main panel but in depth through time. Highlighted regions denote the three pre-eruptive lower crustal earthquake swarms (purple), the April shallow-crustal swarm (turquoise), and the collapse of Pu‘u ‘Ō‘ō (orange). (Inset bottom right) The same extent shown in Fig. 1

Full size image

(MP4 102,940 kb).

Following the March 7 swarm, beginning on approximately March 13, there was increasing positive radial tilt (~ 0.2 μrad/day) at station UWD (colocated with seismic station UWE; Figs. 1, 2 and 3) and an increase in the distance between the across-caldera GPS station-pair CRIM-UWEV (~ 0.5 cm/week; colocated with seismic stations RIMD and UWE, respectively). These indicated increasing rates of inflation at Kīlauea’s summit. Similar increases were seen at Pu‘u ‘Ō‘ō (tiltmeter POC, ~ 0.7 μrad/day; GPS stations NPOC-JCUZ, 0.2 cm/week; Figs. 1 and 3). Inflation was accompanied by high lava-lake levels at both locations. At the summit, high lava-lake levels between April 21–24 produced only the third episode of lava flowing over the crater rim in the past decade—with each successive overflow representing a new high (Patrick et al. 2018). Continuing through the April overflows, there was an additional swarm of shallow-crustal earthquakes in the connection between the caldera and ERZ (Fig. 2 and Movie S1). This shallow seismicity was typical of previously observed summit inflation and high lava-lake levels and likely reflected increased pressurization in the upper ERZ, also referred to as the East Rift connector (Patrick et al. 2018; Swanson et al. 2019).

Fig. 3
figure 3

Two years of relative seismic velocity variations (dv/v), radial tilt, and GPS. (Top three) Representative station-pair records for each of the three ANI groups (black) padded by their uncertainty (orange). Colored areas (blue green, yellow) are the one-sigma ranges for the stack of all station-pair records in that subgroup. Highlighted windows qualitatively outline portions of the records possibly associated with (1) earthquake induced dv/v changes (gray) or (2) the 2018 eruption time range (red). All dv/v records are for the ZZ component correlations in the 0.4–1.0-Hz frequency band using a 5-day running stack and a ten-day running boxcar smoothing filter. (Bottom) Pre-eruptive mean radial tilt changes for Pu‘u ‘Ō‘ō (light green; POC, azimuth = 308°) and the summit (light blue; UWD, azimuth = 327°), and changes in GPS line lengths for a Pu‘u ‘Ō‘ō (dark green; NPOC-JCUZ) and a summit (dark blue; CRIM-UWEV) GPS station pair. The UWD tilt and UWEV GPS station are collocated with the UWE seismometer (Fig. 1). Tilt station POC is located ~ 200 m SSE of GPS/seismic station NPOC (Fig. 1)

Full size image

On April 27–28, Pu‘u ‘Ō‘ō was filled nearly to its crater rim and high rates of inflation of the vent continued. On April 30, at ~ 2:15 pm Hawaii Standard Time (HST), the Pu‘u ‘Ō‘ō vent floor collapsed. This collapse was immediately followed by a down-rift propagating increase in seismicity and deformation along a large section of the ERZ as well as a gradual deflation at the summit. The first small ground cracks in and adjacent to the Leilani Estates neighborhood in LERZ, ~ 23 km away from Pu‘u ‘Ō‘ō, were observed on May 2. The following day (May 3), at 10:30 am HST, a Mb5.1 earthquake caused additional collapse of Pu‘u ‘Ō‘ō. Approximately 6.5 h later, at ~ 5:00 pm HST, an eruption started in the LERZ, with fissures opening in the Leilani Estates neighborhood (Neal et al. 2018). At 12:32 pm HST on May 4, a Mw6.9 earthquake occurred on the decollement beneath the south flank of the volcano—the largest earthquake to occur on the island in 43 years.

For the 2 years prior to the 2018 lower East Rift Zone eruption, deformation data indicated that the summit and middle East Rift Zone magma plumbing systems were likely inflating with magma (Fig. 3). However, the cause of the rapid increase in inflation in the 2-month period preceding this new eruption remains unclear. Possible causes include (1) an increase in magma supply into Kīlauea’s magma reservoirs, as suggested by the three pre-eruptive lower-crustal summit-earthquake swarms; (2) a decrease in the eruptive efficiency and eruption rate from Pu‘u ‘Ō‘ō, which in turn promoted up-rift magma accumulation—similar to processes observed in 1992 and 2011 (Kauahikaua et al. 1996; Orr et al. 2015; Patrick et al. 2019); or (3) some combination of both. Equally important, but unknown, are the mechanisms that promoted down-rift dike propagation into the LERZ, opposed to continued eruption at the 35-year-lived Pu‘u ‘Ō‘ō vent (Orr et al. 2015). Here, we use ambient noise seismic interferometry to perform a retrospective analysis focused on addressing these questions.

Methods

Source data

Within days following the start of the first LERZ fissure eruption, we implemented a pseudo-real-time seismic velocity monitoring network using ambient noise interferometry (ANI) (Lecocq et al. 2014). The motivation was to determine whether seismic velocity changes in the LERZ were indicative of propagating intrusions and potentially additional or continued fissure eruptions. At the time of the eruption, the seismic network in the LERZ was limited to three stations: two stations on opposite ends of the Leilani Estates neighborhood, separated by ~ 4 km (PHOD and KLUD), and one station 11 km further up-rift (JOKA) (Fig. 1). Only the up-rift station (JOKA) is a three-component broadband sensor (Nanometrics Trillium 120). Station PHOD, located near the Puna Geothermal Venture, northeast of the first fissure eruptions (fissures 1–15), was a three-component accelerometer (Kinemetrics Episensor, 4 g setting). Station KLUD, on the SE border of the subdivision, is a vertical short-period geophone (Mark Products L-4). To compare to synchronous changes further up-rift (~ 22 km), a fourth station was also included (JCUZ, Guralp CMG-40T).

By late May, a second ANI network was implemented for Kīlauea’s summit (15 stations, Fig. 1), and a third for the SWRZ (5 stations; Fig. 1). For retrospective analysis, following the cessation of eruptive activity, the LERZ network was expanded to include an additional five upper and middle ERZ stations. All stations are maintained by the US Geological Survey’s Hawaiian Volcano Observatory and their data were available in real time from the Incorporated Research Institutions for Seismology (IRIS). We should note, the LERZ station PHOD was destroyed by a lava-induced forest fire on August 1, ~ 90 days after the start of the eruption. Additionally, during the eruption, a large crack opened directly beneath station KLUD. Three summit stations were also temporarily or permanently lost between May 25 and June 1 from caldera collapse (NPT) and power fluctuations (NAHU and PUHI). Although not used in this analysis, six additional seismic stations were established in the LERZ during the eruption (ERZ1-4, P07, 2816 and KIND), and were critical in real-time hazard analysis and earthquake relocation.

Relative seismic velocity change (dv/v) measurement

Noise cross-correlation functions (NCFs) were computed using two independent processing workflows, both implemented using the MSNoise software package (Lecocq et al. 2014). In the first workflow, NCFs were computed for ~ 3 years of continuous data prior to May 2018 using a data-processing scheme similar to Brenguier et al. (2008b), Taira et al. (2015) and Taira and Brenguier (2016). (1) First, instrument responses were removed from 1-day-long records, band-passed between 0.025 and 2 Hz, down-sampled to 20 Hz, split into 30-min-long segments clipped to 3 rms, and spectrally whitened in the frequency range of 0.0666–1 Hz (i.e., 1–15 s). (2) NCFs were computed for these 30-min-long segments for all possible station pairs and varying channel combinations, for a lapse time between ± 120 s. For the ERZ and SWRZ, all channel combinations were computed (i.e., ZZ, TT, RR, ZR, ZT, TR), while for the summit only self-similar combinations were calculated (ZZ, TT, RR). Daily NCFs were obtained by stacking the 30-min NCFs, which were then subsequently stacked again with a sliding window of 3, 5, 10, and 30 days. (3) Delay times were estimated between these windowed NCFs and reference NCFs, using a band-limited moving window cross-spectral (MWCS) approach (Clarke et al. 2011). Reference NCFs were defined as the 1-year average station-station NCFs from May 1, 2015, to May 1, 2016. Five MWCS frequency bands, spanning from 0.07266 to 1.0 Hz and overlapping at the approximate − 1.8-dB level (i.e., two-thirds power level), were investigated (Flinders and Shen 2017; Flinders et al. 2018). Delay time estimates were limited to the coda portions of the NCFs, defined as a window between twice the NCF-dependent minimum interstation shear-wave travel time (~ 1.5 km/s at 500 m) and the subsequent 30 s (Saccorotti et al. 2003). (4) Relative velocity changes (e.g., dv/v) were then calculated as the slope of the MWCS delay times versus time lags (e.g., dt/t). Only MWCS measurements that had a minimum coherence of 0.8 were used in the dv/v estimate.

The second processing workflow was similar to the first, but focused on vertical-component (ZZ) NCFs computed since 2017, using a different initial bandpass filter (0.01–8.0 Hz) and spectral-whitening frequency range (0.3–1 Hz). Unlike the 1-year reference functions used in the first workflow, reference NCFs were defined as the 6-month average station-station NCFs from January 1, 2017, to June 1, 2017, and a single MWCS frequency band was used (0.3–1 Hz). Additionally, constraints used to compute dt/t were varied following the procedure described in Donaldson et al. (2017); only points with an error of less than 0.1 s, a coherence of more than 0.65, and a dt of less than 0.1 s were accepted. Again, delay time estimates were limited to twice the NCF-dependent minimum interstation shear-wave travel time, which was varied from 1.5 to 2.0 km/s. The effect of this velocity produced negligible differences.

Results

Southwest rift zone

At Kīlauea’s Southwest Rift Zone (SWRZ), we observe a noticeable decrease in relative seismic velocity (dv/v) associated with a Ml5.3 earthquake on the decollement underlying the volcano’s south flank on June 8, 2017 (Fig. 3). This decrease is most pronounced in a high-frequency band (0.4–1.0 Hz) equivalent to a Rayleigh surface-wave (ZZ) depth of ~ 0.5–1.5 km (Fig. S1). This ~ 0.5% relative velocity decrease is observed on all SWRZ ZZ station-pairs and all component combinations for the across-rift pair DESD-HLPD (Fig. 3).

During the 2018 LERZ eruption, dv/v changes in the SWRZ are again the most notable in the high-frequency band of DESD-HLPD and are visible to as low as 0.25–0.43 Hz, or 1.5–4-km depth for ZZ correlations (Rayleigh) and 0.25–1.75-km depth for TT (Love; Fig. 4 and S2). Despite the high eruption-related seismic noise, NCFs for all three component combinations remained largely coherent relative to the pre-eruption time period (Fig. S2). In the high-frequency band (0.4–1.0 Hz), the dv/v component combinations ZZ and RR decrease (~ 0.4–0.75%) several weeks before the onset of the eruption, with the decrease accelerating in response to the MW6.9 earthquake on May 4 (Neal et al. 2018; Figs. 3 and 4, Movie S1). The TT combination increases for ~ 10 days following the MW6.9 then subsequently decreases (Fig. 3). Decreases on all component combinations last throughout May and could be attributed to coseismic microcrack opening. The initial opposite behavior of the TT-correlations may be explained by short-lived shallow crack closing in the across-rift direction (TT) and accompanying preferential opening in deeper vertical along-rift (RR) and horizontal (ZZ) planes. This same behavior is observed for the three component combinations for the across-rift Pu‘u ‘Ō‘ō station-pair JCUZ-NPOC (Fig. 5) and is consistent with strain patterns revealed by interferometric synthetic aperture radar (InSAR) and the direction of earthquake slip (Neal et al. 2018).

Fig. 4
figure 4

Eruption-related relative seismic velocity variations (dv/v) in the Southwest Rift Zone. Relative seismic velocity variations for the station pair DESD-HLPD from 3 months preceding the eruption until 2 months after. Each component correlation—ZZ (red), TT (blue), and RR (gray)—is padded by its uncertainty. DESD and HLPD are the only three-component seismometers in the SWRZ. Highlighted windows qualitatively outline portions of the record possibly associated with (1) earthquake-induced dv/v changes (green) and (2) possible anisotropic effects (orange) such as along-rift crack opening. Also shown are the May 4 MW6.9 earthquake (focal sphere) and dates of noted eruption-related events. All records are for the 0.4–1.0-Hz band using a 5-day running stack and a 10-day running boxcar smoothing filter

Full size image
Fig. 5
figure 5

Eruption-related relative seismic velocity variations (dv/v) at Pu‘u ‘Ō‘ō. Relative seismic velocity variations for the station pair NPOC-JCUZ from 3 months preceding the eruption until 2 months after. Each component correlation—ZZ (red), TT (blue), and RR (gray)—is padded by its uncertainty. Highlighted windows qualitatively outline portions of the record possibly associated with (1) earthquake-induced dv/v changes (green) and (2) possible anisotropic effects (orange) such as along-rift crack opening. Also shown are the May 4 MW6.9 earthquake (focal sphere) and dates of noted eruption-related events. All records are for the 0.4–1.0-Hz band using a 5-day running stack and a 10-day running boxcar smoothing filter

Full size image

One week after the onset of large, near-daily, caldera collapses on May 29, and accelerated summit subsidence, all three dv/v component combinations rapidly increase (Fig. 4). This increase begins in the first week of June, lasts for approximately 2 weeks, and coincidentally for the RR and TT correlations which are more sensitive to shallow changes, and ends days before the caldera collapse expanded to include a larger portion of the summit (~June 24). The cause of this increase is unclear but a possible mechanism could be poroelastic contraction from a decrease in up-rift summit magma-static pressure and consequential withdrawal of SWRZ magmatic fluids (void closing). However, this same behavior is also observed for the MERZ station pair at Pu‘u ‘Ō‘ō, JCUZ-NPOC (Fig. 4).

Throughout the rest of the eruption (July–August), dv/v for the TT and ZZ correlations of DESD-HLPD remain elevated. However, dv/v for the RR correlations (Fig. 4) was decreased, mirroring in polarity the changes seen on TT/ZZ. This anisotropic behavior could be due to differences in strain in the along-rift and across-rift directions. An interferogram spanning July 7 to August 7, 2018 (Fig. 1), suggests subsidence in an elongated pattern, with a north-east major axis parallel to the SWRZ, supporting an interpretation of an anisotropic strain field (data acquired by Sentinel-1, European Space Agency). While it is unclear whether any macroscopic faulting occurred between DESD and HLPD, changes could be related to microscopic strain changes on vertically oriented cooling joints in previous emplaced lava flows (Duffield 1975). Coincidentally, this behavior ceases 7–10 days after the last caldera collapse event (August 2).

Kīlauea Caldera and the east rift zone

The decrease in dv/v associated with the June 2017 earthquake in the 0.4–1.0 Hz band is also observed throughout the East Rift Zone and possibly at the summit (Fig. 5 and 6). Additionally, prior to the May 2018 eruption, beginning in early March, a small amplitude gradual increase in dv/v is observed (Figs. 5 and 6), coincident with the first summit earthquake swarm and several days prior to increased rates of summit inflation (UWD; Fig. 3). This change culminated on ~April 20–25, then reversed, decreasing until May ~ 1, immediately preceding the collapse of Pu‘u ‘Ō‘ō. Following the first week of the eruption, the northern and southern sections of the summit changed inconsistently (Fig. 6). In workflow 1, a significantly larger amplitude and increasing change in dv/v to the north compared to the south is observed. In workflow 2, the opposite is observed. These inconsistencies may be due to differing process parameters, e.g., minimum coherency limits, lower bound frequencies, or stations used. Additionally, the consistency of the NCFs during the eruption appears to be highly path-dependent with some station pairs more affected by eruption-related seismicity (Fig. S4).

Fig. 6
figure 6

Eruption-related relative seismic velocity variations (dv/v) at the summit. Vertical-vertical (ZZ) relative seismic velocity variations at the summit from 3 months preceding the eruption until the 2 months after. Measurements are subdivided into two caldera subgroups (north/south). For each subgroup (north/south), we measure velocity changes using two independent workflows (see “Methods”). Workflow 1 (red/blue) stacked ten ZZ station-pair records between five stations for both the north and south subgroups (see Fig. 1 for station locations). Colored areas (light blue/red) are the one-sigma ranges for each subgroup stack. The dashed lines (August/September) indicate high signal decorrelation (unreliable). Workflow 2 (green/violet) used the same stations as workflow 1, excluding SBL in the north, and RIMD in the south. Each subgroup in workflow 2 was a stack of six ZZ station-pair records. The one sigma ranges for workflow two are given by vertical error bars. Also shown are the approximate date ranges of the summit-swarm coincident pre-eruptive dv/v increase (blue), the pre-eruptive dv/v decrease (red), and dates of noted eruption-related events. All records are for the 0.4–1.0-Hz band using a 5-day (workflow 1) or 1-day (workflow 2) running stack and a 10-day running boxcar smoothing filter

Full size image

By the end of May and the onset of near-daily caldera collapses, inconsistencies in the summit responses ceased, and dv/v throughout the summit decreases. This consistent decrease is short-lived, lasting for ~ 10 days, after which an increase in dv/v continues throughout June, culminating on ~June 24 when the caldera collapses expanded to include a larger portion of the summit (Fig. 6). This cycling between decreased and increased dv/v could reflect progressive pressurization/depressurization of the shallow summit reservoir. Between June 25 and August 2, coincident with the end of caldera collapse events, dv/v decreases in the north of the summit and increases in the south. Later, dv/v estimates become too scattered for meaningful interpretation.

Following the first fissure outbreak and the collapse of Pu‘u ‘Ō‘ō, we also observe decreases in dv/v in the upper and middle ERZ (Figs. 5 and 7). The decrease is largest at Pu‘u ‘Ō‘ō itself, with the next largest changes occurring up-rift (upper ERZ), and no detectable change in the middle ERZ. At Pu‘u ‘Ō‘ō, ZZ and RR dv/v remain decreased until the onset of daily caldera collapses, when dv/v begins to increase rapidly (Fig. 5). A corresponding slight increase in dv/v might also be interpreted in the UERZ (Fig. 7). Increased dv/v at Pu‘u ‘Ō‘ō reached a maximum on ~ July 15, before decreasing over the next 2 weeks (Figs. 5 and 7). An additional cycle of increase/decrease occurs again in early August, culminating close to when caldera collapse events end. Due to the lack of pre-eruption station density in the LERZ, station quality, and high-rates of seismicity, dv/v changes throughout the LERZ were highly erratic and largely indistinguishable from the long-term average background variation.

Fig. 7
figure 7

Eruption-related relative seismic velocity variations (dv/v) in the upper/middle East Rift Zone. Vertical-vertical (ZZ) relative seismic velocity variations for three East Rift Zone subgroups from 3 months preceding the eruption until the 2 months after. The upper East Rift Zone (UERZ) subgroup (orange) is a stack of nine ZZ station-pair records from five stations up-rift of Pu‘u ‘Ō‘ō. The middle East Rift Zone (MERZ) subgroup (green) is a stack of four ZZ station-pair records from stations between the Pu‘u ‘Ō‘ō and JOKA (Fig. 1). Both the UERZ and MERZ subgroups include stations NPOC and JCUZ but exclude the NPOC-JCUZ pair. The Pu‘u ‘Ō‘ō subgroup is the NPOC-JCUZ station-pair record. Colored areas (light red/green/orange) are the one-sigma ranges for each subgroup stack. Also shown are dates of noted eruption-related events. All records are for the 0.4–1.0-Hz band using a 5-day running stack and a 10-day running boxcar smoothing filter

Full size image

Discussion

Pre-eruptive summit inflation

At the summit, the pre-eruptive increase in dv/v between March and April 2018 correlates with modest increased rates of summit inflation (UWD, UWE-CRIM; Fig. 3). This increase was also observed by Olivier et al. (2019) and is consistent with the model proposed by Donaldson et al. (2017) in which increasing pressure in the 1–2-km depth shallow summit magma reservoir results in localized compression and closing of cracks in the near-surface (Donaldson et al. 2017; Olivier et al. 2019). Olivier et al. (2019) also showed a separate pre-eruptive decrease in dv/v 10 days preceding the first fissure outbreak, which we also observe at both the summit and the SWRZ (Figs. 4 and 6). The decrease could result from inflation of the deeper summit magma reservoir and wider spread summit dilation, or as proposed and preferred by Olivier et al. (2019) accumulating damage induced by magma reservoir pressure exerted on the surrounding edifice.

At the summit, the pre-eruptive dv/v increase begins following the first lower-crustal (6–12 km) earthquake swarm on March 7 (Figs. 2 and 6). These swarms could be indicative of short-term increases in magma supply (Wright and Klein 2006; Dzurisin and Poland 2019) which in turn increase inflation at the summit and the MERZ near Pu‘u ‘Ō‘ō. The correlation between the increase in dv/v and the first summit swarm implies rapid pressure increases in the shallow summit reservoir following an increase in magma supply into the base of the volcano’s deeper reservoir (3–5 km). The six-day delay between the first earthquake in the first lower-crustal summit swarm and the beginning of increased MERZ inflation equates to a transport speed of ~ 0.05 m/s. Assuming the volcano’s summit reservoir systems are hydraulically connected to the MERZ (Montagna and Gonnermann 2013), this transport speed is significantly slower than the propagation of magma-static pressure, implying the physical movement of magma down-rift.

Alternatively, the rapid increase in the summit and MERZ inflation rate could have resulted from a constriction between Pu‘u ‘Ō‘ō‘s shallow magma storage system and its actively erupting outflow vent. Again, assuming the MERZ and summit plumbing systems are hydraulically linked, both Pu‘u ‘Ō‘ō‘s eruptive vent efficiency and the vent’s elevation would influence summit pressure (Patrick et al. 2019). Patrick et al. (2019) proposed that large differences between the elevation of the summit and the Pu‘u ‘Ō‘ō lava lake, and/or between the Pu‘u ‘Ō‘ō lava lake and its outflow vent, can imply localized decreases in ERZ vent efficiency. A constriction between Pu‘u ‘Ō‘ō’s shallow magma storage system and its actively erupting outflow vent could cause up-rift magma accumulation and inflation, and decreased effusion rates (Patrick et al. 2019). This style of rift zone–driven summit inflation has likely occurred previously (Kauahikaua et al. 1996; Orr et al. 2015).

While Pu‘u ‘Ō‘ō‘s then erupting flow field (episode 61 g) had been waning since late 2017, any relationship to the three lower-crustal earthquake swarms (5–12 km) is ambiguous. If a constriction near Pu‘u ‘Ō‘ō had been driving pressurization in the summit and MERZ, inflation would likely precede the first earthquake swarm, with the first swarm occurring in the upper ERZ (East Rift connector). Similarly, an overall deepening in the average swarm location might also be expected. However, the opposite is observed (Movie S1 and Fig. 2). The deepest swarm on March 7, preceded the change in summit and Pu‘u ‘Ō‘ō vent inflation by approximately 6 days. Additionally, while the first summit earthquake swarm was quickly followed by an increase in summit dv/v, there was no corresponding change in dv/v near Pu‘u ‘Ō‘ō, despite increased inflation (Fig. 2). For these reasons, we favor a “bottom-up” mechanism whereby short-term increases of magma supply into Kīlauea’s deep reservoir caused increased inflation at the summit and the MERZ near Pu‘u ‘Ō‘ō. This increase in pressure ultimately leads to failure of the rift zone near Pu‘u ‘Ō‘ō and the initiation of the eruption.

Promoting down-rift intrusion

Almost a year prior to the 2018 LERZ eruption, the June 2017 Ml5.3 earthquake on the decollement underlying Kīlauea’s south flank caused a coseismic dv/v decrease in both the ERZ and SWRZ (Fig. 3, Table 1). For the near Pu‘u ‘Ō‘ō station-pair JCUZ-NPOC, we estimate pre- and post-earthquake 2-month dv/v averages of 0.08 ± 0.14% and 0.38 ± 0.16%, respectively, with the Ml5.3 earthquake inducing a 0.30 ± 0.21% velocity reduction. Ten months later, the pre-eruption 2-month average dv/v was 0.17 ± 0.20%, approximately half its pre-earthquake magnitude (Fig. 3). Despite considerable uncertainty, there appears to be a long-trend post-seismic recovery of ~ 0.03%/month. Stress-induced opening and closing of crack space is a commonly evoked mechanism for coseismic velocity changes (Adams and Williamson 1923; Lockner et al. 1977; Brenguier et al. 2008b, 2014; Lesage et al. 2014; Taira and Brenguier 2016). We attribute the observed coseismic dv/v decrease to increases in crack opening throughout both rift zones. At Kīlauea, this opening and closing could predominantly occur on vertically oriented cooling joints in previously emplaced lava flows (Duffield 1975). While the transient changes in dv/v occurred nearly instantaneously, in the ERZ, the long-term post-seismic recovery and healing of coseismic fracture may have never been fully completed (Fig. 3). Crustal damage caused by stresses from this earthquake could have increased the compliance of Kīlauea’s rift zones, potentially priming the ERZ for intrusion almost a year later.

Table 1 Stress estimate parameters and estimates
Full size table

We estimated volumetric static stress changes (Δσs) for the June 2017 earthquake from dilatation strains (EXX + EYY + EZZ) using the Coulomb 3.3 software (Toda et al. 2011) at five locations: two in the SWRZ (DESD, HLPD), one in the upper ERZ/ERZ connector (PAUD), at Pu‘u ‘Ō‘ō, and one in the LERZ (PHOD). In both Coulomb 3.3 and to convert the computed strain to stress, we used a shear modulus (G) of 25.6 GPa and a Poisson’s ratio of 0.15, based on the temperature dependence of Ryan (1987) assuming T = 140 °C, i.e., the temperature at 1-km depth in the Scientific Observation Hole 4 (SOH-4) borehole (Olson and Deymonaz 1993). The earthquake origin and focal mechanism were taken from the US Geological Survey Hawaiian Volcano Observatory (event ID: hv61792931) and National Earthquake Information Center (NEIC; event ID: us20009lc3) catalogs. The fault length and width were computed with an empirical relationship between earthquake magnitude and faulting style (Wells and Coppersmith 1994). The grid size was set to 200 m.

At these same five locations, we estimated the earthquake-induced dynamic stress changes (Δσd) using peak ground velocities (PGVs), where Δσd = (PGV/VS) G (Hill et al. 1993). For four of the five locations (excluding Pu‘u ‘Ō‘ō), PGV was taken as the average of the seismometer’s northing and easting component as calculated by the NEIC ShakeMap. At the Pu‘u ‘Ō‘ō, PGV was not available for the three closest stations (NPOC, JCUZ, STCD) and was estimated from NEIC ShakeMap’s interpolated PGV grid (Worden et al. 2018). We used the same shear modulus (25.6 GPa) used in calculating the volumetric static stress change, with a shear wave velocity (VS) of 2.3 km/s, i.e., the velocity at 1-km depth (Saccorotti et al. 2003; Lin et al. 2014).

The magnitude of dilatational static-stress changes for the June 8, 2017, earthquake, in both the SWRZ and ERZ, was less than 1 kPa, i.e., less than those from barometric pressure changes and tides. However, dynamic stress changes from this earthquake, inferred from PGV measurements, are estimated to be 0.08–0.76 MPa, with the largest changes occurring near Pu‘u ‘Ō‘ō (see Table 1). These dynamic stress changes are on the order of the “fracture toughness” for basalt, i.e., its resistance to dilatant crack growth (Kemeny and Cook 1988). Using the average flow thickness for Kīlauea of 1.1–1.86 m (Katz and Cashman 2003) as a proxy for crack length, the fracture toughness of basalt, ~ 0.9–2.5 MPa m−1/2, was almost certainly exceeded during the 2017 earthquake, particularly near Puʻu ʻŌʻō (0.91 ± 0.12 MPa m−1/2; Atkinson and Meredith 1987).

Similarly, we estimated changes in fracture toughness from dv/v, assuming they were only due to changes in shear modulus. Near Pu‘u ‘Ō‘ō, the June 2017 coseismic velocity decrease of 0.30 ± 0.21% is equivalent to a change in fracture toughness of ~ 0.47 ± 0.44 MPa m−1/2, comparable to the 0.91 MPa m−1/2 estimated from PGV. Ten months later, the pre-eruption 2-month average dv/v of 0.17 ± 0.20% corresponds to a remaining residual decrease in fracture toughness of ~ 0.09 ± 0.01 MPa m−1/2. This weakening of the MERZ’s resistance to dilatant crack growth combined with increased up-rift magma-static pressure, possibly from an increased summit supply, may have facilitated down-rift dike propagation and the devastating May 2018 eruption.

Conclusion

Earthquake-volcano interactions have long been proposed at Kīlauea (Lipman et al. 1985; Walter and Amelung 2006; Montgomery-Brown et al. 2015). The effects of the 1975 MW7.7 Kalapana earthquake on Kīlauea’s shallow magma plumbing system is thought to have initiated its middle East Rift Zone eruption almost 2 years later (Lipman et al. 1985). Similarly, at Mauna Loa, static stress changes associated with large flank earthquakes (e.g., > M6.6; 1868, 1949–1951) possibly induced subsequent eruptions through decompression of the shallow magma reservoir (Walter and Amelung 2006). Here, we present a new mechanism where long-lived crustal weakening resulting from a moderate magnitude local earthquake can later promote down drift-dike propagation and possibly an eruption. This mechanism is similar to that proposed for earthquake-triggered changes in geyser eruption frequency, where dynamic stress induces changes in permeability and fracture opening and closing (Hill et al. 2002; Manga and Brodsky 2006), and dynamic stress induced changes in volcanic heat flux (Delle Donne et al. 2010). We identify two necessary conditions of this mechanism: (1) the incomplete post-seismic recovery and healing of the coseismic fracture generated during the local earthquake, and (2) a synchronous increase in up-rift magma-static pressure. Dike propagation occurs once magma-static pressure exceeds the weakened volcano’s ability to resist dilatant crack growth. During the 2018 lower East Rift Zone eruption of Kīlauea, both of these pre-eruptive conditions were met. However, whether the increase in pressure up-rift of Pu‘u ‘Ō‘ō reflected an increase in magma supply remains uncertain.

While simultaneous visual and seismic monitoring of Kīlauea has been occurring since 1912, these ~ 100 years represent only a small time window of the eruptive history of Hawaiian volcanoes, and earthquake-modulated volcanism could be a common occurrence. Moderate magnitude local earthquakes could increase a volcano’s susceptibility to intrusion and/or eruption for months-to-years afterwards, i.e., until post-seismic healing is complete. Unknown yet are the effects of the syn-eruptive May 2018 MW6.9 earthquake on Mauna Loa’s and Kīlauea’s magma plumbing systems. However, both conditions required to generate dike propagation through the mechanism described here are observable in real time; increased pressurization through interpretations of tilt, GPS, and InSAR, and coseismic weakening through ambient noise interferometry. Ambient noise seismic interferometry may better enable us to characterize and identify this and other new precursory eruptive signals in near real time, proving invaluable for interpretations of ongoing unrest at these two volcanoes. While this technique has been well demonstrated scientifically, its application to real-time observatory response is still evolving and will rely on a robust seismic network and independent corroborating measurements (tilt, GPS, InSAR, etc.).