A quest for unrest in multiparameter observations at Whakaari/White Island volcano, New Zealand 2007–2018

The Whakaari/White Island volcano, located ~ 50 km off the east coast of the North Island in New Zealand, has experienced sequences of quiescence, unrest, magmatic and phreatic eruptions over the last decades. For the last 15 years, seismic data have been continuously archived providing potential insight into this frequently active volcano. Here we take advantage of this unusually long time series to retrospectively process the seismic data using ambient noise and tremor-based methodologies. We investigate the time (RSAM) and frequency (Power Spectral Density) evolution of the volcanic tremor, then estimate the changes in the shallow subsurface using the Displacement Seismic Amplitude Ratio (DSAR), relative seismic velocity (dv/v) and decorrelation, and the Luni-Seismic Correlation (LSC). By combining our new set of observations with the long-term evolution of earthquakes, deformation, visual observations and geochemistry, we review the activity of Whakaari/White Island between 2007 and the end of 2018. Our analysis reveals the existence of distinct patterns related to the volcano activity with periods of calm followed by cycles of pressurization and eruptions. We finally put these results in the wider context of forecasting phreatic eruptions using continuous seismic records.


Introduction
Volcano-hydrothermal systems generate a variety of seismic signals, ranging from discrete earthquakes to continuous tremor . Opportunities to study the long-term behavior of these signals, and thus better understand the activity and nature of magmahydrothermal systems, are relatively scarce, because years of continuous data are required. Such opportunities, however, exist for New Zealand volcanoes, allowing to better discriminate background activity from a situation of unrest, explore how a system evolves towards an eruption, and ultimately isolate potential precursors.
GNS Science and its predecessor, the Department of Scientific and Industrial Research, has been collecting continuous seismic data for decades, which have been made available to the scientific community through GeoNet (https:// www. geonet. org. nz/). GeoNet data, together with data from temporary seismic deployments on New Zealand volcanoes by numerous organizations on New Zealand volcanoes, have led to pioneering work in the classification of volcanic earthquakes (Latter 1979;1981); the application of standard methodologies (Titzschkau et al. 2010;Jolly et al. 2012a;Park et al. 2019); and the development of new data processing schemes and monitoring strategies, such as ambient noise seismic interferometry (Mordret et al. 2010;Yates et al. 2019), Open Access *Correspondence: corentin.caudron@univ-smb.fr 1 ISTerre, Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, IFSTTAR , 38000 Grenoble, France Full list of author information is available at the end of the article Caudron et al. Earth, Planets and Space (2021) 73:195 shear wave splitting (Gerst and Savage 2004;Godfrey et al. 2014;Johnson et al. 2011;Johnson and Savage 2012;Keats et al. 2011;Godfrey et al. 2014;Castellazzi et al. 2015), self-organizing maps (Carniel et al., 2013a, b), application of the failure forecast method to tremor data (Chardot et al. 2015), artificial intelligence techniques (Dempsey et al. 2020), the Displacement Seismic Amplitude Ratio (DSAR) technique , and the detection of subsurface stress changes from the timedependent response of volcanoes to external processes (e.g., tides; Girona et al. 2018). One of the volcanoes for which long-term data are available is Whakaari (in Te Reo Māori) /White Island (in English), located ~ 50 km off the northern coast of the North Island in New Zealand, in the eastern Bay of Plenty (Fig. 1). The subaerial portion of Whakaari/ White Island volcano is an andesitic composite system comprising two overlapping cones. The crater of the Whakaari/White Island system hosts persistent super-heated fumaroles (Giggenbach 1987), persistent crater floor outgassing of CO 2 (Bloomberg et al. 2014) and elevated temperature springs (Christenson et al. 2017). The crater is also characterized by variations in the level, temperature, and degree of ebullition of a hot acid lake (Edwards et al. 2017), which testifies to the existence of a long-term (decades to centuries) active hydrothermal system . The hydrothermal system of Whakaari/White Island has been shown to be controlled by the percolation of meteoric and marine water through the volcanic flanks and the interaction of this water with magmatic gas and heat (Christenson et al. 2017). This interaction may be responsible for the more or less persistent seismicity recorded at the volcano, including continuous tremor (Chardot et al. 2015) and discrete events ranging from Very Long Period (VLP; Jolly et al. 2017) to Volcano Tectonic (VT) earthquakes (Nishi et al. 1996).
Whakaari/White Island volcano has been in a near constant state of unrest or eruption in recent history (Kilgour et al. 2019;2021). In particular, in the 14-month period from 5 August 2012 to 11 October 2013, Whakaari/White Island experienced 5 ash and steam eruptions (Chardot et al. 2015), numerous small-scale mud/sulphur and ash eruptions (Christenson et al. 2017;, and the extrusion of a small lava dome in September-November 2013 ). In addition, after more than 2 years of quiescence, a multi-pulse phreatic eruption occurred on April 27, 2016 (Kilgour et al. 2019). The most recent eruption occurred on 9 December 2019 and tragically resulted in fatalities (Park et al. 2020).
In retrospect, precursory seismicity of the 2012-2016 unrest may have occurred on 19-21 August 2011 during a swarm of mixed events comprising Long Period (LP, frequencies 0.5-1.1 Hz), Resonant High Frequency (RHF, frequencies 2-5 Hz) and Very Long Period (VLP, frequencies 0.03-0.125 Hz) earthquakes , which we define later. This activity may have prompted in turn the emergence of continuous seismic vibrations referred to as volcanic tremor in late 2011 (Chardot et al. 2015).
In fact, volcanic tremor at Whakaari/White Island has persisted over more than 25 years of monitoring . Tremor features can include harmonic, non-harmonic, banded or spasmodic characteristics (Chardot et al. 2015;Sherburn et al. 1996), mostly radiating energy between 2 and 5 Hz, but predominantly below 3.5 Hz (Chardot et al. 2015), with variable durations (hours to months). Distinguishing between low amplitude tremor and noise can present challenges often simply linked to the signal-to-noise ratio. Sometimes, the existence of tremor can be demonstrated, for example, by the interaction with tectonic events (Carniel and Tárraga 2006). On the other end, strong periods of tremor at Whakaari/White Island can reach amplitudes of ~ 5 μm/s at station WIZ (Chardot et al. 2015; Fig. 1). Using the source-receiver distance of ~ 600 m estimated by Chardot et al. (2015), this corresponds to a reduced displacement of ~ 35 cm 2 (McNutt 1992). Although many models can explain the generation of tremor at different volcanoes and/or at different times (Konstantinou and Schlindwein 2003), according to Chardot et al. (2015), the most plausible source models at Whakaari/White Island are multiple fracturing involving slow rock failure or fluid flow processes that are consistent with vent lithologies captured in eruption ejecta .
In the following, we present a short summary of the activity of Whakaari/White Island between 2007 and 2018. A more complete review, based on observations of lake level, deformation, gas emissions, fumarolic temperatures, discrete seismic events, and tremor amplitude, can be found in the chronology section of the appendices. Then, we explore potential long-term (weeks to years) medium and source changes using four recently developed data processing schemes: (1) the Displacement Seismic Amplitude Ratio technique ), (2) single station ambient noise seismic interferometry (De Plaen et al. 2016), (3) sensitivity to tidal stresses (Girona et al. 2018), and (4) through an analysis of the spectral properties of earthquakes and tremor ). Our new analyses are then combined with the wider geochemical and deformation time-series and the associated visual observations from the volcano. The origin of these patterns in the magma-hydrothermal system are then interpreted within the context of the volcano hydrothermal conceptual model of Christenson et al. (2017). We finally put these results in the wider context to support improved eruption forecasting for phreatic eruptions using multi-parameter data.

Methods
In this section, we describe the data processing approaches applied in this study using the Whakaari/ White Island seismic data recorded between 2007 and the end of 2018.

Spectral properties of earthquakes and tremor
Different data processing approaches have been applied to the continuous seismic records to build spectrograms and thus retrieve the Power Spectral Density (PSD) and the dominant frequency. The spectrograms are built by dividing the records in 30-min time windows, overlapping by 50%. The PSD of each window is then calculated using MSNoise (Lecocq et al. 2014), which relies on the Probabilistic Power Spectral Densities (PPSD) implementation of ObsPy (Beyreuther et al. 2010;Krischer et al. 2015). The PSDs are calculated using Welch's method, which reduces noise in the power spectra at the expense of reducing the frequency resolution because of frequency binning (Welch 1967). The PSD is binned by 1.25% of an octave (default 12.5%, or 1/8) and is smoothed over 2.5% of an octave (default 100, or 1 octave) around the central frequency of each bin. We have applied less smoothing than in the classic McNamara and Buland (2004) computations (smoothing width of 1 octave with a step of 1/8th octave) to better capture the temporal and spectral changes. The PSDs are transformed from power to amplitude (i.e., square root) and then aggregated to 6-h means (Fig. 2a). We finally normalize each spectrum by the range between minimum and maximum values in each time window (Fig. 2b). This version of the spectrogram loses information about the absolute amplitude time variations but better highlights the often-subtle time evolution of the relative importance of different frequency bands (Carniel 2014).
The main frequency of the continuous recordings is calculated over 5-min, non-overlapping windows. After high-pass filtering the daily seismic data (above 0.5 Hz to prevent any influence by oceanic microseisms), we integrate them to displacement and high-pass filter them once again (above 0.5 Hz). The 5-min windows are then detrended (by subtracting a linear function and removing the mean) and tapered (10% cosine taper). The spectra are then computed using a Discrete Fourier Transform, smoothed in the spectral domain (mean of 40 points) and then normalized by the maximum amplitude. Finally, the main frequency is defined as the mean frequency of the 20 most energetic peaks (i.e., with highest spectral amplitudes). It, therefore, differs from the common definition of dominant frequency as the central frequency of the spectrogram bin, where the maximum is found in a given time window (Carniel 2014).

Displacement Seismic Amplitude Ratio (DSAR) and Real-Time Seismic Amplitude Measurement (RSAM)
Path effects are explored using the Displacement Seismic Amplitude Ratio (DSAR) approach, recently developed by Caudron et al. (2019). Initially designed to study changes in seismic attenuation prior to gas-driven (i.e., phreatic or hydrothermal) eruptions, we test the applicability of DSAR in the context of both magmatic and phreatic eruptions. We begin by pre-processing ) and integrating the raw seismic velocity data. Then, we compute the seismic amplitude in two frequency bands (4.5-8.0 and 8.0-16.0 Hz) and compute the ratio between them (i.e., low over high frequency bands). A change of this ratio may, therefore, indicate a fluctuation in the path region sampled by seismic waves continuously propagating from the source region to the seismic recording site. As shown below, source effects mostly concern frequencies below 4.5 Hz (i.e., we assume negligible source effects above 4.5 Hz).
In addition to the DSAR approach, we computed the RSAM (Real-Time Seismic Amplitude Measurement (Endo and Murray 1991)) typically filtered between 2 and 5 Hz at Whakaari/White Island to focus on volcanic tremor activity (Chardot et al. 2015). In contrast to Chardot et al. (2015)'s approach, we integrate the velocity seismic records after pre-processing them, before bandpass filtering the data. For purposes of consistency with previous studies (Chardot et al. 2015), we use the term RSAM rather than the SSAM naming (Rogers and Stephens 1995).

Single station seismic noise interferometry
We use seismic noise interferometry to estimate relative seismic velocity variations (dv/v) and decorrelation (cross correlation waveform changes) every day using the MSNoise software (Lecocq et al. 2014), focusing on the single station WIZ (Fig. 1). The set of pre-processing parameters described in De Plaen et al. (2016) was applied (clipping to 3 root-mean squared (rms), bandpass pre-filtering between 0.01 and 8.00 Hz) to compute the Cross Correlation Function (CCF), but further tests were made building on the results of Yates et al. (2019). Four different frequency bands were explored, but results above 1 Hz show potential contamination by tremor and were, therefore, discarded to focus on changes between 0.1 and 1 Hz. Both the Moving Window Cross Spectral (MWCS) and the stretching approaches provide similar dv/v estimates (Additional file 1: Figure S1), except that the stretching method gives larger velocity variations than the MWCS approach. The dv/v results were estimated using a reference corresponding to a stack computed using data between 2009 and 2011 when the Each tick on the x-axis corresponds to 1 January of each year volcano was not in unrest. The use of different reference periods has a minor influence on dv/v results.

Tidal modulation of shallow tremor
Variations of the path properties are also explored by analyzing the Luni-Seismic Correlation (LSC), i.e., the evolution of the temporal correlation between seismic amplitude and the fortnightly modulation of the tidal stresses (Girona et al. 2018). In particular, the LSC approach consists of the following four main steps: (i) we calculated the daily median seismic amplitude ( y sam (t) ) and applied a 20-day, 10 corners, Butterworth filter to remove potential long-term variations.
(ii) We modeled the fortnightly modulation of the tidal stresses as a sinusoidal function emulating lunar cycles ( y lun = −cos(2π(t − t low )/T ) , where t low is a reference day of the calendar with neap tide-i.e., quarter moonand T = 14.7653 days is the average time between full moon and the next new moon). (iii) We calculated the Pearson product-moment correlation coefficient ( ρ ) in moving backward windows of 300, 350, and 400 days for comparison. (iv) We calculated the probability of obtaining by chance the observed, or more extreme, values of ρ for each moving window. This approach allows us to explore the response of active volcanoes to tides, which is expected to vary with time depending on the level of ambient stress in the subsurface. The details of the methodology can be found at Girona et al. (2018).

Analysis of discrete earthquakes
Discrete seismic activity at Whakaari/White Island can be categorized into distinct event types, including LP, RHF and VLP events (see Additional file 2: Figure S2 for examples of waveforms and corresponding spectra). Although some of the literature considers 'high frequency events' to be a synonym of volcano-tectonic (VT) earthquakes (e.g., McNutt 2005), we follow the nomenclature used by Jolly et al. (2017) and associate the term RHF activity to its spectral content only, thus removing reference to particular source processes.
We are specifically interested in near-vent discrete seismic events. Hence for VT earthquakes, we select events with broadband spectra (Additional file 2: Figure S2) and distinct P and S phases. These events are probably related to double-couple fault stick slip source processes and a wide range of S-P times are observed at Whakaari/ White Island. For our purposes, we require S-P times of less than 1 s (at station WIZ) for the catalogue to assure that the events are near the volcanic edifice and, therefore, probably related to near vent volcanic activity. LP and RHF activity are also distinguished: LP events have spectral peaks from 0.5 to 2 Hz, while RHF-type earthquakes range from 2 to 5 Hz (Additional file 2: Figure S2).
Note that RHF seismicity can be distinguished from VT seismicity by the strongly peaked spectra of the former, which are distinct from the broad spectra of VT event types. VLP events have the lowest frequencies for local Whakaari/White Island seismicity (~ 0.1 to 0.03 Hz, Park et al. 2020).
We used two approaches to detect the activity at Whakaari/White Island. For VT, LP and RHF seismicity, we used an STA/LTA (Short Term Average/Long Term Average) algorithm, where STA = 1 s, LTA = 5 s, and the threshold STA/LTA > 3.0 for passband filtered data. The specified passband is 0.03 to 0.1 Hz for VLP, 0.2 to 2 Hz for LP, and 2-5 Hz for RHF, consistent with the spectral content described above. This algorithm produced a significant number of false detections, hence visual inspection of each waveform recorded between 2007 and 2018 was required.
For VLP event types, we used a waveform semblance approach using the two permanent GeoNet seismic stations (WIZ and WSRZ). This method searches VLP signals that occur beneath the vent area using the waveform coherence between seismograms. Previously, the approach was applied to the 2016 eruption episode by Jolly et al. (2018), while we extended it here to the period April 2013-December 2018. The detected events show two repeating waveform patterns; hence, we performed a clustering analysis, whereby we identify two VLP earthquake families. We then used a template matching method to detect the families that occurred prior to April 2013 when only one station (WIZ) was operated at Whakaari/White Island. The families are characterized by 'mirror image' reverse waveform polarities and, therefore, possibly genetically related. The details of this analysis are presented in Park et al. (2020) (this issue).

Results
Two eruptive periods occurred in Whakaari/White Island between 2007 and 2018. The 2012-2013 explosive eruptions included lava dome extrusion , steam and ash plumes, as well as ballistic ejecta and pyroclastic flows (Kilgour et al. 2019). In hindsight, the first precursor might have been a swarm of earthquake that occurred on 19-21 August 2011, followed by an overall increase in RSAM and SO 2 and H 2 S gas flux . In contrast, the 27 April 2016 eruption was a multi-phase isolated phreatic eruption that generated ballistics and a pyroclastic surge within less than an hour. Few short-term precursors could be detected prior to this event even in hindsight Dempsey et al. 2020).
A detailed description of the key periods of quiescence, unrest and eruptive activity can be found in the chronology section of the appendices (Additional files 3, 4: Figures S3 and S4). Below, we focus on potential longterm (weeks to months) precursors rather than shortterm (hours to a few days) precursors, which have been recently detected for some eruptions at Whakaari/White Island using machine learning approaches (Dempsey et al. 2020).

Spectral properties of earthquakes and tremor
To get an overview of the typical spectral features recorded at Whakaari/White Island, we computed the PSD using data from station WIZ (Fig. 1). The first clear change in seismic noise properties appeared after the earthquake swarm in August 2011 , which led to a new regime that persisted until the 2016 eruption ( Fig. 4), with most of the energy radiated below 5 Hz. The signal of this new regime is hereafter referred to as tremor, and probably reflects the activation of one or more sustained and coherent seismic sources that remained active between 2011 and 2016.
We also computed the dominant frequency of tremor between 1 and 5 Hz (Fig. 3). From 2007 to 2011, the main frequency of tremor coincides with our lower frequency bound (1 Hz, Fig. 3), hence most likely reflecting oceanic microseism activity. Scattered peaks in the main frequency between 1.5 and 3.5 Hz appeared after the 19-21 August 2011 earthquake swarm, until the 2016 eruption ( Fig. 3). Gliding spectral lines were observed from 2012 until the beginning of 2013, coincident with the dome extrusion period . After the beginning of 2013, the dominant frequency became more scattered until the April 2016 phreatic eruption (red line in Fig. 3). Following the eruption, the situation recovered to the pre-2011 situation with spectra probably dominated by oceanic activity.

Seismic noise-medium/path effects DSAR
A pronounced increase in DSAR from around the range of 2 to 5 occurred between May 2012 and the 2012 eruption (orange curve in Fig. 4a). The DSAR values remained elevated until the beginning of 2013 when they suddenly decrease to 2.5-3-still remaining above the pre-eruptive The 1-day data are shown as black dots. We also smoothed the 1-day data with a Hendrick Prescott filter, which prevents temporal aliasing (blue line (parameter = 1000) and orange (parameter = 10 8 ) lines). Note that the y-scales are different. See Fig. 4 for the description of the vertical lines corresponding to key time periods. Each tick on the x-axis corresponds to 1 January of each year values ( Fig. 4). Unlike standard seismic summary parameters (RSAM and earthquake activity (Additional file 4: Figure S4), the DSAR values remained elevated (~ 3-4) between 2013 and a few days before the 2016 phreatic eruption, as do the SO 2 emissions (Additional file 3: Figure S3c). DSAR values then dropped below 2 right after the eruption. Additional file 5: Figure S5 displays the time evolution of DSAR and the individual frequency bands (LF (4.5-8 Hz) and HF (8-16 Hz)).
No obvious precursor has been observed prior to the 2012 eruption using the DSAR approach, in contrast to case studies at Ruapehu, Tongariro and Kawah Ijen ), but the main change of trend occurred at the time of the 2011 swarm of earthquakes (orange curve, Fig. 4a). A sudden increase of DSAR also occurred a few months before the 2012 eruption, but similar short-term changes had already been observed previously [e.g., in 2007 ( Fig. 4a)]. As highlighted by Caudron et al. (2019), short-term changes (a few days) in DSAR values should, therefore, be considered with caution.

Luni-Seismic Correlations (LSC)
A previous study on Ruapehu volcano (Girona et al. 2018) suggests that a changing luni-seismic correlation (i.e., correlation between seismic amplitude and the fortnightly modulation of the tidal stresses, related to lunar cycles) can inform of changing stress conditions in the subsurface. Our analysis on Whakaari/White Island data reveals the consistent emergence of a highly significant (above 5σ confidence level) sensitivity to tidal stresses in the 2015-2016 period (Fig. 5), reaching negative correlation coefficients of up to ~ 0.50, ~ 0.43, and ~ 0.38 for 300-day, 1-year, and 400-day moving windows, respectively (Additional file 6: Figure S6). The correlation coefficient, in contrast, remains below ~ 0.20 for the rest of the period analyzed (Fig. 5).

Relative seismic velocity variations
The long-term increasing dv/v trend between 0.1 and 1 Hz described and interpreted by Yates et al. (2019) is also observed in this study between 2008 and 2016 ( Fig. 6a). Following the 2016 phreatic eruption, seismic velocity increased but with greater values in the 2017-2019 time period and the CCFs became highly decorrelated ( Fig. 6b). We do not interpret the results during the 2012-13 eruptive sequence as they could be contaminated by eruptive tremor sources (Gómez-García et al. 2018). The results for the different combinations of components (EN, EZ and ZN, where Z is the vertical component) are provided in Additional file 7: Figure S7.

Discrete earthquakes
We summarize below the main findings as some aspects have been presented in the chronology above and in the appendices. From the discrete earthquake analyses (Additional file 4: Figure S4), the primary features include: (1) VTs (S-P < 1 s) are not related to eruptive activity (Additional file 4: Figure S4c); (2) generation of LPs usually stops after an eruption sequence and may relate to low levels of tremor as mentioned in the chronology   Figure S4d); (3) RHF activity appears to be related to rapid lake loss and may reflect heating episodes (compare Additional files 3, 4: Figures S3a and S4a); and (4) VLP were related to the onset of unrest in mid-2011 and are associated with eruptions (Additional file 4: Figure S4b).

Discussion
Significant changes in seismicity were recorded at Whakaari/White Island between 2007 and 2018, providing further insights into the volcanic processes at play. We will first discuss the advantages and limitations of each parameter as well as their relation to possible triggering mechanisms. We will then revisit the chronology in light of these events, and discuss the values of these parameters for long-term forecasting. Finally, we will discuss our results in the broader context of long-term forecasting of phreatic eruptions at volcanoes.

Tremor source changes
We only detect tremor following the 2011 earthquake swarm and it remains clearly active until the 2012-13 eruptions. Tremor does not disappear completely after the 2012-13 sequence of eruptions, but its main frequency shows some variability until the 2016 eruption when it abruptly stops (Fig. 2). A violin plot illustrates the variability associated with both unrest and eruption periods (Additional file 8: Figure S8). Excepting the dome building episode , the tremor is generally broadband with the energy mostly distributed between 2 and 5 Hz. Chardot et al. (2015) noted that the main frequency shifts from lower (2.0-2.5 Hz) towards higher frequencies (3.0-3.5 Hz) during the 2012 dome growing episode ). Our study shows that this situation persisted when the dome growth stopped until the 2016 eruption and is, therefore, not a diagnostic feature of dome building episodes. It could instead reflect a permanent structural change and gravitational loading. Chardot et al. (2015) also argued that the low-frequency tremor signal (2-5 Hz) could be used as a proxy for fracturing. Due to the relatively large source-receiver distance (~ 600 m), higher frequencies (> 5 Hz) are expected to be more attenuated before reaching station WIZ. Our study is in line with their discussion of possible models (fracturing vs resonance/magma wagging/ boiling), and cannot better constrain either of these possibilities. Kennedy et al. (2020) also recently highlighted the role of fractures, and efficient fracture opening and closing in response to changes in effective pressure as a control on fluid flow and in the generation of geophysical and geochemical signals at Whakaari/White Island volcano.
Unlike the method developed by Chardot et al. (2015), our data processing approaches would not be useful for short-term forecasting purposes, but can highlight periods of elevated activity or dome building episodes.

Medium changes DSAR
The overarching hypothesis to interpret DSAR variations as resulting from path or medium effects relates to the absence of source effects above 4.5 Hz at the scale of days, weeks, months or years ; an assumption seemingly valid at Whakaari/White Island excepting during the 2012-13 eruptive sequence and a few weeks before the 2016 phreatic eruption ( Fig. 2 and Chardot et al. 2015). Results at the WSRZ station (Fig. 4b) show similar trends to WIZ station ( Fig. 4a) with elevated values until the 2016 eruption, followed by a significant drop and re-increase in 2017. These similar trends suggest that site effects seem to be of secondary importance. The DSAR parameter may, therefore, reflect medium changes in the shallow part of the hydrothermal system.
Elevated DSAR values are thought to reflect high attenuation in the shallow subsurface, perhaps through crustal gas accumulation due to enhanced supply of volatiles from depth, or due to enhanced sealing or crack closing leading to gas accumulation, as also suggested by relative velocity increase derived from ambient noise interferometry (Yates et al. 2019). Recent laboratory experiments have measured high attenuation in gas-saturated samples consistent with hydrothermal systems, where fluid saturation and overpressure occur (Clarke et al. 2020). Importantly, high attenuation associated with gas-saturated tuffs produces frequency spectra containing lower frequencies, consistent with elevated DSAR ratios.
Low gas emissions and seismicity were recorded between November 2013 and February 2015, which is characterized as a quiescent period (see chronology and Additional file 3: Figure S3c, d). In contrast, the period of strong gas emission in October 2015 could indicate less accumulation of gas in the subsurface, consistent with the DSAR decrease from 5 to 3 until before the 2016 phreatic eruption. Hence, as shown elsewhere , this parameter appears predominantly sensitive to enhanced attenuation due to the accumulation of volatiles in the subsurface.
Among the complementary parameters available, lake level percentages derived from InSAR after 2015 (Additional file 3: Figure S3a) show a remarkable agreement with DSAR values (Fig. 7). Periods of crater floor inflations correlate with lake level decrease, likely resulting from hydrostatic loading distributed through the crater floor (Christenson et al. 2017). Such correlation also suggests that DSAR variations at Whakaari/White Island are primarily driven by changes in the shallow hydrothermal system rather than deep-seated processes. This is an interesting finding as other seismic parameters, such as RSAM, do not correlate with lake behavior (Additional file 3: Figure S3e and Christenson et al. 2017).
In the case of Kawah Ijen, Ruapehu and Te Maari volcanoes , the DSAR variations were interpreted as due to enhanced attenuation in the shallowest parts of the volcano-hydrothermal system (< 1 km below the craters). DSAR absolute values are closer to the ones observed at Kawah Ijen volcano prior to phreatic eruptions (4-6), than at Tongariro and Ruapehu volcanoes (1-2) . Jolly et al. (2012b) found very strong attenuation (Q < 10) at Whakaari/ White Island for both consolidated and unconsolidated rocks. Q values at Kawah Ijen are not known but high attenuation (low Q values) is likely.

Luni-Seismic Correlation
The possibility that tides can influence volcanic activity has been questioned (e.g., Neuberg (2000)) although recent studies support that tides can at least modulate the dynamics of volcanoes that are erupting or in a critical pre-eruptive state (Girona et al. 2018;Dumont et al. 2020). The LSC approach reveals that shallow volcanic tremor at Whakaari/White Island was sensitive (beyond the 5σ confidence level) to tidal stresses during 2015, i.e., for several months before the 2016 phreatic eruption (Fig. 5). Similar results were found prior to the 2007 phreatic eruption of Ruapehu volcano (Girona et al. 2018), thus suggesting that LSC increases when the stress conditions of the near subsurface remain at critical values for months or years before gas explosions. Girona et al. (2018) interpreted the pre-eruptive increase of LSC at Ruapehu as a result of the partial or total sealing of gas pathways, which decreases the effective permeability of the shallow crust and hinders the escape of gases from the volcano-hydrothermal system. Sealing was shown through a new mechanistic model (Girona et al. 2018 to be able to change the resonant properties of gasleaking subsurface cavities, causing resonance and thus tremor to become sensitive to external perturbations, such as tidal stresses, when the effective permeability of the shallow crust decreases below a given threshold. We, therefore, propose that a similar sealing mechanism is responsible for the increase of LSC during 2015 at Whakaari/White Island. This is consistent with the increase of DSAR during the same time period, although we cannot reject the possibility that other mechanisms play an important role in the observed signal (e.g., tidal stresses may enhance bubble coalescence in bubble columns; Dinger et al. 2019).

dv/v and decorrelation
The relative seismic velocity overall increased by 0.6% between 2008 and 2019, and at a much higher rate from 2017 until 2018 (Fig. 6). Seismic velocity may increase through the closing of cracks through pressure increase at depth (Donaldson et al. 2017;Yates et al. 2019). A pressure increase can, however, produce both a dv/v increase and decrease in different areas depending on the source geometry and depth, as well as the locations of the sensors (Budi-Santoso and Lesage 2016; Donaldson et al. 2019Donaldson et al. , 2017. As shown by Yates et al. (2019) using the same frequency band, dv/v increased prior to the 2016 phreatic Fig. 7 Comparison DSAR (green) and percentage of lake level (blue) relative to full, i.e., elevated values correspond to higher lake levels (Hamling 2017). See Fig. 4 for the description of the vertical lines corresponding to key time periods. The 1-day DSAR data have been smoothed with a Hendrick Prescott filter, which prevents temporal aliasing (parameter = 10 4 ) eruption during a period of uplift at shallow levels, implying increased pore pressures in the shallow hydrothermal system. Our results tend to confirm this finding and show even more elevated dv/v between 2017 and 2019, coinciding with uplift and refilling of the lake (Hamling 2020).
Increases in seismic velocity prior to eruptions have been less frequently reported than drops (e.g., Brenguier et al. (2008)), consistent with crack closure and sealing as hypothesized by Kennedy et al. (2020). Using hydrothermal breccia ballistics, Kennedy et al. (2020) recently showed that permeability within the conduit breccia is more confining-pressure-dependent than within the surrounding edifice. The pre-existing cracks in hydrothermally mineralized zones, therefore, more easily and efficiently close and open when the advection of fluids overcome confining pressure and increase pore pressure. The efficiency of this process shown in hydrothermally altered conduits would be manifested in dv/v changes. The processes altering seismic wave propagation may, therefore, be extended to similarly 'wet volcanoes' , where weak alteration minerals allow efficient crack opening and closing. Finally, the dv/v for the EN cross-components shows larger variations, which may be caused by different waves contributing to the wavefield (dominated by Love waves) or anisotropy (Machacca-Puma et al. 2019). This will be the topic of future studies.
A puzzling decorrelation also commenced shortly after the April 2016 eruption (Fig. 6). Waveform decorrelation of CCF has generally been related to structural changes, such as displacements of scatterers or opening/ closing of cracks (Obermann et al. 2013;Machacca-Puma et al. 2019). Thery et al. (2020) recently observed greater decorrelation, while water level was increasing. Their results indicated that structural changes at the grain scale were more important than mechanical changes (e.g., fracturing (Obermann et al. 2013)). They concluded that decorrelation was sensitive to pore filling, in contrast to seismic velocity, which is weakly impacted by the presence of fluids.
In conclusion, dv/v increases prior to phreatic eruptions may relate to weak alteration minerals that allow efficient crack opening and closures. Whereas small dv/v increase preceded the 2016 eruption, the dv/v increase between 2016 and 2018 was much more significant and accompanied by a decorrelation consistent with pore filling by alteration minerals.

2007-2018 chronology
Several notable changes in seismic noise-related parameters occurred between 2007 and 2018. Following a period of quiescence, the August 2011 earthquake swarm marked the onset of persistent tremor radiation (Fig. 2) and the end of an 8-year decrease of DSAR (Additional file 9: Figure S9). The DSAR values dramatically increased shortly before and during the 2012-13 eruption and energetic tremor was detected throughout the time period from August 2011 until the 2016 phreatic eruption. The DSAR values were particularly elevated during the 2013-2016 time period, while the LSC only correlated with RSAM in 2015-2016 (Fig. 9). We finally observe a decorrelation, DSAR and dv/v increase in mid-2017 (Fig. 9). Below, we investigate the origin of these changes.

until 2011
Whakaari/White Island has been monitored using a wide range of sensors and parameters. Deformation of the main crater floor has been measured for decades at Whakaari/White Island and correlates remarkably with lake level (Christenson et al. 2017). Modelling efforts have shown that deformation results from shallow hydrothermal pressurization due to hydrostatic loading on the fumarolic conduits at 10 s-100 s meters depth (Peltier et al. 2009;Fournier and Chardot 2012;Christenson et al. 2017).
We processed all the existing seismic data using the DSAR approach and found values around 1 in 1989, but only based on a few days of recordings, and a longterm decreasing trend between 2001 and 2011 (from 5-6 to 1-2, Additional file 9: Figure S9). This declining trend mimics the inflation from 2003 onward (Christenson et al. 2017) driven by a 100-m depth shallow source (Fournier and Chardot 2012) and ascribed to hydrostatic loading of the lake basin and fumarolic system (Christenson et al. 2017), rather than the build-up of volatiles (e.g., Werner et al. 2008). A minor thermo-elastic component has been attributed to the deformation as well (Christenson et al. 2017). The DSAR decreasing trend would preclude gas accumulation as the dominant source of the deformation signal as it would lead to higher attenuation in the shallow system ) and thereby elevated DSAR values.

until the 2016 eruption
The first important tremor change occurred in 2011 following a multi-phase earthquake swarm (Additional file 4: Figure S4). The VLP swarm in 2011 may have marked the onset of mobilization of fluids within the upper portions of the magma filled conduit system at ~ 1 km. The perturbation of this system also modified the overlying hydrothermal system. Tremor may have been triggered by enhanced fracturing through the newly mineralized carapace, possibly reflecting magma intrusion at depth. This period also marks the end of longterm DSAR decrease ( Fig. 4 and Additional file 8: S8).
No obvious long-term (weeks/months) precursor is found before the onset of the 2012-13 eruptive sequence. Most of our approaches are not suited to dissect eruptive sequences; eruptive sources would pollute noise-related approaches designed to probe the medium. Instead, the methods outlined here may sample characteristics of unrest and we note that seismic velocity stopped increasing at the end of the eruptive sequence (Fig. 6a). We do not interpret this further as CCFs could be polluted by eruptive tremor thereby impacting the reliability of our dv/v estimates.
The quiescence period between 2013 and mid-2015 is associated with elevated DSAR values (> 4, Fig. 9) and stable dv/v (Fig. 9) which may reflect gas accumulation in the shallow subsurface or more input from depth. A notable change in gas compositions and temperature of the F0 fumarole occurred in mid-2015, wherein temperatures increased from 154 °C to 180 °C (Additional file 3: Figure S3b). The CO 2 /SO 2 ratio variations suggest low-level scrubbing of SO 2 . Strong gas emissions and seismicity (RHF and LP activity (Additional file 4: Figure  S4)) and a tiny dv/v increase (Fig. 6) were also registered in mid-2015. Taken altogether, these observations suggest pressurization in the main vent. Interestingly, this period of pressurization corresponds to the only statistically significant (5σ confidence level) correlation between tremor and lunar cycles (Fig. 5), as observed prior to the 2007 phreatic eruption at Ruapehu (Girona et al. 2018). Together, this indicates that the volcano subsurface was critically stressed in the months preceding the 2016 eruption (Christenson et al. 2017). While LPs have been noted earlier to relate in general to a low tremor period before the 2012-2013 eruptions, one exception exists. Prior to the 2016 eruption higher tremor is accompanied by LP seismicity, and this activity suddenly stops coincident with the 2016 eruption. In this case, it is surmised that depressurization of the hydrothermal system must also relate to LP production. The hypothesized pre-eruption conditions may support Dempsey et al. (2020)'s suggestion that the April 2016 eruption had a different triggering mechanism from prior eruptions.

Post-2016 phreatic eruption period
Tremor vanished below the detection level after the 2016 phreatic eruption (Figs. 2, 3). In contrast to tremor, seismic interferometry revealed greater decorrelation at the end of 2016-beginning of 2017, while the water level was increasing (Fig. 8) and seismic waves were propagating faster (Fig. 9). The decorrelation could indicate pore filling by minerals rather than the typical structural changes, such as displacements of scatterers or opening/closing of cracks (Machacca-Puma et al. 2019 ;Obermann et al. 2013). Around the same time, the DSAR increased again and uplift resumed, indicating a period of enhanced pressurization (Fig. 9).
While most observations indicate changes in the very shallow portions of the subsurface (≤ 1 km of depth), the decorrelation may reflect the injection of fluids at greater depths. For example, stress changes due to the movement of magma have been hypothesized as early precursors 10-5 months prior to eruptions at Ruapehu volcano (Hurst et al. 2018;Kilgour et al. 2014). Fluid movement was invoked to explain changes in anisotropy and b values before the 2006 and 2007 eruptions at Ruapehu  (Keats et al. 2011), whereas progressive seal formation between 2006 and 2007 might have pressurized the gas column beneath the seal. This period of decorrelation at Whakaari/White Island interestingly coincides with a substantial and sustained CO 2 degassing event that commenced in early 2017.

eruption
The most recent eruption occurred in December 2019. A swarm of VLP events occurred in 2018 and evolved into elevated RSAM beginning in 2019 (Park et al. 2020). The RSAM unrest was then followed by the eruption on 9 December 2019. Dempsey et al. (2020) designed an automated precursor recognition method that retrospectively detected the unrest at least 4 h before the eruption based on the RSAM.
Our new set of approaches shows interesting long-term changes prior to the December eruption (Fig. 9). The DSAR kept increasing until the eruption. Following the most significant long-term dv/v increase, a sharp drop occurred until the eruption. The most intriguing change is the progressive decorrelation that initiated shortly after the 2016 eruption. These observations might indicate the ingress of fluids in the upper portion of the volcano after the 2016 phreatic eruption. The sharp DSAR increase, dv/v decrease, further decorrelation, and absence of tidal modulation during the few months preceding the eruption could reflect increased gas inputs rather than permeability-controlled pre-eruptive activity. This is consistent with the observed uplift (Hamling 2020), increased gas emissions, earthquakes and RSAM prior to the eruption in December 2019 (Park et al. 2020). If complemented by other observations, these new techniques could help interpreting subsurface processes and forecasting phreatic eruptions. Table 1 summarizes the main results for the periods covered in this study. We specifically focus on long-term (weeks to months) behaviors rather than short-term behaviors [hours to days (Dempsey et al. 2020)].

Forecasting unrest and eruptions at Whakaari/White Island
The 2012-13 eruptive sequence had been preceded by a swarm of earthquakes that occurred on 19-21 August 2011, followed by overall increase in SO 2 and H 2 S gas flux and RSAM ( Additional files 3, 4: Figures S3, S4 and Jolly et al. 2017). The analysis of PSD and dominant frequency of tremor complements the RSAM (Fig. 2). In addition, the PSD allows us to detect gliding during the dome extrusion period, whereas the violin plot (Additional file 8: Figure S8) quickly highlights unrest. The other less conventional parameters (DSAR, LSC, dv/v, decorrelation) do not significantly change before the 2012-13 eruptive sequence with the exception of an increase in DSAR (to a value of 5) 2 months before the eruption (Additional file 9: Figure S9) following an 8-year decrease. As a result, prior to eruptions including magmatic emissions, these parameters are overall low.
In contrast, most of them are unusually high (DSAR ~ 3-4, significant LSC in 2015, dv/v increase) prior to the 2016 phreatic eruption, and also in 2017-2018. This new set of parameters indicate that the medium was highly pressurized for months to years prior to the phreatic eruptions. Elevated periods of DSAR and dv/v in particular have been always followed by eruptions, then abruptly decreased. This is valuable as phreatic eruptions are notoriously hard to forecast. Dempsey et al. (2020) tested their forecasting approaches using the DSAR parameter, but its implementation for long-term forecasting is the topic of future work. Combining both long-term and short-term proxies holds great potential for forecasting phreatic eruptions. The decorrelation which started in 2017 and only preceded the 2019

Comparisons with other volcanoes
Following the pioneering reviews of Barberi et al. (1992) and Browne and Lawless (2001), Stix and de Moor (2018) recently compiled and examined existing data associated with various phreatic eruptions and highlighted two endmembers: eruptions fed by deep hydrothermal systems, where magmatic gases are sealed creating overpressure (type 1) and open-vent shallow degassing systems, where vaporized liquid water drive eruptions (type 2), with both systems potentially undergoing sealing at various timescales and both type of eruptions being possible at the same system (e.g., Poás volcano). Phreatic eruptions can be triggered by magma intrusions (Nakamichi et al. 2009), injection of gas (Christenson et al. 2017) or infiltration of meteoric water (Gaete et al. 2019). Although not restricted to these processes, hydrothermal sealing and failure has recently been recognized as pivotal in driving phreatic eruptions as it can lead to decreased vent porosity and permeability, and pressurization followed by sudden decompression (Christenson et al. 2010;Mayer et al. 2015;Rodgers et al. 2015). Yet these processes remain hard to detect. Promising results to detect formation of seals, presumably made of sulphur and sulphate phases, have been shown using geochemical ratios prior to phreatic eruptions (Battaglia et al. 2019;de Moor et al. 2016de Moor et al. , 2019. Increase in CO 2 / SO 2 and reduced surface S degassing have been used to track sealing at Rincón de la Vieja (Battaglia et al. 2019). However, results at Poás volcano show precursory trends toward more S-rich gas compositions (de Moor et al. 2016). These differences have been interpreted as indicative of lower magmatic gas input or more effective sealing at Rincón de la Vieja (Battaglia et al. 2019). Alternatively, an increase in CO 2 /SO 2 and H 2 S/SO 2 for type 2 systems could indicate an expansion of the liquid-dominated hydrothermal systems.
Geophysical exploration of sealing and opening processes has been rather limited. At Telica, high frequency event swarms appear as precursors to some phreatic eruptions (Rodgers et al. 2015) and may be indicative of fracture associated with increased pressurization (Roman et al. 2019). Both changes in shear-wave splitting and the ratio of compressional to shear wave speeds (Vp/Vs) were associated with degassing prior to a summit eruption at Kilauea Volcano, Hawaii (Johnson and Poland 2013). If sealing occurs, changes in path effects may be inferred from DSAR patterns, although DSAR absolute values seem influenced by local attenuation conditions (high DSAR with low Q values). Short-term fluctuations at Tongariro and Ruapehu were much more scattered, and so are difficult to compare with Kawah Ijen and Whakaari/White Island. Yet Kawah Ijen and Whakaari/ White Island DSAR patterns suggest sealing over timescales of weeks to months, which may relate to rapid mineral precipitation (sulphur, sulphates, silica) or sufficient pore space to accommodate precipitation (Christenson et al. 2017). A future comprehensive study should be performed, since sealing is largely influenced by pH, temperature, redox conditions, pressure, pore structure, and crystal and glass content of original lithology (Christenson et al. 2010;de Moor et al. 2019;Gaete et al. 2019;Mordensky et al. 2018;Roman et al. 2019;Stix and Moor 2018). Little is known about the timescales of these processes that may vary from decades down to days (Stix and de Moor 2018). Due to low fragmentation threshold in hydrothermally altered material, Kennedy et al. (2020) recently showed that mineralized zones (alunite, anhydrite) at Whakaari/White Island may open/close more rapidly, in seconds to hours, when pressurized fluids travel towards the surface. The physical properties of a hydrothermal system will vary in both space and time, and at any one time the deeper part of the hydrothermal system may be behaving differently to the shallow part. For example, minerals may be precipitating at certain locations densifying the rock (increasing seismic velocity, reducing attenuation and strengthening the rock) (e.g., Heap et al. (2019)) while dissolving at other locations reducing the density of the rock and providing space for fluids (decreasing the seismic velocity, increasing attenuation) (Clarke et al. 2020;Mordensky et al. 2018). In addition to this, both the minerals, e.g., clay vs silica ) and the fluid (Clarke et al. 2019) will affect the elastic and attenuation properties, thereby simultaneously driving opposing changes in seismic attenuation at different levels, and affecting source trigger mechanisms. Systematic comparison with seismic velocity variations, as done in this study, and with tremor-related parameters or the levels of RHF vs LF seismicity (Rodgers et al. 2015) is desirable, especially when complemented by gas and deformation data, which would provide critical insights into the changes occurring at different levels in these systems. In addition, a new methodology to analyze the release of heat emissions using satellite-based infrared data has detected precursory signals before the phreatic eruptions of Ontake and Ruapehu volcanoes (Girona et al. 2021). The changes have been attributed to enhanced hydrothermal activity. This approach could, therefore, provide an additional observation to probe volcano-hydrothermal dynamics. Stix and de Moor (2018) also noted that VLPs and tremor appear to be reliable indicators of pressurization, but with variable timescales. However, some exceptions exist, e.g., a phreatic eruption at Merapi in 2014 that occurred without any VLP (Métaxian et al. (2020)). Our study and Park et al. (2020) highlighted that VLPs at Whakaari/White Island were also closely related to periods of elevated activity. However, our study stresses the need to further explore and process continuous seismicity with novel methodologies. Seismic interferometry can be carried out using the MSNoise software (Lecocq et al. 2014), whereas PSD/spectrograms, DSAR and sensitivity to lunar cycles could become standard products, some of them being already plugged as modules to the MSNoise software. The time evolution of all these parameters could be jointly analyzed also using a supervised or unsupervised machine learning approach (Carniel and Guzmán 2020; Dempsey et al. 2020). We note, however, the need to have stations located nearby the locus of phreatic activity as phreatic-related processes are typically subdued and are easily attenuated by poorly consolidated media, resulting in low signal-to-noise ratios at greater distances.

Conclusions
Our study brings a new perspective into the Whakaari/ White Island volcano-hydrothermal system. By taking advantage of an exceptionally long continuous time series, we processed seismic data from a single station using seismic noise-based approaches. The combined analysis of many different parameters has great potential to provide insights into volcanic activity (Fig. 9).
Tremor generally radiated seismic energy between 2 and 5 Hz, and was closely associated with the period of unrest ranging between 2011 and 2016, probably driven by magma intrusion. Our study concurs with the Chardot et al. (2015) and Kennedy et al. (2020) discussions of models, with the possibly predominant role of fracturing and mineral precipitation in the generation and path modification of geophysical and geochemical signals at Whakaari/White Island volcanoes. Although our data processing approaches would not be useful for short-term forecasting purposes, they can highlight periods of elevated activity or dome building episodes and can serve as a technique to assess volcano alert levels, complementing the short-term forecast provided by the approaches of Chardot et al. (2015) or Dempsey et al. (2020).
Increases in seismic velocity are ascribed to gradual mineralization and pressurization in the subsurface over months to years, between 2008 and 2016 and to a much larger scale between 2016 and 2018. Waveform decorrelation only initiated after the 2016 eruption, concurrent with increased CO 2 emissions, and could reflect unusual injection of fluids at greater depth.
Further work is, however, needed to fully explore this intriguing observation and reliably interpret it.
Provided that seismic sources are generally not radiating above 5 Hz, the DSAR parameter could serve as a key proxy to detect periods of enhanced attenuation in the subsurface. Long-term trends are similar to deformation patterns derived from ground and satellite-based observations. These trends also show that system sealing was probably not efficient enough to develop large columns of gas prior to the 2012-13 eruptive sequence, while higher DSAR values prior to the 2016 eruption suggest that the subsurface was critically stressed.
At Whakaari/White Island, the 2013-2015 quiescence was associated with elevated DSAR values, while overall seismicity and gas emissions were particularly low. We interpret these observations as due to fluids and mineral precipitation accumulating in the subsurface. When pressurization was registered at fumarole F0 in 2015, the correlation between tremor and tides emerged for the first and only time, indicating that the volcano was critically stressed over long-term periods (months-to-years).
Finally, we could study for the first time different eruptive styles at the same volcano using noise-and tremorbased approaches. Similarly to phreatic eruptions at Ruapehu, Tongariro and Kawah Ijen , elevated DSAR values have been recorded prior to the 2016 phreatic eruption and between 2017 and 2018, but very low values were recorded prior to the 2012-2013 phreato-magmatic sequence. This observation stresses the need to multiply observations at different volcanic systems using stations deployed in the near field that provide continuous timeseries with a good signalto-noise ratio. This would be of the utmost importance to better understand the processes at play as well as their timescales.

Data and materials Chronology
This section compiles the overall observations related to volcanic activity during our period of analysis. These observations include lake level and deformation (Additional file 3: Figure S3a), fumarolic and lake temperatures (Additional file 3: Figure S3b), SO 2 (Additional file 3: Figure S3c) and CO 2 and CO 2 /SO 2 ratios (Additional file 3: Figure S3d) emissions, and seismic tremor amplitude (Additional file 3: Figure S3e and Fig. 3e) and the earthquakes (Additional file 4: Figure S4a-d) between 2007 and the end of 2018. When available, they are presented in this order for each period. For each eruptive period, we also provide a small summary presenting the key observations. At White Island, the deformation is typically derived from crater floor levelling surveys (Fournier and Chardot 2012;Peltier et al. 2009) and more recently via InSAR (Hamling 2017). Until 2015, we report elevations observed directly by observers on the main crater floor (Christenson et al. 2017). The percentage of lake level is then provided using the method described in Hamling (2017) using InSAR data.
The temperatures are either measured at fumarole F0 ( Fig. 1) or in the lake. SO 2 emissions are tracked with Cospec and using DOAS. CO 2 is collected using airborn platform (Werner et al. 2008). The terminology used to classify and analyse the earthquakes is presented in the Methods section.

Quiescent period from 2007 to 2009 (period ɪ)
From its highest ever recorded level (~ 1 m below the overflow; (Christenson et al. 2017), the lake had been declining until the start of our observation period in 2007, when it stood at ~ 5 m below overflow (Additional file 3: Figure S3a). Spring discharge on the main crater floor ceased by April 2007, and by September 2007, the lake was reduced to isolated pools in the lake basin (some 25 m below overflow). Spring discharge recommenced on the main crater floor in late 2007, leading to a new filling cycle and a rapid increase of the lake level to within 7 m of overflow by mid-2008 (Additional file 3: Figure S3a).
The temperature for fumarole zero (F0, blue triangles, Additional file 3: Figure S3b) was not available during this period. The lake temperature had been steadily decreasing (Lake, cyan triangles, Additional file 3: Figure S3b).
Cospec SO 2 emissions (magenta) varied between ~ 100 and ~ 400 t/d, with only rare excursions to about 600 t/d (Additional file 3: Figure S3c). CO 2 emission from the volcano, on the other hand, was pulsatory during the 2007-2009 period, ranging between 1000 and 3000 t/d (red circles, Additional file 3: Figure S3d), although the longer term trend points to slowly ramping emission of CO 2 during this interval (red line, Additional file 3: Figure S3d). CO 2 /SO 2 mole ratios decreased to ~ 5 in mid-2008 (green circles and line, Additional file 3: Figure S3d).
Seismicity during this period generally consisted of low amplitude tremor (RSAM, Additional file 3: Figure S3e), steady VT, RHF and VLP seismicity and slowly increasing LP earthquakes (Additional file 4: Figure S4).

Quiescent period from 2009 to August 2011 (period ɪɪ)
Water level in the lake was largely static over the mid-2008 to mid-2010 period (Additional file 3: Figure S3a). Despite this increased precipitation, lake levels declined between mid-2010 to mid-2011 (Additional file 3: Figure S3a). This pattern is consistent with uplift and was modelled as part of a shallow inflation source (Fournier and Chardot 2012), and later shown to include probable hydrostatic modulation by the crater lake level (Christenson et al. 2017).
Fumarolic temperatures at fumarole F0 ( Fig. 1) abruptly increased in late 2009 by some 40 °C to an all-time high of 218 °C (Additional file 3: Figure S3b), but slowly declined thereafter to a temporal low of 155 °C by early 2014. In contrast, the lake temperature increased (Additional file 3: Figure S3b). SO 2 emissions were relatively static, ranging from just ~ 100 to ~ 300 t/d, with a single emission rate of ~ 700 t/d measured in mid-2009 (Additional file 3: Figure S3c). By comparison, CO 2 emissions remained pulsatory (ranging between 800 and 3000 t/d), but the long-term trend continued to increase through to mid-2010 (Additional file 3: Figure S3d). CO 2 /SO 2 ratios were erratic (Additional file 3: Figure S3d).
Seismic tremor was generally very weak during the pre-eruptive period (Additional file 3: Figure S3e). The main identified discrete event types, especially discrete LP events, may be more evident during the low tremor period but are seen less commonly during later high tremor activity (Additional file 4: Figure S4). The period also included waning VT seismicity, and generally low numbers of RHF and VLP events.

Unrest to first eruption: 2011 to August 2012 (period ɪɪɪ)
Lake level slowly receded during this period (Additional file 3: Figure S3a), leaving two muddy pools by late July 2012. However, on 27 July, local tour operators noted a rapid rise in lake level of 3 to 5 m from the prior day's tours. This lake rise corresponded to a sharp increase in tremor (Additional file 3: Figure S3e) and ebullition in the lake, i.e., discharge of water from the aquifer beneath the lake into the lake basin. This ebullition persisted for several days and on 3 August. On 5 August 2012 04:54 UT, an explosive eruption was observed, the first significant event at Whakaari/White Island since 2000 (Chardot et al. 2015). Low level venting from this period persisted until about 13 August, after which low level steam and gas emissions were observed from vents within a small tephra cone, which had grown to become a small island within the shallow lake.
During this time, however, emission temperatures steadily declined, from 218 °C to 160 °C. This behavior continued through to latest 2012 (Additional file 3: Figure  S3b). SO 2 emissions slowly increased over values recorded during the 2007-11 period (Additional file 3: Figure  S3c), CO 2 emissions during this period continued to be pulsatory but did not exceed 2000 t/d (~ 30% lower than peak pulses measured during 2007-11 (Additional file 3: Figure S3d)), and CO 2 /SO 2 mole ratios were 10 or less (Additional file 3: Figure S3d). Importantly, from earliest 2012, F0 chemistry displayed strong departures from its previous behavior, wherein CH 4 /CO 2 ratios declined by one order of magnitude, and CO/CO 2 ratios increased by ~ 2.5 orders of magnitude (Additional file 10: Table S1).
In late July 2011, a sequence of mixed frequency earthquakes, including VLP (Additional file 4: Figure S4b), LP (Additional file 4: Figure S4d), and high frequency (Additional file 4: Figure S4a) components were observed ). These events were followed by a slow rise of RSAM (Additional file 3: Figure S3e, Fig. 3e) from background levels.

eruptive period and Dome emplacement: August to December (period IV)
Short description of the eruption here. Then: This period was remarkably quiet and the morphology of the crater and the lake were generally stable ). On November 24, tour operators observed a small irregular dome exposed within the tephra cone (Chardot et al. 2015;Jolly et al. 2020) but it was uncertain exactly when the dome was emplaced.
In the immediate post eruptive period, two small, shallow lakes occupied the eruption crater complex. No lake level and only two temperature measurements could be collected.
Post eruption RSAM levels reduced after the 5 August eruption, remaining low but still well above quiescent values (Additional file 3: Figure S3e). RSAM increased again rapidly on 25 August and persisted until 2 September when vigorous mud geyser activity, small explosions and ash venting were generated from the tuff cone vent. By 4 September, RSAM levels returned to stable but persistently elevated levels and the surrounding lake was generally very quiet.
In a retrospective analysis, Jolly et al. (2020) analyzed the persistent tremor recorded between September and December 2012, and observed very slowly evolving gliding spectral lines that started around the time of the late August ash venting period and persisted until after the dome emplacement.
The monitoring data was interpreted as reflecting a period of shallow magmatic degassing accompanying the dome emplacement ).

Dry lake unrest: January to July 2013 (period V)
Summary of the eruptive activity. From the time of dome emplacement, first observed on 24 November 2012, which provided a strong heat source, the main crater lake began to progressively evaporate, turning into a set of small isolated ponds by early January 2013. This activity led to the onset of discrete mud/molten sulphur eruptions (Christenson et al. 2017;Edwards et al. 2017;). This activity lasted until July 2013. The mud/ molten sulphur eruptions persisted intermittently from February to early April 2013.
As mentioned above, the lake nearly completely evaporated. Fumarolic data from F0 showed cycling between magmatic and hydrothermal vapour compositions (supplementary table), probably reflecting pressure transients associated with pulsatory transfer of magmatic vapour through the vent (Edwards et al. 2017).
This period was also marked by some of the strongest RSAM activity of the 2012-2016 eruption period, including bursts from mid-January to late February, mid-March to early April, and from late July to a steam eruption on 19 August 2013 (Additional file 3: Figure S3e).
The monitoring data may reflect a shift from gas emission being sourced from the conduit/vent towards the SE vent; this was possibly due to quenching/crystallization of the extrusion that reduced permeability in the western vent area.

Active eruption period: August-October 2013 (period VI)
No lake measurement is available for this time period. Renewed eruptive activity occurred 19 August when a small steam eruption generated an ash cloud reaching to ~ 4 km above the main crater rim. SO 2 emission around this eruption oscillated between 300 and 550 t/d (Additional file 3: Figure S3c), while CO 2 emission rates varied by a factor of two, with at least 2 days of strong CO 2 emission (> 2400 t/d on 20 August 2013 and 23 August 2013; Additional file 3: Figure S3d). CO 2 /SO 2 ratios oscillated between 6.4 and 10.6 (Additional file 3: Figure S3d) during this period. Interestingly, an exceptionally large SO 2 emission of 2075 t/d was observed on 20 August 2013 following the steam eruption. Deformation data were not available during this time period.
The 19 August eruption occurred in the evening (21:23 UTC) and was composed of steam ejection, and only minor mud fountaining within the mostly empty crater lake basin.
After episodes of strong tremor (Additional file 3: Figure S3e), three additional eruption sequences occurred on 4 October, 8 October and 12 October 2013 (Chardot et al. 2015). The 4 October eruption has been analyzed in detail by Caudron et al. (2018). This eruption was initiated by a Very Long Period (VLP) event located between 700 and 900 m depth below the crater, possibly triggered by the release of steam and gas that was visible from the mainland. Harmonic tremor was recorded ~ 25 min. Excepting SO 2 fluxes measured with DOAS increasing from 200 to 440 t/day (A. Mazot, unpublished data), very few data are available for this time period (Christenson et al. 2017). The final eruption of this sequence occurred on 12 October at ~ 20:09 Local Time.
The monitoring data probably reflect intermittent but ongoing periods of strong gas and tremor emission associated with the eruptive events.

Quiescent interlude: November 2013 to February 2016 (period Vɪɪ)
A crater lake was re-established in December 2013 and lake level surveys recommenced at that time (Additional file 3: Figure S3a). The last measurements available were collected early 2015 (cyan line, Additional file 3: Figure S3a). The first percentage of lake level measurements were available around the same time (blue line, Additional file 3: Figure S3a). The lake level progressively increased and by February 2016 covered the previously active vent(s).
Volatile emissions declined over the period between November 2013 to August 2015, with CO 2 levels falling to levels not seen since mid-2012 (Additional file 3: Figure S3d). However, a period of strong gas emission commenced in early October 2015 and lasted until 4 February 2016, with CO 2 emission rates ranging between 1800 t/d and ~ 2300 t/d. CO 2 /SO 2 ratios varied from 5 to 9.7 over this time (Additional file 3: Figure  S3d). This period was marked by low tremor amplitudes (Additional file 3: Figure S3e), a modest increase in number of discrete LP events (Additional file 4: Figure S4d) and low numbers of RHF, VLP and VT earthquakes (Additional file 4: Figure S4a-c). However, a swarm of LP and RHF earthquakes was recorded in mid-2015, coincident with the gas influx into the system at that time.

February to April 2016: eruption (period VIII)
Description of the eruption. In the evening of 27 April 2016, a multi-phase eruption occurred that consisted of at least six discrete explosive phases and a pyroclastic surge up and over the inner crater wall and out onto the main crater floor. Pyroclastic surges were likely fed by the collapsing jets ejected during each discrete phase; however, the inner crater wall likely formed an efficient barrier to all but the largest currents. The energy of each discrete explosion generally progressed from relatively weak to strong jetting over a 35-min period. A combination of seismo-acoustic analysis and ballistic ejecta mapping was used to infer that the eruption was sourced from several, partially inclined vents within the lake basin (Walsh et al. 2017;Jolly et al. 2018;Kilgour et al. 2019). Post-eruption observations and photogrammetry surveys found strong degassing and circular depressions distributed across the crater basin. At least one of these vents coincides with the SE vent; a frequently active vent that was the dominant one during the 2012-2013 eruption sequence.
The lake level stabilized in early March 2016, then abruptly declined just 6 days before an eruption on 27 April 2016 eruption (Additional file 3: Figure S3a).
Although we do not have lake temperatures at this time, the lake was vigorously steaming along its western margin, suggesting that evaporation was the major process responsible for the decrease of lake level. Interestingly, photographic records reveal that the major gas venting was limited to the western margin of the lake. Gas emissions and seismicity had been decreasing prior to the eruption (Additional files 3, 4: Figures S3c, d and  S4).
Pre-eruptive unrest indicators were generally lacking, and the eruption occurred with little effective warning. Nevertheless, a retrospective assessment noted the existence of discrete and extended VLP seismicity from about 2 h before the eruption ). This activity was located approximately 800-1000 m below the crater floor and may have been related to the progression of magmatic gases towards the surface.

Post eruption quiescence: April 2016 to December 2018 (period IX)
Shallow ponds were present in the eruption crater basin through to the start of 2018, when a proper lake started to accumulate again. By the end of 2018, the lake level had increased (Additional file 3: Figure S3a). The availability of InSAR data provided insights into the deformation. Deformation was largely concentrated near the crater lake and adjacent south-west crater wall (Hamling 2017). Downslope displacements of the crater wall exceeded 200 mm/year in 12 months after the eruption (Hamling 2017), but stabilized with the gradual increase in lake level. At the same time, areas at the back of the crater lake and much of the crater lake basin, which had previously been submerged, showed signs of rapid subsidence at rates of ~ 150 mm/year (Hamling 2017). Through 2017, the rate of subsidence gradually slowed, and by mid-2017 had switched to uplift (Hamling 2020). The change from subsidence to uplift coincided with the refilling of the crater lake through 2018 (Additional file 3: Figure S3a) and observations of renewed uplift within the main crater floor, in the vicinity of Donald Mound. Line-of-sight displacements indicated uplift at rates of ∼50 mm/year consistent with shallow inflation beneath the lake floor (Hamling 2020).