Relating gas ascent to eruption triggering for the April 27, 2016, White Island (Whakaari), New Zealand eruption sequence

Abstract The April 27, 2016 eruption sequence at White Island was comprised of 6 discrete eruptive events that occurred over a 35-min period. Seismicity included three episodes of VLP activity: the first occurring ~ 2 h and a second occurring 10 min prior to the first eruption. A third larger VLP event occurred just prior to the fourth eruption. A VLP source depth of 800–1000 m below the vent is obtained from an analysis of the waveform semblance, and a volumetric source is obtained from waveform inversion of the largest VLP event. Lag times between VLP occurrence and eruption onsets provide an opportunity to examine gas migration and stress transfer models as potential triggers to the eruptive activity. Plausible lag times for a deep gas pulse to the surface are obtained by application of a TOUGH2 computational model which suggests propagation times of 0.25–1.9 m/s and are informed by previously measured White Island rock porosities and permeabilities. Results suggest that pre-eruption VLP may be plausibly linked to advection of gas from the VLP source at a magmatic carapace located ~ 800–1000 m depth. Alternatively, the large VLP that occurred just prior to the fourth eruption may be linked to a quasi-dynamic or quasi-static stress perturbation. Graphical abstractXXXX

section_content: White Island (Fig. 1 ) displayed two periods of volcanic unrest from July 2012 to April 2016 separated by a quiescent period from October 2013 to October 2015. The end of the second unrest period was marked by a sequence of 6 discrete eruptions beginning April 27, 2016, at ~ 21:37 NZST (9:37 UTC and all times hereafter in UTC) with the last discrete burst occurring about 10:12 (Fig. 2 ). The duration of the eruption sequence was about 35 min and the inter-event periods included elevated persistent seismic tremor and spasmodic bursts but no distinct acoustic activity. The eruption deposits included mud, ash and aggregate hydrothermal materials from the crater lake but did not include juvenile magmatic material (GeoNet Volcanic Alert Bulletin, 2016/14: https ://www. geone t. org. nz/vabs/4jWBX nKmn6 sWMs6 KUKGi aK). Hence, the eruption may be regarded as a purely phreatic sequence. In addition to the observed eruption record, VLP earthquakes were recorded on a small broadband seismic array and may represent the flux of magmatic fluids (gas and liquid without magma melt) to the surface. The timing between observed VLP earthquakes and the surface eruptive activity provides a unique opportunity to examine possible links between subsurface VLP events to the eruptive activity and to constrain these relationships to the range of possible eruption trigger mechanisms. 

White Island (Whakaari in Maori) is an andesitic composite cone volcano located in the south-eastern Bay of Plenty (Fig. 1 ). The horseshoe-shaped eastward opening crater has an active central subcrater whose base was near sea level at the time of the eruption. The greater crater system has persistent super-heated fumaroles (Werner et al. 2008), elevated temperature springs (Christenson et al. 2017 ) and an active central crater characterized by highly variable lake levels, temperatures and lake ebullition, which testify to the existence of a long-term (decades to centuries) active hydrothermal system. The active crater has a lake (Fig. 1 ) whose level has varied from isolated shallow pools to filled to the upper scarp walls over the past 15 years, a result of the strong interplay between the ingress of sea and meteoric water from above and juvenile water, gas and heat from below. Persistent long-term degassing at White Island requires extant magma in the upper km or two of the vent system. This might occur via the pervasive injection and crystallization of magma to shallow depth (e. g., Cole et al. 2000) or alternatively via the presence of a shallow convecting and degassing magma system (e. g., Kazahaya et al. 1994; Stevenson and Blake 1998). For the former, it is uncertain if the required magma supply rates would yield sufficient degassing without strong deformation. For the latter, it is uncertain if degassed and viscous magmas can convectively resupply fresh magma to the hydrothermal system. 

As a result of the shallow magma supply and active hydrothermal system, the volcano has also produced common volcano tectonic (VT) (Nishi et al. 1996), long-period (LP) (Sherburn et al. 1998 ) earthquakes and tremor (Sherburn et al. 1996) over several decades of monitoring. The installation of broadband seismometers after 2007 allowed discrimination of very-long-period (VLP) earthquakes which may be generated at a relatively shallow depth of ~ 1 km below sea level (Jolly et al. 2017). The VLP source is thought to represent intermittent failure of a magma carapace which may act as a leaky obstruction which traps and releases magmatic gases that subsequently propagate into the overlying hydrothermal system (Jolly et al. 2017). Infrasound monitoring was first established at White Island in November 2002 and has been applied to confirm eruptive activity ranging from larger phreatomagmatic eruptions, to small repetitive phreatic activity from a shallow mud/molten sulfur pool (Jolly, et al. 2016; Edwards et al. 2017; Schmid et al. 2017; Christenson et al. 2017). White Island has also produced explosive magmatic, and strombolian activity (e. g., Houghton and Nairn 1991; Cole et al. 2000) as well as a dome forming eruption (Chardot et al. 2015) The primary goals of this paper include the constraint of subsurface VLP source location/mechanism and to establish the range of plausible links to the surface eruption activity. To these ends, we first examined the timing of VLP occurrence and its relationship with surface eruptive activity. Next, we determine its location using a waveform semblance approach, and examine the source properties via a waveform inversion. Finally, we modeled the lag-time predicted from TOUGH2 models (Pruess 1991) for a gas pulse propagating from the known VLP source region to the surface. The modeling allowed us to examine two plausible end-member eruption triggers, namely gas advection (mass transfer) from depth to the surface, or stress transfer from either a static stress pulse or via the dynamic stresses from elastic waves. 

Observations suggest that the April 2016 eruptions may be related to both mechanisms. While acknowledging the difficulty in rigorously linking subsurface earthquake activity to any surface observation, the methods outlined here provide a framework for identification and interpretation of pre-eruption subsurface activity, and how they relate to subsequent eruptions. This information can inform unrest and eruption monitoring in volcanic environments.
section_name: Introduction

section_content: The permanent volcano monitoring network included two permanent GeoNet seismic stations (WIZ and WSRZ) during the eruptive period (Fig. 1 ). In addition, GNS Science operated four portable seismic stations (WI01, WI02, WI04, and WI13) as part of the re-occupation of a wider portable network from 2011 (Jolly et al. 2017). The GeoNet seismic stations consisted of Guralp 3ESP broadband seismometers with Quantrerra Q330 digitizers. In addition, GeoNet operated Setra model 270-600 and InfraBSU pressure transducers located near the permanent seismic stations. The portable seismometers were composed of both Trillium compact seismometers and 24-bit Nanometrics Taurus digitizers. All seismic data are sampled at 100 Hz and have a corner frequency of 120 s. 

The eruptive activity consisted of 6 small discrete eruptions which occurred at about 9:37, 9:43, 9:48, 9:54 10:05 and 10:12 UTC based on the acoustic record (Fig. 2b ). Subsequent evaluation of the eruptive deposits suggest that the event was phreatic in nature and had pre-eruptive shallow vent conditions consistent with localized hydrothermal mineralogical sealing over numerous lake vent fumaroles (Christenson et al. 2017). A large VLP earthquake was observed within the eruption sequence (Fig. 2a, c ) and preceded the fourth eruption by about 40 s. Analysis of the map and section particle motions for the onset of this event (Fig. 3 ) show strong rectilinear motions that point to a source position near the volcanic crater and at a depth similar to earthquakes observed in 2011 (Jolly et al. 2017). Prior to the eruption, we observed a sequence of VLP earthquakes that had low signal-tonoise levels and occurred around 7:35 UTC and between 9:25 and 9:35 UTC (Fig. 4 ).
section_name: Data and observations

section_content: We used the waveform semblance as described in Kawakatsu et al. (2000) for the analysis of the VLP activity. The eruption sequence was recorded on a subset of seismic stations that were also utilized for a VLP swarm occurring in August 2011 (Jolly et al. 2017). From Kawakatsu et al. (2000), coherency among seismograms can be characterized using:

where L is the length of the seismogram, N is the number of stations and u i,j(i) ≡ u i (t i + jΔt) is a seismogram of the ith station at the j(i)th time sample at the start time t i. 

For a stationary isotropically radiating point-source, seismic waves will propagate as compressive body waves that may be rectified into a linear particle motion (Kawakatsu et al. 2000). In the radial component (the direction of a station from the source), the seismogram particle motion is given by u i,j(i) ≡ R i,j(i). Differences in the multi-station travel-times are accommodated with adjustments for these time delays (t i ), which should improve the radial waveform coherence. The mutually perpendicular components H i,j(i) and V i,j(i) for an isotropic radiating source, should have negligible contributions. In this case, the waveform semblance (Kawakatsu et al. 2000) may be applied to emphasize the radial excitation from the isotropic source and penalize departures from rectilinear behavior via:

where V is the component in the direction perpendicular to R within the vertical plane, which contains both source and receiver, and H is in the horizontal component perpendicular to both R and V (V may not be vertical). L is the number of time samples and N is the number of stations. We take the scaling factor D as:

(2)

and S 3 measures the rectilinearity of particle motion pointing via the signal coherency in the radial direction. 

As in Kawakatsu et al. (2000), we normalize each seismogram so that the RMS amplitude of each signal at each station becomes unity:

which removes the amplitude information from the original seismograms giving equal weight to each station during the coherence computation. We performed the semblance analysis in a moving time window approach with each time step including a grid search for the semblance from volcanic VLP activity. Because the semblance utilizes a waveform normalization, we present the data with only minimal pre-processing, establishing robust time synchronization, and then filtering the data within the VLP passband of 0. 03-0. 1 Hz (Fig. 2d ). Note that although the corner frequency of the instrument response is 0. 083 Hz (120 s), the frequency band between the corner frequency and 0. 03 Hz was heavily contaminated by the lateral component noise and had to be filtered out. The upper limit of the filter was chosen in order to isolate the coherent VLP waveform across the network. The semblance is computed within 1-min non-overlapping time steps with an initial coarse grid analysis. The coarse grid had spacing of 500 m over a cube of 3 × 3 × 3 km centered on the eruption vent. 

After the coarse grid search, the semblance was recomputed for high semblance (S 3 > 0. 6) time steps yielding precise locations estimates for the pre-and syn-eruptive VLP. To obtain location error estimates, we used methods outlined in Almendros and Chouet (2003) who established a criterion semblance level S L = (1 -δS)S MAX, based on the signal-to-noise ratio (SNR) where δS = 0. 062 SNR -1. 54. In our case, we measured the noise for a 5-min period prior to the first VLP event (Fig. 2a ) and utilized this to obtain error bounds from the semblance level S L. The dense 3D grid search was completed over a 2. 0 km grid centered on the crater lake and from sea level to 2. 0 km below sea level. The grid spacing was 25 m. The grid search encompassed the full local network, and hence offered good resolution of potential activity related to the eruption. At each grid point, we compute the source to station distance, azimuth and incidence angle as well as rotate the E, N, Z seismograms for each station into R, H and V components having an appropriate travel-time-delay based on the independently determined 1D velocity model (Jolly et al. 2017). Straight line ray-paths are assumed for rotated components which is probably valid for all stations on consolidated volcanic rocks (WI01, WI02, WI04, WIZ, WSRZ) based on prior active source experiments (e. g., Jolly et al. 2012). For each time step the waveforms are normalized using Eq. ( 5 ) and the waveform semblance is computed with Eq. ( 6 ). The datum for all subsequent analysis is sea-level in New Zealand Map Grid (NZMG) coordinates and are expressed in meters. 

The waveform semblance computation used in this paper was rigorously tested using both synthetic waveform data and a swarm of 25 VLP events that occurred at White Island in August 2011 (Jolly et al. 2017). The earlier synthetic tests incorporated uncertainties in velocity structure and source geometry and monitored departures between the computed and retrieved locations. The synthetic tests allowed uncertainties for the 2011 swarm to be robustly assessed and were about ± 0. 2 km both laterally and in depth for locations computed on 12 broadband stations. Given the sparse 6 station array available for the 2016 eruption, we re-performed the prior synthetic analysis for the case of a VLP earthquake having an isotropic source excitation at 1 km depth and noting that one station (WI13) was not available in 2011 and hence, had no synthetic computation. The revised test shows that the synthetically derived location was recovered to within a few hundred meters by our waveform semblance algorithm (Table 1 ), similar to location results obtained for the earlier dense array synthetic relocation. The results suggest that the sparse network had sufficient spatial distribution and azimuthal coverage to recover VLP locations within the target region. 

Results of the present semblance analysis show that locations are clustered from 600 to 1000 m depth (Fig. 4 ) and that the largest VLP (Fig. 4a and b ) have the strongest semblance (Fig. 4c ) and smallest error bounds (Fig. 4d ). This was also observed in the earlier 2011 VLP swarm (Jolly et al. 2017 ) and suggests that the shallower VLP solutions (blue symbols in Fig. 4d ) may indeed suffer from strong noise overprints that cause systematic mislocations to shallower depths. If so, then the full VLP sequence may represent a source excitation at a depth of (800-1000 m) (red symbols in Fig. 4d ). The VLP
section_name: Waveform semblance

section_content: The synthetic waveforms are computed at a position 1 km below sea level. The recovered solution is within 200 m of the true source position and is similar to locations derived from 12 stations available in 2011 (see Jolly et al. 2017). The synthetic location is determined using the same method as outlined in the methods section, 3 ) is well located at ~ 1 km depth beneath the southern crater wall (Fig. 4 ). It lies within the same location distribution found for the 2011 earthquakes (Jolly et al. 2017 ) and has similar location uncertainties (Fig. 5 and Table 2 ). Given strong coherence of waveforms (Fig. 3 ) reasonable signalto-noise (Fig. 4 ), and location (Fig. 5 ) for the sun-eruption VLP, we proceed next to an analysis of the moment tensor and source-time function for this specific event.
section_name: Table 1 Summary of semblance location for synthetic waveforms from an isotropic source

section_content: Like in the previous study of the VLP events from 2011 (Jolly et al. 2017), we analyzed the source mechanism of the August 2016 VLP signals using a full waveform moment tensor inversion. The Green's functions (GF's) were calculated with the Spectral Element Method (SEM) code EFISPEC3D (http://efisp ec. free. fr; De Martin 2011), using a numerical domain based on a 10 × 10 × 7 km 3 rectangular cuboid with a free surface on top and absorbing boundaries at the bottom and the sides of the cuboid. 

The free surface was then altered combining a Digital Elevation Model (DEM) and a bathymetry model of the surrounding seabed (Prasetya and Wang 2011), in order to reflect the topography of the volcanic edifice (for further details see Jolly et al. 2017 ). The source location obtained in the previous section was used for the computations. The velocities in this model were homogeneous (V P = 2. 2 km/s, V S = 1. 272 km/s and the density ρ = 2120 kg/m 3 )-an assumption justified by the long VLP wavelengths in comparison with experiment dimensions. As tilt ground motion can have a strong effect on long-period seismograms (e. g., Rodgers 1968), we also included tilt in the GFs, following the work of Maeda and Takeo (2011) and Van Driel et al. (2015). The gravitational acceleration affecting the i-th horizontal component can be expressed using the spatial derivatives of the ground motion around the station positions:

(5) where g is the gravitational acceleration, and u i are the horizontal translational displacements. The spatial derivatives in Eq. ( 5 ) were calculated in the wavenumber domain directly in the EFISPEC3D code, rather than being approximated by a central differences method as in Maeda and Takeo (2011). The obtained tilt traces ( a i were numerically double integrated with respect to time and added to the corresponding displacement GFs, thus obtaining the 'apparent' GFs. In the frequency domain, it can be written as:

where GF i (ω) denotes the i-th horizontal component of the apparent Green's function, GF i (ω) is the translational Green's function, g is the acceleration of gravity, and ω is the angular frequency. As the tilt contribution (the second term in Eq. 6) in the near field of the source leads to a large static step at the end of the time window for which a GF is calculated (see Van Driel et al. (2015) for details), the inversion in frequency domain becomes unstable. In order to remove this step and stabilize results, Maeda and Takeo (2011) suggest the convolution of the GFs with the instrument response, while Van Driel et al. (2015) perform double temporal differentiation of the data and GFs prior to the inversion. Instead, we perform the inversion of band-pass-filtered waveforms, where possible tiltinduced trends/steps are removed by band-pass-filtering, as shown in the next section. It is also important to mention that tilt can significantly affect the seismic records, (6

even when narrowly filtered within the VLP spectral band. A good illustration of this is given in Fig. 5a of Van Driel et al. (2015), where the influence of the tilt is clearly seen at some stations for the frequencies 0. 03 and 0. 1 Hz.
section_name: Moment tensor inversion

section_content: After the instrument response was removed, the observed waveforms were integrated to displacement. 

While the vertical component of motion showed a clear static offset related to VLP waveform, the horizontal components were dominated by tilt and long-period noise. This is illustrated by Fig. 6, showing the unfiltered instrument-corrected displacement waveforms at station WI04. Note that ultra-long-period signal on the horizontal components makes it impossible to recover the static offset on horizontal components. Hence, we band-passfiltered the horizontal components between 0. 03 Hz and 0. 1 Hz. In order to preserve the static offset on the vertical components, we low-pass-filtered them below 0. 1 Hz and applied a 100-s-long median filter in order to stabilize the observed displacement step (see Thun et al. 2015 Thun et al., 2016)). The filtered displacement waveforms are shown in Fig. 7. 

The inversion was carried out in two stages: in the first stage we determine the source mechanism by performing the standard full waveform source inversion of the all components of the band-pass-filtered velocity waveforms (0. 03 Hz-0. 1 Hz) in the time-domain; in the second stage, we invert for the broadband source-time function (STF) by using the low-pass-filtered vertical components (shown in the last row of Fig. 7 ), and fixing the source mechanism to the one obtained above. In this way, we utilized the whole bandwidth of the vertical waveforms, even at very low frequencies (< 0. 03 Hz) where the horizontal signals were unusable. The procedure is similar to the work of Maeda et al. (2015a), where the authors invert the band-pass-filtered waveforms first, and then use the recovered STF to derive the "true" STF by inverse filtering. However, in their procedure, the assumption about the shape of the true STF is needed, while in this study it is obtained directly from the vertical components of the data. 

Single forces were not included in the inversion, just like in Maeda and Takeo (2011). We were inverting for the simplest possible source model which explains the data, and achieved excellent results without the need for single forces (see next section). In addition, we wanted to avoid spurious single forces, often arising due to the source mislocation and velocity mismodeling (De Barros et al. 2013). 

The misfit between the data and synthetics associated with the solution was calculated by the following equation: where d s,c t and s s,c t denote the data and synthetic time samples, respectively, recorded at instant t, at the station s and component c.
section_name: Data and inversion strategy

section_content: The observed waveforms for the VLP recorded just prior to eruption P4 (Fig. 2 ) were inverted, showing the misfit between the real and synthetic data of R = 0. 07. The obtained fit of velocity waveforms is shown in Fig. 8. It shows the excellent match between all the waveforms, except N component of station WI02. Close inspection for this station/component reveals some discrepancy compared to all other waveforms across the network. The reason is likely due to the presence of a strong tilt component on the waveform, which cannot be accurately modeled due to the source location error. Another possibility 2017 ), the moment tensor obtained here cannot be explained by a single tensile-shear crack: the inversion for the orientation of the crack, the direction of the crack opening and the Poisson's ratio (see Vavryčuk 2001) results in the non-physical values of the Poisson's ratio for any mode of a tensile crack opening (from a pure shear to a pure tensile crack). Instead, the solution can be explained by a non-planar source with the main tension axis in the NNE direction. 

The moment tensor from Eq. ( 8 ) can be decomposed into an isotropic (ISO), a compensated linear vector dipole (CLVD) and a double-couple (DC) component (e. g., Vavryčuk (2015) ):

where M = 2. 7 • 10 13 Nm denotes the seismic moment, and C ISO = 52%, C CLVD = 32% and C DC = 16% are percentages of the ISO, CLVD and DC components, respectively. It is important to mention that only the ISO component is uniquely determined from the trace of the matrix M, while the rest-called the deviatoric part of moment tensor-can be decomposed in an infinite number of ways. As a consequence, many different conceptual models could relate to the same moment tensor. Thus, the moment tensor obtained here could also be attributed to the source involving expansion and fluid movement between ellipsoidal cavities (see Mizuno et al. 2015).
section_name: Source inversion results

section_content: The recovered source-time function (STF) is shown in Fig. 9. It consists of a waveform, followed by a static offset. It is important to mention that the waveforms recorded across the network are directly proportional (10) to the STF, that is, they are mainly a combination of the near-field and intermediate-field wavefield terms (see Lokmer and Bean 2010). In other words, the network is sitting on the source and directly sensing its movement, while the source radiation is a minor contribution to the VLP signals. 

In order to determine the volume change in the source, we extract the isotropic part of the moment tensor:

Note that the value of the isotropic seismic moment, M ISO 0 = 1. 43 • 10 13 Nm, corresponds to the maxi- mum of the STF shown in Fig. 9. In order to estimate the permanent volume change, we set the value of the seismic moment to 40% of its peak value (see the value of permanent offset in Fig. 9 ), i. source types between these two end members. Hence, we use both Eqs. ( 12 ) and ( 13 ) in order to estimate the volume change in the source. We estimate the elastic parameters in the source region from a P-wave velocity ranging between 1. 3 and 2. 3 km/s, a density of 2160 kg/ m 3, and a Poisson's ratio between 0. 25 and 0. 35. For such a choice of parameters, and using Eq. ( 12 ), the estimated permanent volume change is ranging between 374 and 1560 m 3. For the same choice of parameters and using Eq. ( 13 ), the permanent volume change is between 560 and 2340 m 3. On the other hand, the total volume of material potentially fluxing through the system could be determined from the peak volume change which is 2. 5 times larger. Using similar approaches, Maeda et al. (2015b) found a volume change of 75 m 3 for the 2014 Ontake eruption sequence, while Ohminato et al. (1998) found that the modeled volume change from Kilauea was about 3500 m 3. Hence, the VLP event analyzed here can be regarded as an upper bound VLP source as expressed by volume change. Its moment magnitude is equivalent to a M W = 2. 9 earthquake. 

Extending the waveform inversion of the largest VLP to the earlier activity was precluded by low signal-to-noise levels, especially on lateral component seismograms and the limited number of seismic stations, which do not allow for robust source modeling. Therefore, it is difficult to draw conclusions about the source process of the preeruption activity or to state if they result from the same general source excitation.
section_name: Source-time function and volumetric change in the source

section_content: As stated in the data and semblance sections, we observed a sequence of small pre-eruption VLP earthquakes that occurred in the hours prior to the onset of the first eruption (P1) (Fig. 4a, b ). An examination of the time series using the waveform semblance allows separation of the pre-eruption VLP from non-volcanic seismic noise (Fig. 4c, d ). The first VLP occurs about 2 h prior to the first eruption at ~ 7:36 UTC. Additional VLP occurs just prior to the first eruption between 9:28 and 9:35 UTC (Fig. 8a ) and is summarized in Table 2. 

While the lack of juvenile materials in the eruption deposits strongly point to the link between VLP excitation and gas discharge, it does not exclude the possibility that gas charged magmas, and not just the exsolved gases, propagated at depth. Indeed, the longer-term offset observed on the vertical component seismograms attest to a substantial 'permanent' deformation just prior to eruption P4. While the lack of magma within deposits demonstrate that magma likely did not propagate deeply into the overlying hydrothermal system, we cannot exclude the possibility that the VLP source represents a rapid 'propping' of the magma carapace which in turn promoted gas slug migration. 

The largest VLP event initiates at 2016-04-27 09:53 UTC and eruption (P4) is observed about 40 s later. If we contend that the P4 eruption was caused by advection of gas from the VLP onset ~ 40 s prior and 1000 m deeper, then the propagation velocity would be ~ 25 m/s which seems implausible. Alternatively, the two events may be completely unrelated, may be related to the earlier VLP earthquakes or may relate via transit of a stress response to the system. This example illustrates the importance of constrained advective gas transport models as a tool to interpret the link between subsurface VLP and the trigger mechanism for eruption activity. To this end, we next determine reasonable gas pulse rise times specific to White Island conduit conditions.
section_name: VLP as a precursor to eruption

section_content: To examine the importance of advective mass transfer at White Island, we have undertaken an integrated finite difference simulation of the vent/lake environment using TOUGH2 (Pruess 1991). The model has radial symmetry with 2400 reservoir elements and a single atmospheric block of large volume (10 50 m 3 ). The mesh consists of 80 elements of 10 m length along the z axis (800 m total depth), and 30 elements along the x axis with a scaling factor of 1. 14 and radial distance of 300 m (2400 elements). This numerical model is consistent with a VLP source being constrained at a ductile magma carapace located between the magma conduit and the shallower hydrothermal system (Jolly et al. 2017). The carapace would act as a gas trap, and failure of the trap (the VLP) could be driven by over-pressurization at the top of the magma column (Fig. 11 ). Subsequent migration of the gas pulse through the overlying hydrothermal system would be driven by differential pressure of the gas pulse compared to the surrounding medium at its equilibrium conditions and their attendant permeabilities and porosity. 

The crater rocks have model porosities of 0. 30, and permeabilities of 1 -13, 1 -13 and 1 -11 m 2 in the x, y and z directions, respectively. These values are in line with measured porosity (0. 01-0. 7) and permeability (~ 10 -19 to ~ 10 -11 m 2 ) values at White Island (Heap et al. 2017). A higher permeability along the z direction was chosen to simulate the effects of steep concentric fractures which are likely present in eruptive crater environments (Gundmundsson et al. 1997). The center four elements of the radial grid have a combined horizontal (i. e., x-y) surface area of 54 m 2, and represent the conduit with a diameter of ca. 8 m. These elements were assigned a porosity of 0. 35, and permeabilities of 10 -12 m 2 in the x and y directions. We modeled three scenarios, varying only the vertical permeability parameter between the 10 -11 and 10 -9 m 2 in order to ascertain whether permeability has an appreciable effect on CO 2 transport through the conduit. 

The equation of state module used here (EOS2) solves state properties for H 2 O and CO 2 mixtures. Each scenario was initially run to steady state at 25 °C to establish cold hydrostatic pressure conditions in the model aquifer. A guiding assumption for the modeling is that the magma conduit is hot enough to sustain at least 2-phase vapor-liquid or single-phase vapor conditions in the column overlying the magma carapace (cf. conceptual model of the vent system of Christenson et al. (2017) ). To achieve these conduit conditions, we inject sufficient high enthalpy water vapor (3 × 10 6 J/kg) to create either single-phase vapor, or 2-phase vapor-liquid conditions in the conduit. Minimum vapor saturations of 0. 5 are here deemed appropriate to readily facilitate vapor phase transport of the CO 2. 

Into these pre-conditioned conduits in the three scenarios, we introduce a pulse of CO 2 (4 kg/s of both CO 2 and H 2 O) into each of the four inlet elements, which corresponds to the representative emission rate of ca. 1380 T/d CO 2. CO 2 arrival times (s) into the inner-most lake element are plotted as functions of permeability in Fig. 10. The introduction of gas into the four inlet elements may be regarded as the onset of a VLP transient. 

The modeling shows that gas transfer rates increase coherently with permeability in environments having high levels of gas saturation. For the lowest permeability case (10 -11 m 2 ), the model CO 2 transfer time of 28,000 s equates to a velocity of ~ 0. 03 m/s, whereas the rate for the highest permeability considered here (10 -9 m 2 ) is ~ 1. 9 m/s. 

The permeability characteristics of the actual conduit which connects the inferred magma chamber at 800 m depth to the crater floor of White Island are unknown. Although the maximum measured permeability for White Island lithologies is ~ 10 -11 m 2 (Heap et al. 2017), noting the dynamic nature of the conduit environment which consistently delivers between 1000 and 3000 T/d of CO 2 to the surface (Christenson et al. 2017), the bulk vent permeabilities are likely to be higher than this maximum measured value. A permeability of 10 -10 m 2 yields a transfer velocity of ~ 0. 25 m/s, which is comparable to the crack propagation speed determined from migrating VT earthquakes prior to the 2014 Ontake eruption (Kato et al. 2015 ). In addition, numerical simulations of VLP deformations associated with gas slug ascent (O'Brien and Bean 2008) yield ascent times of 0. 3-0. 8 m/s. While the ascent of gas slugs is dynamically different from the gas pulse migration modeled here, the alternative modeling yields results that are in the lower range of gas ascent times from our modeling. Whereas no VT events were recorded immediately before the 2016 White Island eruption, it is likely that gas transfer was less hindered, and probably travelled at a faster rate than that calculated for Ontake. 

It is worth noting that it proved impossible to maintain vapor-static conditions in the model conduit considered here with permeabilities higher than 10 -9 m 2. Above this value, vapor saturations fall to between 0. 1 and 0. 3 which precludes independent transport of a vapor phase, and resulting transport velocities which are slower than observed for permeability of 10 -9 m 2. This suggests that CO 2 transfer in the conduit probably falls within the range of 0. 25-1. 9 m/s. For White Island VLP earthquakes at ~ 800-1000 m depth, this would imply surface expression lag times between ~ 7 min to ~ 1 h for plausible White Island permeabilities (Fig. 11 ). If lower permeabilities exist, then time lags of several hours may also be plausible.
section_name: Transfer rates of CO 2 pulses through the conduit

section_content: The mass flux times obtained by modeling using realistically constrained vent permeabilities allow us to examine the subsurface VLP events in context with eruptive activity. The first pre-eruption VLP occurs approximately 2 h prior to the first eruption. Given the source depth of 800-1000 m, the VLP may not be related to the advective transfer of gas from the inferred magma supply zone (e. g., Christenson et al. 2017; Jolly et al. 2017 ) into the quasi-stable shallow hydrothermal system (Caudron et al. 9 Normalized time-dependent seismic moment. Note that it is not possible to say if the static offset remains at the 40% of the peak amplitude or slowly decreases toward zero at later times, thus forming a sawtooth function 2018). The next burst of VLP is at a similar source depth, begins about 10 min prior to the first eruption (P1) and stops just before the eruption onset (Figs. 4, 5). This VLP tremor seems plausibly related to the first eruption via an advection mechanism and may also relate to later eruptive events (P2-P6). The large VLP event occurring just prior to the P4 eruption is unlikely to be linked to the P4 eruption via advection, and hence, other mechanisms for the eruption should be examined. 

We consider next the concepts of a static stress change as an eruption trigger mechanism as well as dynamic triggering from passing seismic waves. The former is regarded as a primary trigger of near rupture aftershock sequences from large earthquakes (e. g., Gomberg et al. 2003), while for the latter, the passage of surface waves from large earthquakes are noted to perturb seismicity in geothermal and volcanic environments (e. g., Hill et al. 1993 Hill et al., 1995) ) at great distances. In our case, dynamic versus static stress perturbations would be difficult to disentangle because we are dealing with a triggered site (the eruption vent) that is very close to the VLP source region. For the large VLP earthquake, there is very little far-field radiation from the source and hence, the entire waveform can be regarded as a deformation pulse including both a near-field oscillatory transient and an apparent static offset (see Fig. 7 noting that the static offset is measured from a 120 s seismic instrument and may also represent an unrecovered sawtooth-shaped deformation pulse having a longer duration than the sensor corner frequency). While surface wave excitations from large earthquakes at teleseismic distances (e. g., Peng et al. 2018 ) may trigger a near surface response in a quasi-dynamic manner, it is uncertain if near-field VLP may generate dynamic stresses in the same manner. Alternatively, the eruption trigger could be due to a quasi-static near-field deformation from either the oscillatory VLP or the apparent static offset (Fig. 7 ). In this case, the perturbation and failure of a shallow hydrothermal seal (by either advective pressurization or a stress perturbation) could induce localized failure and explosive ejection of the shallow seal materials (e. g., Caudron et al. 2018). In these cases, the transfer of stress could propagate at elastic wave velocities and time lags might be very short, while an advective transfer mechanism may allow much longer eruption lag times (Fig. 11 ). We here regard both options as plausible eruption trigger mechanisms at White Island. 

While all near-field ground shaking (tremor and shaking from prior eruptions) may have triggered the April 2016 eruption sequence, the modeled VLP just prior to eruption P4 would be the only event where quasidynamic seismic waves and a quasi-static deformation pulse are temporally coincident to the eruptive activity (Fig. 7 ). Hence, we regard the P4 eruption as plausibly related to a stress change associated with VLP activity but emphasize that the eruption may also be driven by advection from earlier VLP or could be completely unrelated to the subsurface observations. Ideally, a more detailed analysis of the VLP source mechanism for pre-eruption activity could be applied to further constrain the trigger for each specific event in the eruption sequence. As stated before, the pre-eruptions VLP are only weakly recorded and their waveform characteristics are obscured, especially for their lateral components. This, in turn, impacts the source location and waveform features, making detailed earthquake source mechanism assessments problematic for the earlier preeruption seismicity.
section_name: Triggering mechanisms for volcanic eruptions

section_content: We have conducted a detailed examination of the preeruption and eruption seismic sequence, obtaining a record of VLP excitations. For the largest VLP earthquake, we obtain the best fit location at ~ 1000 m depth from a waveform semblance analysis. We invert the waveforms for the source-time function and source Fig. 10 CO 2 travel-times plotted as functions of permeability. The gas pulse originates at 800 m depth which marks the VLP source location which is inferred to be the boundary between the magma supply zone and the hydrothermal system geometry and recover an excitation of a mostly volumetric source. Analysis of the vertical component records via a median filter approach shows that the VLP contributes to an ~ 100 micron uplift on all stations consistent with a small inflation at the source region. 

The VLP source process (established via well-determined location, source mechanism and deformation analysis) may be related to surface eruption events by assessing the observed earthquake to eruption lag times. These lag times may then be compared to advection and stress transfer models as an aide to interpretation. Future work may bring opportunities to evaluate longer time series and establish links via robust statistical tests. The approach utilized here may be generally applicable to other volcanic systems worldwide. 

Fig. 11 Cartoon depicting the inferred volcanic conduit system at White Island. Magma is thought to reside at a depth of approximately 1 km depth and provides the source of gas and heat to the shallow hydrothermal system. These are transferred to the surface as a single-phase gas (aquamarine), a two-phase vapor and liquid (tan) and single-phase liquid (blue). The region surrounding the VLP source comprises a ductile low permeability zone (brown) which may intermittently rupture causing high gas discharge and VLP activity. The inset box depicts two potential mechanisms whereby subsurface activity from VLP might be linked to surface activity. Stress transfer might occur almost instantaneously, while mass advection from depth is inherently slower and is dependent on the balance of mass driving and resisting forces
section_name: Conclusions

section_name: Acknowledgements

section_content: This study incorporated GeoNet data. The research was supported by New Zealand Strategic Science Investment Funding (SSIF) from the New Zealand Ministry of Business and Innovation (MBIE) to AJ, IL and BC. Geoff Kilgour and Jerome Salichon are thanked for reviews of the manuscript prior to journal submission. Three anonymous journal reviewers and editor Stephanie Prejean are thanked for their efforts to improve the manuscript.
section_name: Acknowledgements

section_name: Authors' contributions

section_content: AJ completed the semblance computation and completed the pre-eruption and eruption analysis. IL and JT completed the inversion of the VLP. BC completed the TOUGH2 analysis. All authors contributed to the conceptual model development and supported the writing. All authors read and approved the final manuscript.
section_name: Authors' contributions

section_content: The corresponding author (Arthur Jolly) is a member of the editorial board and has been blinded from the review process. The other authors of this submission declare that they have no competing interests.
section_name: Competing Interests

section_content: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
section_name: Publishers Note