Detecting pre-eruptive magmatic processes of the 2018 eruption at Kilauea, Hawaii volcano with ambient noise interferometry

A sequence of eruptions in the lower East Rift Zone (LERZ) of the Kilauea, Hawaii volcano commenced in early May 2018 and caused serious damage and residential evacuation. The post-eruption drainage and collapses of the summit lava lake and caldera suggest a well-connected magmatic plumbing system along the rift zone. How and when the pre-eruptive magmatic processes in such connecting system occurred are, however, still unclear. For that, we apply ambient noise interferometry with seismic data from January 2017 to June 2018 from 12 broadband seismometers to investigate spatiotemporal seismic velocity changes (dv/v) of the upper crust in the Kilauea area. The dv/v variations in three frequency bands (0.3–0.6 Hz, 0.6–0.9 Hz, and 0.9–2.0 Hz) show distinct responses to strong earthquake ground shaking and deep magmatic intrusion processes. Earthquake-induced dv/v drops mainly in the higher two frequency bands imply shallow mechanical changes within the uppermost 1 km of the crust. In contrast, the magma-related dv/v changes can be characterized into three periods of activity: from the December 2017 to March 15, 2018, a dv/v excursion only seen in the lowest frequency band indicates the magmatic intrusion processes taking place at the depth range of 1–4 km, consistent with the proposed depth of the magma reservoir-dike system in the rift zone. The spatial dv/v distribution suggests that the magma may intrude to the deeper summit magma reservoir and in the upper East Rift Zone (UERZ) at this time. From March 15 to April 17 in 2018, the summit inflation recorded by tiltmeters causes the dv/v increases in the higher two frequency bands. After April 17 to the eruption, the accumulating damage of the edifice together with the stronger intrusion activity in the UERZ result in dv/v decreases in all three frequency bands around the summit and UERZ areas. Our observations highlight that the ambient noise interferometry analysis provides an opportunity to image and understand the pre-eruptive processes of reservoir-dike magma system and could be a useful supplement to current volcanic monitoring systems.


The lower East Rift Zone eruption of the Kilauea volcano
The 2018 eruption of the Kilauea volcano began in early May and created a series of new fissures in the lower East Rift Zone (LERZ, Fig. 1). The Hawaiian Volcano Observatory (HVO) issued notices of magma pressurization at Pu'u'Ō'ō on April 17 and lava overflowed from the Halemaʻumaʻu crater on April 26. Before the end of April, Pu'u'Ō'ō crater collapsed and a burst of seismicity occurred, migrating eastward from the middle East Rift Zone (MERZ) to the LERZ. On May 2, the first fissure opened near Leilani Estates followed by a series of new fissures opening in the LERZ. Until September 2018, lava flows have covered 35.5 km 2 of land and destroyed over 700 homes (Neal et al. 2019;Patrick et al. 2019). This protracted eruption drained magma from the summit reservoir and induced a series of caldera collapses ~ 40 km away, suggesting a well-connected magmatic plumbing Feng et al. Earth, Planets and Space (2020) 72:74 system along the East Rift Zone (Neal et al. 2019;Wu et al. 2020). However, when and how such connecting plumbing system formed, stored and transported magma from the summit to its lower flank prior to the eruption remains poorly understood.

Ambient noise interferometry (ANI) for monitoring volcanoes
The ambient noise interferometry (ANI) has become a promising tool for monitoring volcanic unrest and magma pressurization at depths in the past decade (Brenguier et al. 2008. The pioneer study (Brenguier et al. 2008) at volcano Piton de la Fournaise (PdF) found consistent decreases in seismic velocity a few weeks before eruptions over several cycles from 1999 to 2000 and suggested that precursory velocity decreases are caused by the inflation of volcanic edifice. These inflation-related velocity decreases prior to eruptions have been observed in later studies at PdF (Duputel et al. 2009;Obermann et al. 2013a) and other volcanoes (Mordret et al. 2010;Bennington et al. 2015Bennington et al. , 2018Machacca-Puma et al. 2019); However, velocity increases related to pre-eruptive inflation have also been found at volcanoes as Miyakejima volcano (Anggono et al. 2012), Merapi volcano (Budi-Santoso and Lesage 2016), and Kilauea volcano (Donaldson et al. 2017). Such contradiction can be reconciled with a model based on strain theory proposed by Donaldson et al. (2017). They showed whether the edifice undergoing extension or compression depends on the depth of inflating magma reservoir. The modeling of shallow magma reservoir at Kilauea (1-2 km below the surface) suggests that the most of the edifice is dominated by compression and hence results in increases of seismic velocity. Investigations on the variations in seismic velocity before eruptions can thus provide valuable insights into the strain state of the magmatic systems in response to the intrusion processes (Brenguier et al. 2008;Bennington et al. 2018).
For the May 2018 eruption of the Kilauea volcano, a recent study of Olivier et al. (2019) has observed seismic velocity increases and subsequent decreases in response to the inflation recorded by tiltmeters beginning from March and the accumulating damage of the volcano edifice 10 days before the eruption. These observations are, Fig. 1 Distribution of broadband stations used in the Kilauea area. The volcanoes, used seismic stations, and earthquakes with M ≥ 3.5 (from United States Geological Survey earthquake catalog) are shown by the red triangles, black inverted triangles, and circles color-coded with focal depth. Red dash line marks the fissures opened during the 2018 eruption. The blue diamond marks the location of the most active fissure vent 8. Two shaded pink zones indicate the East and Southwest Rift zones (referred to Poland et al. 2014). The inset map shows Hawaiian Island chain and the location of our study area (red box) however, more focusing on short-term precursory activity around the summit. To further understand the magma connectivity and interactions between the summit and LERZ, we conduct ANI analysis in broader space and time to monitor and detect possible mechanical or structural changes in seismic velocity (dv/v) associated with pre-eruptive magmatic intrusion processes. Using the vertical component of daily seismic recordings from 12 selected broadband seismic stations distributed around the summit and along the ERZ from January 2017 to June 2018 (Fig. 1), we measure dv/v in three distinct frequency bands of 0.3-0.6 Hz, 0.6-0.9 Hz, and 0.9-2.0 Hz. The frequency-dependent responses of dv/v clearly show different mechanisms related to earthquake ground shaking and magmatic intrusion processes in different depth ranges. Spatiotemporal dv/v patterns that mainly show velocity decreases along the upper East Rift Zone (UERZ) and increases around the summit and surrounding areas suggest that magmatic intrusion processes beginning from December 2017 until the May 2018 eruption reflect dike-reservoir inflation and edifice accumulating damage in three stages.

Background of Kilauea volcano
The Kilauea volcanic area consists of a main summit and two rift zones, the East Rift Zone (ERZ) and the Southwest Rift Zone (SWRZ). In the past 30 years, major eruptions at Kilauea have occurred along the ERZ between the summit and Pu'u'Ō'ō ( Fig. 1) and built up a series of cones along the linear fissures of the rift. The last eruption happened in 2008 followed by several months of increased sulfur dioxide emissions and seismicity and formed a lava lake within Halemaʻumaʻu crater at Kilauea's summit. The active and connected magmatic plumbing system beneath the Kilauea volcano includes a shallow magma reservoir at ~ 1 km beneath the Halemaʻumaʻu crater (Halemaʻumaʻu reservoir), a deeper magma reservoir at a depth of ~ 3 km below the southern side of the Kilauea caldera (Southern Caldera reservoir), and two dike systems extending from the summit to the SWRZ and ERZ (Heliker et al. 2003;Poland et al. 2014;Pietruszka et al. 2015;Chen et al. 2019).
Two types of volcanic tremors have long been recognized beneath the region of Kilauea summit: (1) the first is related to the degassing and spattering of the lave lake surface (Patrick et al. 2016a, b), which can been seen in between 0.3 and 0.8 Hz in the amplitude spectrogram of station UWE (Fig. 2a); and (2) the second is related to the fluid-rock interactions in the conduits, which can be found in the band 0.01-0.1 Hz as very long-period tremor signals (Patrick et al. 2011;Dawson and Chouet 2014).
Geodetically, the spreading of the ERZ is primarily unilateral according to global positioning system (GPS) observations (Owen et al. 2000). In contrast to the relatively inactive northern flank, the southern flank displaces seaward on a north-dipping basal decollement at a depth of approximately 7-9 km owing to the effects of magma injection into the rift zones and the topographic load of the island (Dieterich 1988;Denlinger and Okubo 1995;Dieterich and Cayol 2003;Heliker et al. 2003;Poland et al. 2017).

Ambient noise interferometry (ANI) to measure dv/v
ANI is used to probe the time-lapse changes of subsurface medium in the Kilauea volcanic area. By cross-correlating and stacking sufficiently long time series of seismic data, the noise correlation functions (NCFs) between two stations will in principle approximate to Green's functions as propagating from one station to the other (Lobkis and Weaver 2001;Snieder et al. 2002), including both direct ballistic waves and multiple-scattering coda waves (Sens-Schönfelder and Wegler 2006). Because of multiple-scattering, coda waves have much longer paths than direct waves and are more sensitive to changes in medium and less sensitive to changes in the noise source (Colombi et al. 2014). Based on the assumption of homogeneous velocity changes in the medium, the time shifts (dt/t), if present, between coda waves of two NCFs on different days can be used to characterize the average seismic velocity change (dv/v) of the medium via the equation (Snieder et al. 2002): We analyze 1.5 years of continuous vertical-component seismic records (January 2017 to June 2018) from 12 broadband seismic stations ( Fig. 1) distributed around the summit and along the ERZ. We also include two additional stations (MLOD and HPLD) which are located in northwest and southwest about 10 km away from the summit for a better station coverage. In data preprocessing, we cut the continuous data into 1-day-long segments and remove the instrumental response for each trace. The daily waveforms are demeaned, detrended, tapered, and decimated from 100 to 25 Hz. Amplitude spectrograms are calculated to show the frequency content of the data recorded by stations from the summit toward the MERZ over time ( Fig. 2 and Additional file 1: Figure S1). A clear seasonality that shows lower amplitude in the summer (June-August) is observed in frequencies between ~ 0.1 and 0.3 Hz from the oceanic secondary microseism as previously reported (Donaldson et al. 2017). For the higher frequency range, the volcanic tremor signals from the Halemaʻumaʻu and Pu'u'Ō'ō craters are most (1) dt/t = −dv/v. Feng et al. Earth, Planets and Space (2020) Ballmer et al. (2013). Similar spectrum features can be seen in Station HAT, RIMD, and SDH near the Halemaʻumaʻu crater and Station JCUZ and NPOC near the Pu'u'Ō'ō crater ( Fig. 1 and Additional file 1: Figure S1). We use signals in a frequency band of 0.3-2.0 Hz (black boxes in Fig. 2) for the computation of dv/v. While the tremor signals are inevitably present in our analyzing band, compared to secondary microseism band they are generally stable over time except the February 2017 (Fig. 2c, d and Additional file 1: Figure S1b-e, f ). The effect of tremor sources will be further discussed in "Noise sources" section.
The cross-correlation functions between stations are calculated with Welch's method (Seats et al. 2012) which cuts 1-day waveforms into 5-min moving time windows (75% overlap) to improve the quality of correlation functions. A threshold of 1.5 root-mean-square ratio (RMSR, root mean square of amplitudes in 5-min time window divided by that in the entire 1-day period) is set to remove those RMSR > 1.5 windows that contain unfavorable energetic signals (e.g., earthquakes, instrumental irregularities, and non-stationary transient signals). The rest of the time windows are spectrally whitened, calculated for 100-s time lag cross-correlation functions, and then stacked into a single representative NCF for the day. After daily NCFs are calculated, a further 10-day backward stacking for each day is applied to gain better coherence between NCFs. Stacking more days would increase NCF coherence but lower the time resolution of dv/v measurements. Coherence tests using a 30-s coda window at 6 s after the Rayleigh wave arrival with different stacking lengths are conducted for all three frequency bands and show that almost all coherence curves reach their plateau after 10 days of stacking (Additional file 1: Figures S2, S3). As an example, Fig. 3 shows the waveform coherency of the final NCFs for station pair DEVL-STCD across the entire analyzed period in three frequency bands. Except for February 2017 at 0.9-2.0 Hz, the NCFs generally exhibit stable and high coherence until the eruption in May 2018 (green solid lines). A clear loss of waveform coherency can be seen after the eruption due to structural changes of magmatic system. The waveform coherency of all other station pairs is shown in Additional file 1: Figures S4-S34.
Next, the stretching method (Sens-Schönfelder and Wegler 2006) is adopted to measure dv/v, by which the reference NCF (NCF ref ) is either stretched or compressed to obtain the best fit cc(ε) with the current NCF (NCF cur ) via the equation where ε is a stretching factor, and t 1 and t 2 define the used time windows. By grid-searching ε in a range of − 5% to 5% with a 0.005% interval (2000 trials), the dv/v is determined by the ε which gives the highest correlation coefficient cc between the NCF cur and NCF ref . A reference NCF is constructed by stacking daily NCFs over the entire period of 2017. We compute dv/v in causal and anti-causal part of NCFs individually and average the results. For the station pairs with unstable dv/v fluctuations are excluded in the following analysis (e.g., Additional file 1: Figures S6, S9, S10, S21-S25).

Depth sensitive kernel
Assuming that the early codas of the NCFs are mainly composed of Rayleigh wave energy (Obermann et al. 2013b(Obermann et al. , 2016, we estimate the depth sensitivity of the three analyzed frequency bands using the computer program Surface85 (Herrmann 1987). Given a 1-D velocity model (Fig. 4a) averaged from a recent 3-D velocity model of the Kilauea area (Lin et al. 2014), the sensitivity kernels for shear wave velocity changes in different frequencies between 0.3 and 2.0 Hz are plotted as in Fig. 4b. The sensitivity of the lowest (0.3-0.6 Hz), middle (0.6-0.9 Hz), and the highest (0.9-2.0 Hz) frequency bands peaks at the depths around 1 to 2 km, 500 m to 1 km, and shallower than 500 m from surface, respectively. Weaver et al. (2011) have proposed a theoretical formulation to estimate the root mean square (rms) of the apparent stretching factor ε measured on two waveforms, which can be used as an error indicator of the estimated dv/v: where cc is the correlation coefficient value between the reference and current NCFs of station pairs, T is the inverse of the measured frequency bandwidth, ω c is the central frequency of the measured frequency band, and t 1 and t 2 are the beginning and the end time of the processing coda time window. The errors are computed for causal and anti-causal part of NCFs individually and averaged ( Fig. 5b and Additional file 1: Figures S4-S34).

Noise sources
A stable distribution of background noise sources is one of the most important requirements for using ANI to estimate relative velocity changes (Campillo and Paul 2003;Hadziioannou et al. 2009). In the Kilauea area, the volcanic tremors generated by lava lake spattering are well-recognized signals constantly present in the frequencies of 0.3-0.8 Hz at the Halemaʻumaʻu crater and above 0.8 Hz at the Pu'u'Ō'ō crater (Patrick et al. 2016a, b;Ballmer et al. 2013), as observed in Fig. 2 and Additional file 1: Figure S1. But first, except for few stations around Pu'u'Ō'ō crater showing relatively stronger energy fluctuations in the frequencies above 0.8 Hz ( Fig. 2c and Additional file 1: Figure S1c, e), these tremor signals seem persistent and stable for most of the stations over the study time period. Second, these tremors are generated within the top part of the magma column from lake surface to ~ 1 km depth (Patrick et al. 2016a, b). Since this migration distance is small compared to the length scale of the station pair distance we analyzed, the tremor sources can be reasonably treated as stationary sources. They therefore contribute to the cross-correlation as stable energy sources in the region and unlikely bias the dv/v measurements for the long-distance station pairs used mostly in this study. The effect of tremor migration is probably more significant for short-distance pairs so the crossing pairs around the Halemaʻumaʻu and Pu'u'Ō'ō craters are removed from later dv/v analysis. The seasonality of ocean noise could be another potential source of dv/v errors. However, the seasonal variations in noise energy are seen strongest in 0.1-0.3 Hz ( Figs. 2 and 5d). In the analyzed frequency band of 0.3-2 Hz, the temporally stable noise energy (Fig. 2), high coherency of NCFs over time ( Fig. 3 and Additional file 1: Figures S4-S34), and estimated dv/v not showing related seasonable patterns (Figs. 5 and 6) suggest that the influence from ocean noise seasonality is not significant. Also, Eq. 3 should in principle account for such errors if present, and the estimated errors are generally lower than the dv/v variations discussed later in the study (Fig. 5b and Additional file 1: Figures S4d-S34d).

Estimated seismic velocity variations (dv/v)
The time evolution of estimated dv/v results and corresponding cc values are shown in Fig. 6. For the lower two frequency bands (Fig. 6c-f ), the cc values are stable and above 0.85 most of the time (Fig. 6c, e); for the highest one, the cc values are relatively lower but still above 0.7 except for the February 2017 (Fig. 6a). The variations of dv/v show similar trend for all three frequency bands in 2017 with two velocity drops in March and June (Grayish shaded zones in Fig. 6). These dv/v drops are most obvious in the frequency bands of 0.6-0.9 Hz and 0.9-2.0 Hz. On the other hand, clear dv/v excursions either negative or positive, with slightly decreasing cc values, begin at the end of 2017 for the lowest frequency band (0.3-0.6 Hz). The dv/v in the higher two frequency bands wait until the mid-March 2018 to start to gradually increase. At the end of April, the dv/v in all frequency bands drop rapidly to the eruption on May 2.
To investigate the mechanism of these distinct dv/v features, we use the station pair DEVL-STCD as a representative to demonstrate and compare the estimated dv/v with daily seismic amplitude RMS of each station in different frequency bands, daily peak ground velocity (PGV), and daily earthquake number (Fig. 5). Results for all station pairs are displayed in Additional file 1: Figures S4-S34. In Fig. 5c-e, the 0.01-0.1 Hz and 0.3-2.0 Hz are the frequency bands related to the first and second type of volcanic tremors (Patrick et al. 2011;Dawson and Chouet 2014;Patrick et al. 2016a, b) and the 0.1-0.3 Hz represents the secondary microseism frequency band. We also mark the origin times (grayish dashed lines in Fig. 5) of 16 M 3.5 earthquakes (focal depth less than 10 km) occurring in the Kilauea area prior to the M w 6.9 earthquake on May 4 ( Fig. 1 and Additional file 1: Table S1). It is clear that the two seismic velocity drops in 2017 appear immediately after the earthquakes No. 1-3 (they occurred closely within one day, see Fig. 1 for the locations) and No. 6 (the largest event in the Kilauea area in 2017). These two sets of earthquakes are also the ones that produce the largest PGV (Fig. 5f ) and amplitude RMS in our studied frequency band (Fig. 5c), indicating a correlation to the strong ground shaking of Fig. 4 The 1-D velocity model and corresponding Rayleigh-wave depth sensitivity kernels. a The P-wave (blue) and S-wave (red) 1-D velocity models of Kilauea area averaging from the 3-D velocity model (Lin et al. 2014). b Rayleigh-wave depth sensitivity kernels for shear wave velocity changes in different frequencies between 0.3 and 2.0 Hz earthquakes. Figure 7 further shows the spatial distribution of the two dv/v drops, where the station pairs closer to the epicenters (gray dots) show greater drops of dv/v (warmer colors). For Earthquake No. 6 ( Fig. 7d-f ), the much greater dv/v drops occur in the higher two frequency bands of 0.6-0.9 Hz and 0.9-2.0 Hz, suggesting that the shaking-induced medium changes occur mostly in the uppermost 1 km of the crust based on the depth sensitivity kernels shown in Fig. 4b. This may result from in situ cracks of the shallow crustal layer being temporally increased or opened during dynamic ground shaking and gradually recovered by the crack closure and re-compaction over time Taira et al. 2015Taira et al. , 2018. Similarly, for Earthquake No. 1-3 (Fig. 7a-c), the dv/v reductions are also larger in the higher two frequency bands. For the dv/v excursion that occurs in late December 2017, we confirm that no significant changes in either noise signals (Figs. 2b, c and 5c-f ) or seismicity (Fig. 5g) observed around the time so the excursion is likely related to the actual medium changes. Also, this excursion that appears most obviously in the lowest frequency band and little in the higher two suggests the medium changes at a depth range of 1-4 km according to the sensitivity kernel calculation (Fig. 4b). This depth range well corresponds to depths of reservoirs at summit and the top of the dike system along the ERZ (Heliker et al. 2003;Poland et al. 2014;Pietruszka et al. 2015), we therefore speculate that probable pre-eruptive magmatic intrusion processes responsible for the dv/v excursion beginning from the end of 2017.

The long-term dv/v excursion and magmatic intrusion processes
The long-term dv/v excursions mainly in the frequency band of 0.3-0.6 Hz can form into two groups: one shows negative dv/v excursions as the station pair DEVL-STCD (Fig. 5) and NPOC-PAUD (Additional file 1: Figure S26); and the other shows positive as the station pair HAT-JOKA (Additional file 1: Figure S10) and HLPD-PAUD (Additional file 1: Figure S16) for reference. Based on the temporal characteristics of dv/v, we divide this excursion time into three periods for further investigation (Fig. 6): (1) Period 1 comprises the first stage of long-term dv/v excursion that is only observed in the frequency band of 0.3-0.6 Hz (December 20, 2018 to March 15, 2018); (2) Period 2 is when the summit tiltmeters start to record inflationary ground deformation to the first pressurization issue at Pu'u'Ō'ō called by HVO (March 15 2018 to April 17, 2018); (3) Period 3 is the time after the first

Phase 1: deep magma intrusion
In Period 1 (yellowish shaded zone in Fig. 6), we plot the dv/v of station pairs in lines color-coded by averaging the dv/v values from January 1 to March 15, 2018. Warm and cold colors indicate negative and positive average values of dv/v (Fig. 8a-c). As shown in Fig. 6, the dv/v variations in higher two frequency bands are small (light color in Fig. 8a, b). The main feature of dv/v in Period 1 appears only in the lowest frequency band where the dv/v decreases mainly for the station pairs from the Kilauea summit to the Pu'u'Ō'ō crater (along the UERZ) but increases for the pairs in the surrounding regions (Fig. 8c). In particular, the station pair DEVL-STCD right along the UERZ decreases most down to − 0.3 (Fig. 5).
Its cc values also start to slightly decrease around the same time from December 2017. While a number of causes can influence the cc values, the coincident timing with the dv/v excursion and relatively stable background noise at the times (Figs. 3 and 5c) may suggest structural changes related to the deep magmatic processes (Obermann et al. 2013a;Bennington et al. 2015). The depth of such processes could be reasonably put at the depths of 1-4 km since the sensitivity of the higher two frequency bands is mainly within top 1 km and that of 0.3-0.6 Hz bottoms at 4 km (Fig. 4b).
The dike-like magma reservoir along the ERZ has been suggested to top at ~ 3 km from the decollement of 9 km in previous studies (Delaney et al. 1990;Cayol et al. 2000;Heliker et al. 2003;Poland et al. 2014). A narrow ribbon of precisely relocated seismicity at ~ 3-km depth was interpreted to mark a highly stressed zone corresponding to the intrusion of dikes (Gillard et al. 1996;Rubin et al. 1998). The limited vertical extent of these seismicity suggests that the stress or strain changes are concentrated around the tip of the deforming dike. The observed strong dv/v decreases for the station pair DEVL-STCD likely reflect a stress/volumetric strain change at such localized zone atop the dike tip caused by the magmatic intrusion activities (3-D illustration in Fig. 9a, c). As a connecting magmatic plumbing system, the intrusion likely also feeds in the deeper magma reservoir at the summit (reservoir SC in cross-section AB) and results in the spatial distribution of dv/v decreases mainly between the Kilauea summit and Pu'u'Ō'ō area; but at the same time the inflation of the dike and reservoir also compresses the surrounding medium and causes the dv/v increases in the surrounding regions of the summit and  (Figs. 8c, 9b, c). The insignificant seismic velocity variations in the 0.6-0.9 Hz and 0.9-2.0 Hz frequency bands imply that such magmatic process mainly caused localized stress/strain changes at limited depth range below 1 km. This may also explain why the magmatic process showed little effect on the surface to be geodetically detected before the mid-March 2018 (Neal et al. 2019).

Phase 2: summit inflation
When the summit starts to inflate from the mid-March 2018 (Period 2, bluish shaded zone in Fig. 6), the dv/v in the higher two frequency bands show gradual increases correspondingly (Fig. 6b, c). For this period, we change to plot the dv/v of station pairs by subtracting the dv/v value on March 15, 2018 from that on April 17, 2018 to emphasize the relative velocity changes so zero values mean that the dv/v remain the same with Period 1. In this way, we see the dv/v in the higher two frequency bands clearly increasing for the station pairs at the summit and surroundings (Figs. 8d, e), whereas the dv/v in the lowest frequency band are small in general and suggest a similar sate of the deep magmatic intrusion processes (Figs. 8f and 9e). This pattern implies that the magma has intruded into and inflates the shallower magma reservoir at this time as recorded by tiltmeter and GPS at summit (Fig. 9d). Moreover, these intrusion/inflation processes of the shallow magma reservoir are also evidenced by the decrease of cc values in two higher frequency bands, a sign of structural changes (Fig. 6a, c).

Phase 3: damage accumulation
After April 17, 2018 when the HVO issued notice of the pressurization at the Pu'u'Ō'ō crater (Period 3, reddish shaded zone in Fig. 6), while a few station pairs in the 0.6-0.9 Hz frequency band exhibit a short increase on first few days, most of station pairs start to decrease monotonically to the eruption time in all frequency bands (Fig. 6). In Fig. 8g-i, the dv/v of the station pairs are again computed by subtracting the dv/v value on April 17, 2018 from that on May 1, 2018 in a relative velocity sense. For the lowest frequency band, the dv/v decreases from the summit to the Pu'u'Ō'ō area and strongest along the rift zone, indicating that the deep magmatic intrusion processes become amplified. Two dv/v increasing pairs in the north and south are then in response to the crustal compression from the dike inflation (Figs. 8i and 9g). In contrast, the dv/v for the station pairs around the summit show decreases but not increases as in Period 1 (Fig. 8c, i). The dv/v in the higher two frequency bands also show general decreases around the summit this time (Fig. 8g, h). This decreasing pattern of dv/v, however, contradicts to the amplification of magmatic intrusion aforementioned and also the continued inflationary deformation recorded by tiltmeters at the summit (Olivier et al. 2019). One possible explanation proposed by Olivier et al. (2019) is the accumulating damage induced by the pressure exerted by the magma reservoir on the edifice (Fig. 9f ). In this weakening process, while the surface deformation keeps inflationary, the edifice strength and magma reservoir pressure both decrease with decreasing elastic moduli (Carrier et al. 2015). As a result, it leads to the decrease in seismic velocity and may have facilitated the eastward transport of magma into the rift zone, as amplified dv/v decrease observed in the UERZ in Fig. 8i.
At the end of April, a burst of seismicity occurred and propagated from the Pu'u'Ō'ō crater toward its lower flank. Most likely, these earthquakes initiate the rupturing of a barrier in the MERZ, which allows substantial volume of magma to move into the LERZ and lead to the eruption on May 2, 2018. All derived dv/v then strongly oscillate after the eruption because of the severe loss of NCF coherency due to strong structural change and the boost of earthquake activity (Figs. 5g and 6).

Precursory changes in dv/v of volcano eruptions
The precursory changes in seismic velocity have been found at some volcanoes in recent decades (Brenguier et al. 2008;Obermann et al. 2013a;Budi-Santoso and Lesage 2016;Olivier et al. 2019). At PdF volcano, for example, the precursor time of seismic velocity changes is about few weeks ahead of eruptions (Brenguier et al. 2008). At Kilauea volcano, here we find a longer precursory phase of several months prior to the LERZ eruption.
However, dv/v precursors are not always observed. At Mt. Enta for instance, recent studies showed no significant dv/v precursors to be detected preceding the eruptions (Cannata et al. 2017;De Plaen et al. 2019). Although the used frequency bands and preprocessing methods also matter, such variability of precursory dv/v observations may imply that different types of volcanic eruptions (i.e., explosive, strombolian, and effusive eruptions) have varying precursor times from little to several months in response to their different magmatic processes. Our results suggest a long precursory phase of about 6 months for the Kilauea basaltic shield volcano.

Conclusions
We analyze one-and-half year seismic data with ambient noise interferometry (ANI) to investigate pre-eruptive magmatic processes of the 2018 Kilauea eruption. With careful inspections on correlation quality, noise source, and error estimation, we find two earthquake-induced dv/v drops in 2017 and a magma-related long-term dv/v excursion beginning from the December 2017 (Figs. 5 and 6). The magma-related process comprises three main periods of activity based on the temporal dv/v characteristics (Figs. 8 and 9). From the December 2017 to mid-March 2018 (Period 1), the magmatic intrusion began to take place in the deeper magma reservoir at summit and in the dikes along the UERZ so mainly influenced the dv/v in the lowest frequency band. From March 15 to April 17, 2018, the intrusion process reached the shallow magma reservoir and caused sensible summit inflation recorded by tiltmeters and increases in dv/v around the summit. After April 17 to the eruption (Period 3), the magmatic intrusion was amplified and caused the accumulating damaging of the edifice, which then results in dv/v decreases around the summit and UERZ. As a result, our results suggest a long precursory phase to the Kilauea volcano eruption beginning from the December 2017, with three distinct stages showing progression and interaction of the reservoir-dike magmatic system. The analysis demonstrates the strength of using ANI with multiple frequency bands to gain better depth and spatial information to imaging pre-eruptive magmatic processes over time and can be a useful supplement to current volcanic monitoring systems for eruption warnings and hazard prevention.
Additional file 1. Additional figures and table.