Long-period seismicity reveals magma pathways above a laterally propagating dyke during the 2014-15 Bárðarbunga rifting event, Iceland

Abstract The 2014–15 Barðarbunga–Holuhraun rifting event comprised the best-monitored dyke intrusion to date and the largest eruption in Iceland in 230 years. A huge variety of seismicity was produced, including over 30,000 volcano-tectonic earthquakes (VTs) associated with the dyke propagation at ∼6 km depth below sea level, and large-magnitude earthquakes accompanying the collapse of Barðarbunga caldera. We here study the long-period seismicity associated with the rifting event. We systematically detect and locate both long-period events (LPs) and tremor during the dyke propagation phase and the first week of the eruption. We identify clusters of highly similar, repetitive LPs, which have a peak frequency of ∼1 Hz and clear P and S phases followed by a long-duration coda. The source mechanisms are remarkably consistent between clusters and also fundamentally different to those of the VTs. We accurately locate LP clusters near each of three ice cauldrons (depressions formed by basal melting) that were observed on the surface of Dyngjujokull glacier above the path of the dyke. Most events are in the vicinity of the northernmost cauldron, at shallower depth than the VTs associated with lateral dyke propagation. At the two northerly cauldrons, periods of shallow seismic tremor following the clusters of LPs are also observed. Given that the LPs occur at ∼4 km depth and in swarms during times of dyke-stalling, we infer that they result from excitation of magmatic fluid-filled cavities and indicate magma ascent. We suggest that the tremor is the climax of the vertical melt movement, arising from either rapid, repeated excitation of the same LP cavities, or sub-glacial eruption processes. This long-period seismicity therefore represents magma pathways between the depth of the dyke-VT earthquakes and the surface. Notably, we do not detect tremor associated with each cauldron, despite melt reaching the base of the overlying ice cap, a concern for hazard monitoring.


Introduction
Long-period seismicity at volcanoes often involves fluid-related processes associated with magma movement (Chouet and Matoza, 2013). Therefore, it is a crucial tool for elucidating the processes occurring in a volcano plumbing system and for eruption forecasting and monitoring (e.g. Chouet 1996;McNutt 2005). Long-period seismicity includes both short-lived events (long-period events or LPs) and longer-duration seismicity, known as tremor, with peak frequencies in the band 0.5-5.0 Hz (Chouet, 1996).  Gudmundsson et al., 2016), eruption fissures shown as orange diamonds; eruption periods shown by orange bands; seismometers indicated by triangles. (For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.) occurring at many places along the dyke path (beneath the ice cauldrons, and at the sub-aerial fissure eruption sites).
We systematically searched for and located both long-period events and tremor during the dyke propagation phase and the first few days of the eruption. We used waveform cross-correlation to detect long-period events, identify similar families of events and to locate tremor. Cross-correlation has an advantage over more traditional onset picking techniques when studying long-period seismicity, which is often emergent in nature (Konstantinou and Schlindwein, 2003). We hope that cross-correlation-based techniques for location of seismic tremor will be incorporated into real-time volcano monitoring.
We investigated how this long-period seismicity is linked to (a) the dyke propagation and associated >30,000 volcano-tectonic (VT) earthquakes (Ágústsdóttir et al., 2016) and (b) sub-glacial eruptions beneath the ice cauldrons formed on Dyngjujökull glacier. To our knowledge, this is the first detailed study of LPs occurring during lateral propagation of a dyke. The excellent seismic network coverage available during this event provided a unique opportunity to investigate in detail whether the LPs and tremor are mechanistically related.

2014-15 Bárðarbunga rifting event
Bárðarbunga is a central volcano on the margin of Iceland's Eastern and Northern Volcanic Zones, underneath Vatnajökull ice cap (Fig. 1A). On 16 August 2014, volcano-tectonic earthquakes started occurring under the ice-filled Bárðarbunga caldera. Over the next 13 days, over 30,000 VT earthquakes occurred at approximately 6 km depth b.s.l., tracking a dyke for 48 km as it propagated laterally north-eastwards. On 29 August, a small, 4-hour eruption occurred in the Holuhraun lava field, reoccupying craters from a 19th century eruption (Hartley and Thordarson, 2013). On 31 August a larger, sustained fissure eruption began, which continued until 27 February 2015, erupting an estimated bulk volume of 1.44 km 3 of lava over an area of 84 km 2 (Pedersen et al., 2017). The eruption on 31 August started with a continuous curtain of fire fountaining along the 1.6 km-long fissure. On 1 September the ac-tivity on the fissure began to localize along specific vents; the main one (shown in Fig. 1) known as Baugur (Pedersen et al., 2017). VT earthquakes continued to occur along sections of the dyke pathway during the eruption. Tremor was detected by the Icelandic Meteorological Office on 23 August at 1120 UTC (during the dyke propagation phase) and for most of the day on 3 September (during the eruption phase) (Icelandic Meteorological Office, 2014). Ice cauldrons were observed on the surface of the glacier during flyovers on 27 August (cauldron DK-02) and 5 September (cauldrons DK-01 and DK-03) .

Data
The microseismicity accompanying the 2014-15 Bárðarbunga rifting event was recorded in unprecedented detail by the University of Cambridge seismic network, which in August 2014 comprised 72 three-component broadband seismometers (see Figs. 1 and S1 in Ágústsdóttir et al., 2016). To complement the network, 14 stations from the national seismic network of the Icelandic Meteorological Office were used, as well as one British Geological Survey station and one University College Dublin station. The network provides good azimuthal coverage, with excellent sampling north of the ice cap. Coverage is less good to the south of the dyke due to difficult field conditions on the ice cap. Primarily, data from stations FLUR (ESP, 30 s-50 Hz response) and DYN (3T, 120 s-50 Hz response) are presented in this study (see Fig. 1A). The data used cover the period from 16 August to 7 September 2014, encompassing the whole dyke intrusion and the first week of the main eruption. All data were recorded at 100 Hz sample rate with a GPS time stamp. This work builds on the earthquake catalogue from Ágústsdóttir et al. (2016), comprising over 30,000 earthquakes.

Methods
In the following section, we first describe the steps taken to investigate long-period events (LPs): an initial search for lowfrequency earthquakes in an existing earthquake catalogue; a second search using a template waveform cross-correlation method; identification of repetitive families of events using waveform crosscorrelation; accurate location of LPs using manually-picked phase arrival times and double-difference relocation; and finally inversion for the focal mechanisms of the LPs. We then describe the method used to locate tremor.

Earthquake frequency content analysis using Frequency Index (FI)
An automatic procedure was used to analyse the frequency content of the earthquakes from Ágústsdóttir et al. (2016). A Frequency Index (FI) was calculated for each event using the method from Buurman and West (2010). The FI is a ratio of an earthquake's high-and low-frequency content, with subjectively defined frequency bands. A logarithm is taken such that earthquakes with equal amounts of energy in the high-and low-frequency bands are assigned an FI of 0, then negative FI values represent earthquakes dominated by energy in the lower-frequency band, and positive FI values the higher-frequency band.
where A is the average amplitude for the frequency band. A correction was applied to compensate for attenuation, which is both distance and frequency dependent. This was approximated by a simple static correction, estimated from FI values versus hypocentral distance (see Supplementary Section S1, Fig. S1, for details). The FI is a more robust analysis than dominant frequency, which does not capture complexities (such as bi-modal distributions) in the spectra and is often biased by noise.
FI values were calculated from the vertical component of six high-quality local stations, over a window from 2 s before the Pwave pick time to 23 s afterwards. The low-and high-frequency bands were visually selected (2-4 Hz and 5-16 Hz respectively) to enable differentiation between manually-identified, representative earthquakes. The low-frequency band was cut off below 2 Hz to avoid any contamination by the strong microseism signal. FI values from station FLUR are presented here (Fig. 1). Absolute FI values varied across stations, in some cases quite significantly, but the observed trends were consistent. This is expected, since complex source geometry and differing station path and site effects typically result in different frequency content recorded for the same earthquake at different stations.
Horizontal components were not used here, being more contaminated by noise, but showed the same trends. While P-wave energy dominates the vertical component, higher amplitude, lower frequency S-wave energy dominates the horizontals, so the horizontal components exhibit lower FI values. Analysis using the vertical component has the advantage of including the clean, higher frequency, brittle onset and also later, lower frequency energy. These features are important for the template waveform method described in Section 3.2, below.

Long-period event detection with template waveform cross-correlation
To obtain a comprehensive catalogue of LPs occurring during the study period (i.e. detecting events that could have been missed in the automatically-generated dyke catalogue from Ágústsdóttir et al., 2016), a template waveform cross-correlation method based on Shapiro et al. (2017) was used. Events were searched for on the vertical component during the period from 16 August to 7 September (inclusive). The absolute value of the trace was calculated then smoothed over a 1 s window. Initial triggers were found in the smoothed trace when the short-term average (STA) in a 1 s window was greater than three times the long-term average (LTA) in the preceding 50 s window. The maximum point in the trace was found between 15 s before and 21 s after this trigger, because the events correlate better when aligned based on their maximum amplitude in the next cross-correlation step. Then the event was cut at 15 s before and 21 s after the maximum. To exclude VTs, the FI was calculated for the event. A cut-off FI of −0.75 was used for station FLUR, visually selected so as to encompass the cloud of lower frequency (low FI) earthquakes shown in Fig. 1 (with station FLUR marked).
Next, master (or template) events were identified by constructing a cross-correlation matrix of all event waveforms. For every waveform, the sum of its correlation coefficients (CCs) with all of the other events (each row of the CC matrix) was computed. The signal with the largest sum of its CCs (CCsum) was then selected as the master event and grouped into a family with all of the other waveforms that had a CC of above 0.33 with the master. A relatively low CC value was used here in order to find as many LPs as possible. Shapiro et al. (2017) also use a relatively low CC of 0.44 at this stage. This procedure was repeated iteratively for the remaining events in the matrix, to identify all families. The master events from the families with more than ten members were used as template events in the next step.
Finally, to produce the catalogue of LPs, the template waveforms were slid along the trace for the entire period and events detected when the correlation coefficient exceeded 0.33. A more accurate onset time was then searched for, so that later clustering of highly similar events -i.e. using a much higher correlation coefficient of 0.7 -was more successful (described in Section 3.3). An STA over 0.8 s and an LTA over 8 s was used, picking a time when the STA/LTA exceeded 5. If this failed, we searched over a window starting 10 s earlier, picking when the STA/LTA exceeded 3. If this too failed the events were discarded. For a final time, the FI of the waveform was measured and only long-period events (FI of <−0.75) were accepted.
This procedure was carried out for six local stations independently, although the resultant LP catalogue comprised events identified on the closest and clearest station, FLUR. We note that the method may not be suitable for studies where very emergent LPs are present because of the use of an onset picker in the initial detection stages; cross-correlation with template events will only find emergent LPs if the template events are themselves emergent. We would instead suggest using an initial detection method based on peaks in the frequency spectrogram.

Identifying clusters of highly similar, repetitive LPs
Clusters of highly similar, repetitive LPs were identified by waveform cross-correlation, using the GISMO toolkit (Reyes and West, 2011). Traces were band-pass filtered between 2 and 16 Hz and cross-correlated over a window 0.5 s before to 6 s after the P wave pick time. Clusters were then defined by a correlation coefficient of 0.7 or greater. These parameters were selected by visual inspection, and are somewhat arbitrary; for the same correlation coefficient a short window around the onset results in larger clusters with highly similar P arrivals but variation in the coda, whereas a long window finds smaller clusters of events highly similar throughout the waveform. A window length in-between these two end members was chosen such that events were highly similar over both the P and S arrival but the coda was allowed to vary. Onset times of the events were corrected by aligning the waveforms based on their maximum correlation in each cluster.

LP manual location
A subset of > 100 randomly-selected LPs were relocated by manually picking phase arrivals. The non-linear location software NonLinLoc (Lomax et al., 2000) was used, with a 1-D velocity model from Ágústsdóttir et al. (2016). Not all events had clear S phase arrivals, and so to check the depth accuracy of location using P phase arrivals alone, a sub-sample of clear, VT earthquakes were located both with and without their S phase arrival picks. The resultant hypocentres were comparable within uncertainty (as described in Supplementary Section S2.2), giving confidence in locations derived using P-phases only. The earthquake locations were further refined by relative relocation with HypoDD (Waldhauser and Ellsworth, 2000). Locations are given in Supplementary Table  S1. Absolute uncertainties are typically 1.1 km in latitude and longitude, and 2.8 km in depth, shown in Supplementary Fig. S2 (with synthetic tests shown in Fig. S3).

LP focal mechanisms
Moment tensor solutions were obtained using MTfit (Pugh et al., 2016), a Bayesian source inversion method. Phase arrival polarities were manually determined at each station across the network and the pick confidence translated into a polarity probability density function (PDF). A double-couple (DC) solution could be imposed, or the fitting could be left unconstrained (non-DC).
Unfortunately, the network configuration results in poor coverage towards the centre of the focal sphere, so it is not clear whether a non-DC component is required. In an attempt to constrain one way or the other, amplitude ratios and SH arrival polarities were explored. S wave arrivals were typically emergent and therefore the polarity picks, and the SH polarity inversions, were of poor quality and unreliable. Including amplitude ratios in the inversions did not appear to influence the solutions or to improve the constraint. Therefore, the presence of a volumetric component could not be ruled out, but where a DC solution could be fitted it was assumed to be the primary failure mechanism.

Tremor detection and location
The freely available program MSNoise (Lecocq et al., 2014) was used to compute cross-correlation functions (CCFs) between pairs of seismic stations. Up to 20 stations were used (see Supplementary Section S3 for details). Only vertical components were used, and were preprocessed individually for each station. An initial bandpass filter was applied (0.01 to 8.0 Hz), and the waveforms were demeaned, tapered and downsampled to 20 Hz. The waveforms were then temporally normalized using 1-bit normalization and spectrally whitened in 30-min windows. After another bandpass filter was applied (0.5 to 2.0 Hz), CCFs were calculated for every pair in the network in 30-min windows for time lags of ±120 s. These CCFs were stacked using phase weighting over the relevant time periods (method from MSNoise, 2017).
Following the method of Ballmer et al. (2013) and Donaldson et al. (2017), a 41 by 41 grid was constructed over the geographical area (approx. 48 by 56 km) and each grid point was considered to be a potential source location. Arrival times were estimated in the CCFs for each pair of stations for each grid point, assuming a lateral seismic propagation velocity of 1.2 km/s. The total absolute amplitude in a 2 s window around the expected arrival time (using the differential interstation distance) was calculated for each pair's CCF and summed over all pairs. This then represents the likelihood of the source being located at that grid point. Summing in this way means that it is possible to locate more than one seismic source at the same time. A range of velocities between 0.6 and 3.6 km/s were tested and it was found that ∼1.2 km/s generally produced the highest amplitudes in the CCFs summed across all pairs (more detail is given in Supplementary Section S3.3). This value is also consistent with other studies of Icelandic volcanoes at similar frequencies (Li et al., 2017, and references therein).
Uncertainty in location is estimated from the 0.75 contour of the normalized interpolation of location likelihood. This uncertainty was tested during periods when the location of seismicity is known: the LPs in cluster 1 and the tremor at the eruption fissure (described in Section 4). For 92% of events in cluster 1, the manually-picked location fell within the 0.75 contour of the location estimate using 10-minute CCFs ( Supplementary Fig. S4). Note that when locating LPs using this technique, the waveforms were normalized before cross-correlating by clipping at the root mean square value multiplied by 3. For the first 10 days of the eruption, the main eruption vent fell within the 0.75 contour of the day-stack of CCFs. Furthermore, jack-knife tests were performed to test the network bias in the location estimates shown in Fig. 3 (Wech and Creager, 2008). The most likely source location was calculated 1000 times for randomly chosen network configurations using three-quarters of the number of stations, for each time period. Over 95% of results located within the 0.75 contours.
To search for tremor, we visually observed the location plots for CCFs stacked over each day as well as over two-hour windows throughout the study period (see Supplementary Movie 1). For periods when the uncertainty in location was reasonably small, the waveforms and spectrograms were checked to identify whether tremor occurred. It was difficult to detect weak periods of tremor without detecting other types of seismicity; because the VTs have some frequency content below 2.0 Hz and were so numerous, these could be picked up in the CCFs, as could large earthquakes at Bárðarbunga caldera. This effect was reduced when using 1-bit normalization, removing all amplitude information from the trace, but then weaker periods of tremor were missed. We have compromised by studying various stack lengths of CCFs and testing both 1-bit and 3 x RMS normalization. Our methodology for detecting and locating tremor is therefore not a fully automatic approach; this will be the focus of future work. Existing studies have suggested automatic criteria for tremor detection based on a normalized value of the likelihood of source location (Droznin et al., 2015) and on network coherency estimated from the network covariance matrix (Seydoux et al., 2016;Soubestre et al., 2018).

Results
The frequency content of earthquakes along the dyke path (during intrusion and the onset of eruption) was analysed and is shown in Fig. 1. A locus of lower frequency events is clearly visible under the glacier, shallower than the majority of the dyke seismicity.
These low frequency events occurred over ∼10 days from 25 August, after propagation phases of the dyke tip (i.e. during stalled phases) at a step between dyke segments, approximately 1.5 km north of the northernmost ice cauldron (DK-03) that formed during the rifting event. A smaller number of low frequency events were also observed further south, occurring in short bursts on the 19 and 20 August. In Section 4.1, we show that many of these low frequency earthquakes make up clusters of repetitive, highly similar events, locations of which are shown in Fig. 2.
Tremor (continuous low-frequency seismicity) was also detected during this period, in three locations shown in Fig. 3. On 23 August, tremor was observed briefly for 1 hour in the vicinity of the ice cauldrons, with the most likely location closest to the central cauldron (DK-02). The tremor had low-frequency content (∼0.5-2.0 Hz) but the spectrum (Fig. 3) at this time is dominated by the VT earthquakes associated with a propagation phase of the dyke (see Fig. 2). Following this, tremor was not detected again until the onset of the main eruption in Holuhraun, when it was located at the eruption vents (represented by 4 September in Fig. 3).
The peak frequency was ∼1.0 Hz with an overtone at ∼3.0 Hz at FLUR. On 3 September, the vent tremor was masked for 18 hours by a stronger tremor signal in the vicinity of the northernmost cauldron (DK-03) and the low frequency earthquakes. Again, the peak frequency was ∼1.0 Hz but now with several overtones visible. The tremor source is likely to be shallow (less than ∼1-2 km below surface). This is because the largest amplitude arrivals in the cross-correlation functions correspond to velocities of ∼1.2 km/s ( Supplementary Fig. S5), consistent with surface waves dominating the tremor.

Clusters of LP earthquakes
Over 450 long-period events were detected, 200 of which make up eight main clusters. Therefore, during the dyke intrusion phase, ∼1% of the earthquakes were long-period. The eight largest clusters are presented here, with cluster 1 containing the largest number of events, cluster 2 the second largest, and so on (crosscorrelation matrices shown in Supplementary Fig. S6).
Cluster 1, comprising 63 events, occurred over ∼10 days from 25 August, coincident with a number of smaller neighbouring clusters. We find that the clusters were offset by 1 km to the west and were (∼2 km) shallower than the majority of the dyke seismicity, 1.5 km north of the northernmost cauldron (DK-03), at the southern end of the main graben (Fig. 2). The earthquake waveforms have a high frequency, brittle onset, an often emergent S arrival, and a low frequency, long duration (>20 s) coda (Figs. 4, 5 and   6). The dominant frequency is ∼1 Hz (Fig. 7), so hereafter they will be referred to as long period events (LPs). Two of the clusters identified (clusters 4 and 6) -those which occurred in short bursts on 19 and 20 August respectively -were located further south along the dyke path, where the network coverage is sparser. These LPs have similar waveform characteristics, although a higher dominant frequency of 3-4 Hz. As with the northern clusters, both were slightly laterally offset from the dyke seismicity. In this case, the LPs appear to have been at similar depths to the dyke seismicity (between 4 and 8 km depth b.s.l.), however synthetic testing suggests that they in fact occurred shallower than this (see Supplementary Fig. S3 B).

LP earthquake failure mechanisms
Focal mechanisms were generated for the manually-picked LPs and are shown in Fig. 8. Double-couple (DC) constrained inversion found that all the earthquakes in all the clusters of LPs could be fitted by reverse faulting, with a strike approximately perpendicular to the dyke. Unconstrained moment tensor inversion suggests a closing crack, with an implosive volumetric component. Since all the earthquakes could be fitted by a DC solution, we consider reverse faulting to be the primary source mechanism. See Supplementary Fig. S14 for the un-constrained solution. For cluster 1 (and neighbouring northern clusters with well-constrained mechanisms -2, 3, 5, 8) the fault plane is assumed to be that which dips toward the northern cauldron, DK-03 (i.e. striking between 180 and 360, coloured in Fig. 8). This WNW-ESE striking reverse faulting is a distinctly different mechanism to that of the dyke VTs, which arise exclusively from strike-slip failure sub-parallel to the dyke (Ágústsdóttir et al., 2016).

Discussion
We have found that at least 450 long-period events, classified by Frequency Index (see Methods), occurred during the Bárðarbunga-Holuhraun intrusion and the first seven days of eruption. Approximately 200 of these events make up eight main clusters of highly similar, repetitive events. For most clusters, peak frequencies were ∼1 Hz and source mechanisms were WNW-ESE striking reverse faulting. The clusters occurred beneath Dyngjujökull glacier -where three ice cauldrons (DK-01, DK-02, DK-03) formed -and were located 1-2 km shallower than the volcanotectonic (VT) earthquakes associated with the dyke propagation.
We have also detected and located long-period seismic tremor during three time-periods: 1 hour on 23 August near the central ice cauldron (DK-02); 18 hours on 3 September near the northernmost cauldron (DK-03); and continuously during the main eruption (31 Aug onwards) near the main vent. Further discussion of the tremor locations is made in Section 5.3.

Fluid process
We suggest that the long-period events were caused by processes related to magmatic fluids (be it the magma itself or exsolved gases). The main evidence for this is that the LPs mainly occurred during phases of dyke-stalling (see Fig. 9), when the dyke was presumably pressurizing before breaking through the crust for the next propagation phase. We infer that during these times of pressurization, magma forced its way upwards above the dyke, generating the LPs. The lack of LPs during propagation phases is not an artefact arising from dominance of VT events, as magnitudes of the LPs are comparable to those of the VTs (and above the magnitude of completion).
A possible alternative explanation is that the resonance of the LPs arises from a hydrothermal fluid, rather than a magmatic one (Petersen, 2007;Matoza et al., 2014). The presence of the overlying glacier and several metres of rifting associated with the dyke intrusion ( Ágústsdóttir et al. (2016) in grey; subset of manually-located LPs in dark grey or coloured by cluster (green -principal cluster, 63 events), with cluster centroid location given by coloured square (determined from average location of manually-refined events in each cluster). Location uncertainties given in Supplementary Section S2.2, Fig. S2. Tremor contours for 23 Aug and 3 Sept (from Fig. 3) shown by brown lines, dashed and solid respectively, with maximum values given by brown stars. Glacier in light blue; graben in dark blue (digitised from Rossi et al., 2016) and surface fractures in navy (Hjartardóttir et al., 2016); ice cauldrons indicated by dark purple diamonds, eruption fissures by orange diamonds; eruption periods shown in orange. Fissure and cauldron locations (not timings) indicated by dashed lines in D). LP locations given in Supplementary Table S1. to percolation of meltwater down through the crust. This seems an unlikely source mechanism though, given that the LPs were located at 4-5 km depth b.s.l. Nevertheless, it may well be that the LPs and tremor were in some way linked to the rifting and opening of the dyke, since the final dyke opening model from Sigmundsson et al. (2015) shows that the clusters of LPs occurred beneath a region of maximum opening (see Supplementary Fig. S15). LPs observed during dyke intrusion events at Krafla volcano in Iceland seem to be correlated with the opening of eruptive or non-eruptive surface fissures (Brandsdóttir and Einarsson, 1992).
It has been suggested that long-period seismicity can be generated by processes not requiring the presence of fluids, for example tectonic tremor (Obara, 2002), although this is less relevant here. In volcanic environments, LPs may in fact represent brittle events which have undergone strong path effects (Harrington and Brodsky, 2007;Bean et al., 2008;Coté et al., 2010). In particular, a long duration, low frequency coda could result from trapped waves in loosely-consolidated, shallow layers of the crust (Bean et al., 2008). We observed an increase in length of LP coda with increasing station distance (see comparison of stations FLUR and DYN in Fig. 5). This has been seen elsewhere in Iceland (Sgattoni et al., 2016), and shows the effects of scattering and path length. However, in our case path effects cannot fully explain the long-duration, lowfrequency coda of the LPs, because the network also records VT earthquakes close by, with very similar paths and without the long coda. Another alternative mechanism is that slow-rupture of the non-consolidated material itself could produce LPs (Bean et al., 2014). Again, the depths of the LPs observed in this study are too great for this to be a viable mechanism. Aki et al. (1977) suggest that jerky extensions of a fluid-filled crack, driven by the excess pressure of fluid, could produce longperiod seismic tremor. Later studies suggest that LPs could be  Waveforms shown for all 3 components on stations FLUR (left) and DYN (right) with P and S arrivals marked in red and blue respectively and followed by a low frequency, long duration energy packet. Typical dyke VT earthquake example given for comparison, with location marked in Supplementary  Fig. S2. produced by the same mechanism (e.g. Fehler, 1983). A crucial observation in this study is that the LPs were highly similar and occurred over periods of up to 10 days, which implies that they were generated by repetitive excitation in the same location (see also Neuberg et al., 1998Neuberg et al., , 2006Petersen, 2007;Richardson and Waite, 2013;Hurst et al., 2014;Shapiro et al., 2017). This rules out the possibility of progressive fracturing of a crack over distances greater than approximately a quarter of the dominant wavelength (∼250 m, Neuberg et al., 1998).

Fluid-filled cavity
Instead, a more likely scenario is that the LPs were produced by resonance of a fluid-filled cavity (Chouet and Matoza, 2013). Such cavities have been suggested to be either planar (dyke or crack) shaped (Chouet, 1996) or conduit (tube) shaped (Neuberg et al., 2000(Neuberg et al., , 2006. Which of these two geometrical models is preferred depends on the compressibility of the fluid; a crack requires gasrich basaltic magma, misty water-steam mixtures, or dusty ash-gas mixtures (Kumagai and Chouet, 2001), whereas a conduit requires bubbly magma (Neuberg et al., 2006).
The trigger (or kick) that initially excites a fluid-filled cavity is the subject of various models. Excitation may arise from transient pressure changes in the fluid, caused for example by fluid heterogeneity or channel geometry (Julian, 1994). Alternatively, the trigger could be caused by brittle failure of magma at the conduit walls (Neuberg et al., 2006). Some heterogeneity in the conduit geometry, such as a bottleneck or kink, is required for this shear-failure to occur at relevant depths in the crust (Thomas and Neuberg, 2012). In this scenario, the event rate is therefore related to the magma ascent rate.
A key observation in this study is that the LPs exhibit P and S arrivals, followed by a higher-amplitude, low-frequency energy packet (see Figs. 4 and 5). We suggest two potential origins of this phase: a surface wave or cavity resonance. In the models of Neuberg et al. (2000,2006), the LP trigger sets up an interface (or tube) wave at the fluid-rock boundary, which travels up the conduit boundary, before 'leaking' at the free surface and radiating energy outwards as surface waves (Neuberg et al., 2000). The particle motion (Supplementary Movie 2) and timing ( Supplementary  Fig. S16) of this later wave packet are consistent with that of a retrograde, elliptical surface wave. Further, the spectrograms in Fig. 5 suggest that the energy may be dispersive. It is important to note that, because any surface waves presumably travel along the rockice boundary, usual assumptions of particle motion and dispersion, which assume that seismic velocity increases with depth, may not hold (since ice introduces a high velocity layer overlying lower velocity upper crust). A constraint in this model is that the top of the tube needs to be within approximately one wavelength (∼1-3 km) of the surface to generate surface waves. The tube would then need to have a vertical extent of a few kilometres for up to 10 days while the LPs are generated.
An alternative is that the later phase in the LPs represents resonance of a fluid-filled cavity. We point out that the LPs do not display simple decaying harmonic oscillations as might be expected in this model (e.g. Fig. 1 of Kumagai and Chouet, 1999). Waveform modelling could address these questions but is not at-tempted in this study owing to a) the complexity of the setting (e.g. a high-velocity ice layer above a lower-velocity rock layer) and b) the non-uniqueness of the model, i.e. unless the compressibility of the fluid or the geometry of the cavity is known then only end members can be constrained. However, the large amplitude and late arrival of the low frequency energy packet, along with its retrograde, elliptical particle motion, is strong evidence of surface waves and therefore we favour the 'tube' model.
Waveform characteristics, spectral content and focal mechanisms are remarkably similar across all clusters of LPs, suggesting some fundamental source process which is distinctly different to the sub-parallel strike-slip faulting of the dyke VTs. The moment tensor inversions show either a reverse shear-fault striking subperpendicular to the dyke, or a closing crack (implosive volumetric component). Both of these possibilities would be surprising if the local stresses associated with the rifting event were controlling the LPs. Instead, we suggest that this is further evidence of the influence of magmatic fluids on the source mechanism. The consistently vertical T axis (minimum instantaneous strain direction, though the minimum stress direction must also be vertical or subvertical) could be due to buoyancy forces due to the presence of magma. It has been suggested that moment tensor inversion (from phase arrival polarities) alone is not appropriate when mass advection processes are involved (Chouet and Matoza, 2013), and that consideration of single forces -i.e. the momentum transfer between source and earth -is required in addition to moment tensor components. In our case, deconvolution of the source (by full waveform inversion) is not feasible without better knowledge of the seismic velocity model. It is possible, though, that the reverse faulting solution arises largely from the (single) upward force of the magma, which could help explain the focal mechanism consistency across spatially-distinct LP clusters. Many studies in Iceland report reverse faulting in the presence of magma, e.g. in Upptyppingar (White et al., 2011), Eyjafjallajökull (Tarasewicz et al., 2012), Vaðalda  and Bárðarbunga (Hudson et al., 2017). White et al. (2011) suggest that reverse and/or normal faulting may arise from fracture within frozen melt. In this scenario, a low melt supply allows freezing of the magma, which increases in strength as it cools through the glass transition , in turn leading to constriction of the flow, high stressing rate, and fracture of the frozen magma. The fault plane solutions would therefore reflect the geometry of the LP trigger mechanism, such as a constriction in the conduit or crack, and the later, low-frequency energy packet is produced by an interface (tube) wave or cavity resonance.

Source models of tremor
We focus on the tremor on 23 August and 3 September. Models of volcanic tremor can broadly be divided into five categories: excitation of fluid-filled cavities, fluid-flow induced oscillations, hydrothermal boiling, repetitive fracturing of unconsolidated material and resonance of large magma bodies (e.g. Konstantinou and Schlindwein, 2003;Eibl et al., 2017b). The tremor source is likely to be shallow, given that the propagation velocity between pairs of stations as measured in the cross-correlation functions is 1.2 km/s (see Section 4). We therefore consider the first four options to be viable mechanisms. Before considering in more detail the possible tremor source mechanisms, we first consider three main pieces of evidence that the LPs and tremor observed in this study are related. We focus on the six northerly clusters of LPs and tremor on 3 September in this section, but the arguments also hold for cluster 6 and the tremor on 23 August. Firstly, within uncertainty, the LPs and tremor occurred in the same lateral location. Secondly, there appears to be a temporal link between the two. The LPs began on 25 August (see Fig. 2) and occurred in swarms (see histograms in Fig. 9). The final swarm, beginning on 2 September, ended with the onset of the 3 September tremor. Only 7 more LPs occurred from these clusters after the tremor stopped. Finally, the spectral content of the LPs and tremor are very similar (Fig. 7), with peak frequencies of ∼1 Hz. This varies slightly between stations, showing either the effect of different paths and sites, or that the source radiation pattern is not uniform.
Therefore, it is possible that the tremor comprises rapid, repeated excitation of LPs, as has been observed in other studies (e.g. Fehler, 1983;Neuberg et al., 1998). The tremor would then represent sustained excitement of the fluid-filled cavities detailed in Section 5.1.2. In this case, picking the P-and S-wave arrivals of the LPs reveals the location of the trigger which excites the cavity, whereas the tremor location represents the energy which escapes at the shallow end of the cavity. The tremor could have occurred during the sub-glacial eruptions associated with the ice cauldrons, when a higher fluid flow rate resulted in continuous excitation. A model of a cavity producing both the LPs and tremor requires a cavity which extends from the depth where the LPs are triggered (4-5 km b.s.l.) to where the tremor occurs (inferred to be less than 1-2 km below surface, see Section 5.1.2) -i.e. with a vertical extent of several kilometres.
Alternatively, it could be coincidental that the spectral contents of the LPs and tremor are similar and, in fact, different mechanisms produce the two types of seismicity. In this scenario, it may still be that the location and timing of the LPs and tremor are related in that the LPs represent melt moving above the dyke towards the surface, which culminates in tremor-producing sub-glacial eruptions. Tremor-producing processes related to the sub-glacial eruptions could include fluid-flow induced oscillations (Julian, 1994), hydrothermal boiling (Leet, 1988) or repetitive fracturing of unconsolidated material (Eibl et al., 2017b). It has been shown that magma-induced tremor from fluid-flow instability requires very narrow cracks and high flow velocities, so hydrothermal or gasrich fluids may be more likely (Hellweg, 2000;Balmforth et al., 2005). Hydrothermal boiling is an attractive option given that subglacial eruptions are known to have occurred. Tremor is generated by bubble collapse (Konstantinou and Schlindwein, 2003), which is expected to produce white noise. However, assuming that the water is contained in a channel, the boiling process may set it into resonance, so that the frequency response of the channel Fig. 8. DC fault plane solutions (FPSs) for LP clusters 1 to 8. Left) average fault plane solution; right) rose diagram. Nodal planes 0-180 in grey and nodal planes 180-360 coloured by cluster, with average strike/dip respectively given above. Average FPSs not calculated for clusters 6 and 7 as solutions inconsistent (and not as well constrained). Strikes, dips and rakes given in Supplementary Table S3. will be convolved with the broadband noise (Konstantinou and Schlindwein, 2003). Finally, Eibl et al. (2017b) suggest that the tremor on 3 September is produced by progressive fracturing of rock in the top ∼2 km of the crust. They model the tremor as repetitive 1 Hz Ricker wavelets, but do not address why brittle failure of dry rock would produce a 1 Hz-earthquake, and Ágústsdóttir et al. (2016) do not locate any brittle earthquakes at these depths.
In Iceland, sub-glacial tremor was observed during a fissure eruption beneath Vatnajökull icecap in 1996 (Einarsson et al., 1997;Konstantinou, 2002) and at Katla volcano in 2011 . Konstantinou (2002) suggests that the tremor is caused by ascending turbulent slugs of magma separated by small intervals of laminar flow, but notes that this is uncertain. Sgattoni et al. (2017) suggest that tremor located near ice cauldrons was caused by volcanic processes, possibly hydrothermal. These studies show that, while volcanic tremor may be a useful tool for eruption monitoring, it is still difficult to constrain exact source mechanisms, even with excellent seismic network coverage.

Tremor location: comparison with existing work
Tremor was detected on 23 August and 3 September. The most likely location on 23 Aug, shown in Fig. 3 and Supplementary Movie 1, is nearest to DK-02, but all cauldrons fall within uncertainty. The location uncertainty is large, likely due to contamination from a large number of VT earthquakes associated with a rapid dyke propagation phase. Cauldron DK-02 must have formed before 27 Aug , when it was first observed, and so we propose that the 23 Aug tremor is associated with DK-02. The tremor location on 3 Sept is closest to DK-03 (although DK-02 is also within location uncertainty). Caudron et al. (2018) use Seismic Amplitude Ratio Analysis (SARA) to locate continuous seismic energy during the dyke intrusion and find that while the SARA location usually tracks the dyke   10. Interpretive cartoon. Long-period earthquakes (LPs) comprise a brittle trigger (P and S arrivals), indicated by coloured stars, followed by a long-duration coda, thought to be resonance or tube waves (see text), shown as white dashes. Tremor associated with the formation of an ice cauldron and the main eruption are labelled. Data is also shown: refined LP cluster centroid locations given by coloured squares, automatic dyke catalogue in grey (labelled as dyke VT (volcano-tectonic) earthquakes). Regions of magma flow shown in orange, and regions of perhaps frozen (or non-flowing) magma in brown. Glacier in blue, bedrock in dark grey.
VT earthquakes, on 23 August and 3 September lower frequency sources (tremor) are detected in the vicinity of the Dyngjujökull ice cauldrons. Eibl et al. (2017b) locate the 3 Sept tremor with array beamforming techniques and report it as occurring beneath the central cauldron (DK-02). However, the location uncertainty spans a wide angle and corrections for array squinting are not applied (as later detailed in Eibl et al., 2017a), so the beam-forming location may in fact align with the northern cauldron DK-03.
We are impressed with the potential of our location method for near-real-time tremor monitoring, given the simplicity of the technique, but recognise that more work needs to be done to definitively pinpoint the tremor to a specific cauldron. It is perhaps surprising that the highest amplitude tremor (3 Sept) seems to be associated with the smallest sub-glacial eruption (cauldron DK-03, by volume, see Reynolds et al., 2017) and that we do not find any tremor associated with the largest sub-glacial eruption (cauldron DK-01). Even given our tremor location uncertainty, we do not detect a period of tremor corresponding to each cauldron. This highlights the fact that movements of magma may not be detected by current methods available to volcano seismology.

Conclusion
We have detected and located long-period events (LPs) and tremor associated with the Bárðarbunga-Holuhraun 2014-15 rifting event. We argue that the clusters of highly-repetitive LPs result from triggers next to a fluid-filled cavity, resulting in waveforms with clear P and S phases followed by a long duration, low frequency coda. The clusters of LPs have been accurately located near to and below three ice cauldrons on Dyngjujökull glacier which were caused by sub-glacial eruptions , at shallower depths than the volcano-tectonic earthquakes associated with lateral dyke propagation. Given this, and that the LPs occurred in swarms during times of dyke-stalling and hence pressurization, we suggest that they represent magma moving upwards sub-vertically, rather than horizontally along the main dyke. The waveforms, spectra and fault plane solutions are remarkably consistent between clusters, which is evidence of the same source mechanism occurring over 10 km apart. We have also located two periods of tremor in the same region as the two northerly ice cauldrons; this tremor and the associated sub-glacial eruptions represent the climax of the vertical melt movement below the glacier.
We find that the timing, location and spectral content of the LPs and tremor are similar, and suggest that the tremor may arise from either sub-glacial eruption processes, or comprise rapid triggering of LPs. A cartoon representing our interpretation is shown in Fig. 10. We do not identify a period of tremor corresponding to the southernmost (and largest) cauldron, which highlights the fact that many movements of magma remain 'silent' to current methods available to volcano seismology. lowship and an F.R.S.-FNRS Chargé de Recherches/Université Libre de Bruxelles fellowship. We thank Sveinbjörn Steinþórsson and others who assisted with fieldwork in Iceland, and Jürgen Neuberg, Andy Hooper, Jean Soubestre and others for their helpful discussions. The Icelandic Meteorological Office, Chris Bean (University College Dublin), and the British Geological Survey kindly provided additional data from seismometers in northeast Iceland, data delivery from IMO seismic database 20151001/01. Earthquake hypocentres can be found in Ágústsdóttir et al. (2016). The raw seismic data are archived at Cambridge University and will be available at IRIS for download from October 2019. The GISMO toolbox is freely available from https://geoscience -community-codes . github .io /GISMO/ and the MTfit software for moment tensor inversion is freely available from https://github .com /djpugh /MTfit. We thank reviewers Nikolai Shapiro and Tony Hurst for their constructive comments. The first two authors contributed equally to this paper. Department of Earth Sciences, Cambridge contribution ESC.4091.

Appendix A. Supplementary material
Supplementary material related to this article can be found online at https://doi .org /10 .1016 /j .epsl .2018 .03 .020.