Closing crack earthquakes within the Krafla caldera, North Iceland

Moment tensor analysis with a Bayesian approach was used to analyse a non-double-couple (non-DC) earthquake ( M w ∼ 1) with a high isotropic (implosive) component within the Kraﬂa caldera, Iceland. We deduce that the earthquake was generated by a closing crack at depth. The event is well located, with high signal-to-noise ratio and shows dilatational P -wave ﬁrst arrivals at all stations where the ﬁrst arrival can be picked with conﬁdence. Coverage of the focal sphere is comprehensive and the source mechanism stable across the full range of uncertainties. The non-DC event lies within a cluster of microseismic activity including many DC events. Hence, we conclude that it is a true non-DC closing crack earthquake as a result of geothermal utilization and observed magma chamber deﬂation in the region at present.

The Krafla central volcano lies within the Northern Volcanic Zone (NVZ) of Iceland (Fig. 1A).The divergent plate boundary is expressed as a series of en echelon volcanic systems, each comprising a fissure swarm and central volcano.The Krafla central volcano has developed a caldera and is transected by an NNE-SSW trending fissure swarm.The last major rifting episode in Krafla occurred from 1975 to 1984, with a series of eruptions and dyke intrusions at Krafla (Björnsson et al. 1985;Einarsson 1991;Buck et al. 2006).
A shallow magma chamber at 3 km depth has been detected in Krafla caldera by the attenuation of shear waves (Einarsson 1978) and compressional wave delays (Brandsdóttir et al. 1997).The magma chamber is approximately 2-3 km in extent northsouth, 8-10 km east-west and is 0.7-1.8km thick (Brandsdóttir et al. 1997).Surface deformation modelling, using geodetic, tilt and levelling data, shows that the magma chamber has been deflating since 1989 (Sturkell et al. 2008 andreferences therein).In 2009, the IDDP (Iceland Deep Drilling Project) drilled into a pocket of rhyolitic melt 2.1 km below the surface (1.6 km below sea level) (Elders et al. 2011).Within the Krafla caldera, there are high levels of microseismicity, mainly associated with magma cooling and geothermal activity (Heimisson et al. 2015;Schuler et al. 2015) within the steam-dominated Krafla geothermal system (Gudmundsson & Arnórsson 2002).Foulger et al. (1989) observed microearthquakes in the Krafla area that had a variety of non-DC mechanisms, including explosive tensile-crack events and implosive events.They hypothesized that this range of source mechanisms indicated a heterogeneous stress regime, following a crustal spreading episode from 1975 to 1984.The physical mechanisms for the implosive events were hypothesized to be due to cavity collapse at depth.Ford et al. (2009) observed examples of non-DC events (M w 3.4-4.0)with a large implosive component (closing crack) in Utah.These events all occurred following mine collapses; hence the source mechanism was known prior to moment tensor analysis.
Several well-constrained, strongly non-DC events were recorded by a dense, temporary seismic array at Krafla during routine analysis of detected microearthquakes (Fig. 1A).We focus on one  (Einarsson 1978) in purple, faults and fissures in brown and eruption sites in red (adapted from Hjartardóttir et al. 2015).Blue/white symbols show locations of seismometers: blue triangles are temporary stations operating from 2009 July, white triangles temporary seismometers operational from 2010 July to 2011 September and blue stars Iceland Meteorological Office permanent stations.Black box shows extent of (B).Inset shows location on a tectonic map of Iceland (Einarsson & Saemundsson 1987) with fissure swarms (yellow) and volcanic centres.(B) Map and orthogonal cross-sections of earthquakes discussed in the text.Events are colour-coded according to source mechanism.Red triangles show the location of the temporary seismometers.Event A is the closing crack event and Event B a reference double-couple event.Colour contours on the cross-section are modelled V p /V s ratios from Schuler et al. (2015), rotated by 8 • to align with orthogonal cross-sections.The low V p /V s ratio region is indicative of the presence of superheated steam above rhyolite melt.event in particular for which we observe only dilatational P-wave first arrivals across the focal sphere, indicating a negative isotropic component to the source mechanism.This microearthquake occurred within a cluster that includes events that show stable, purely DC, moment tensor solutions (Fig. 1B).An example of a wellconstrained DC event, spatially close to Event A, is analysed for comparison.This allows us to be confident that we have detected a true non-DC event, rather than obtaining spurious non-DC solutions resulting from uncertainties in the velocity model or hypocentre locations (suggested by Frohlich (1994) as a cause of non-DC mechanisms).

Seismic network
The temporary deployment by Cambridge University consisted of three-component broad-band seismometers (Guralp 6TD, loaned by SEIS-UK) in and around the Krafla caldera.The temporary network at Krafla was deployed in July 2009 and remained at similar density within the Krafla caldera until September 2011 (Fig. 1A).Data from this period (July 2009 to September 2011) were recorded at 50, 100 or 200 samples per second (sps), depending on the time of year and the memory capacity of the instruments.We analysed waveform data from this array as well as nearby stations in the permanent South Iceland Lowland (SIL) seismic monitoring network operated by the Icelandic Meteorological Office, with data downsampled to 50 sps.

Earthquake locations
An initial catalogue of automatic earthquake locations was obtained using the Coalescence Microseismic Mapping code (Drew et al. 2013), from which the largest earthquakes within a cluster about 1.4 km NNE of the 1975-1989 inflation centre (derived from repeated ground levelling during the Krafla rifting episode; Björnsson 1985;Björnsson et al. 1985) were selected for manual refinement.P-and S-wave arrival times were picked manually for these events, which were located using NonLinLoc (Lomax 2001) and then relatively relocated using HypoDD (Waldhauser & Ellsworth 2000) to obtain final hypocentre solutions (Fig. 1B).
In order to test the robustness of our hypocentre locations and moment tensor solutions under different velocity model assumptions, three different velocity models (Arnott 1990;Arnott & Foulger 1994;Brandsdóttir et al. 1997) as well as an average velocity model, were trialled (see Table S1, Supporting Information).An incorrect velocity model may distort the station positions on the focal sphere, thus affecting the moment tensor solution.
The FIRE refraction profile provided the best constraint on the velocity structure within the Krafla region (Brandsdóttir et al. 1997) revealing a high-velocity mid-crustal dome beneath Krafla caldera, with two narrow (<100 m) low-velocity zones at shallower depths.Earlier velocity models used by Arnott (1990) and Arnott & Foulger (1994) were based on local earthquake tomography.Since the earthquakes are all located <3 km b.s.l., irrespective of which velocity model we use, variation in the shallow velocity structure is not significant in contributing to inferred hypocentre locations.The final hypocentre locations were obtained using the High Velocity Model (HVM), adapted from Brandsdóttir et al. (1997).
The V p /V s ratio is highly variable in the Krafla region (Schuler et al. 2015).Hypocentral locations using V p /V s values between 1.70 and 1.86 were tested, with the 1.76-1.78range providing the smallest misfit residuals and the most internally consistent models.This result is consistent with V p /V s ratios obtained from Wadati plots of earthquakes in this study and from seismic refraction data (Brandsdóttir et al. 1997).

Moment tensor inversion and decomposition
We used a Bayesian approach for moment tensor source inversion, which allows rigorous inclusion of uncertainties, such as location uncertainty, in the resultant probability density function (PDF) using marginalization (Sivia 1996, section 1.3;Pugh et al. 2016).The inversion approach determines the probability of a given moment tensor producing the observed data.We hand-picked the P-wave polarities, but could not use amplitudes for these relatively small events in the source inversion because they were difficult to measure accurately in the presence of noise.
The inversion was run twice, initially constrained to the DC space and then over the full range of possible moment tensor solutions, allowing non-DC components to be constrained.A uniform random sampling search algorithm was used to generate random samples of the source space.Each sample was tested against the data, and the probability of that sample producing the observed data was evaluated.The resultant PDF gives the maximum probability solution along with a range of possible solutions and their associated probabilities.Application of Bayes theory requires the inclusion of an a priori PDF.In this case, a uniform a priori PDF was used, signifying that no prior information about the mechanisms was assumed.

R E S U LT S
A total of 124 events from 2009 July to 2011 September were relatively relocated.Most lie in a cluster 200-500 m southeast of station K070 at 1-3 km b.s.l.(Fig. 1B), with smaller clusters due west.The epicentres are located within the central section of the fissure swarm and follow the northern margin of the shear wave attenuation zones delineated by Einarsson (1978) below 3 km depth.The final hypocentre locations were obtained using the HVM; adapted from Brandsdóttir et al. (1997).We use a V p /V s ratio of 1.78 for the final locations, although, as will be demonstrated, changing the velocity model and the V p /V s ratio within reasonable bounds does not greatly affect the locations or the moment tensor solutions.All locations are tabulated in the Supporting Information.

Moment tensor solutions
The moment tensors of all events have been calculated by inverting picked P-wave polarities.Within the main cluster, one event (Event A, 11:07 on 2009 November 11, Fig. 1B) has dilatational P-wave first-motion arrivals at all stations that could be picked unambiguously.The dilatational arrivals of Event A comprehensively cover the focal sphere (Fig. 2A), which indicates that the event has a large isotropic component.Such events are fairly rare, with only a few examples reported previously (Imai 1982;Foulger et al. 1989;Ford et al. 2009).Hence, this event merits particular scrutiny, to assess the stability of the solution and to ensure that the result is well constrained.Some other non-DC events were also found in the same cluster, although they are less well constrained (fewer clear first arrivals) and therefore are not discussed further.
Moment tensor solutions are plotted on a fundamental eigenvalue lune (Tape & Tape 2012a,b), which can be used to identify possible source types of an event and compare to other events.The moment tensor PDF for Event A indicates that this microearthquake was generated by a closing crack at depth (Fig. 2B).

Assessing the stability of the moment tensor solutions
In order to assess the stability of the closing crack solution, the moment tensor was computed for each velocity model and for a range of V p /V s ratios.These are the main uncertainties in earthquake location, and hence this is where inaccuracies in the location and take-off angles may be introduced.The moment tensor PDF using different velocity models and V p /V s ratios for Event A is shown in Fig. S1 in the Supporting Information.The solutions all plot in the same well-constrained region of the lune, indicating a stable closing crack solution.
A DC event within the same cluster, spatially close to Event A, (Event B 01:11, 2010 July 15, Fig. 2B), can also be shown to have a stable moment tensor solution.The DC events are spatially and temporally close to the non-DC events, therefore excluding systematic errors in take-off angles or seismometer polarity being the cause of a non-DC solution.
Taking into account location uncertainties, the probability of Event A fitting a non-DC closing crack solution is 0.98.The probability of Event A fitting a DC solution is less than 0.002.The probability of Event A fitting a CLVD (compensated linear vector dipole) solution is 0.01.In contrast, Event B has a 0.81 probability of fitting a DC solution and only 0.19 probability of having a non-DC component (Pugh et al. 2016).Event B has vertical and subhorizontal nodal planes, which is unusual for tectonic earthquakes.It is likely that the vertical nodal plane is the fault plane as it is aligned along the rift fabric.Collapse along a vertical plane is likely to be facilitated by a void space beneath it, although in this case, there is no recognized non-DC component in Event B itself.

D I S C U S S I O N
The highly implosive (closing crack) non-DC Event A is of interest because only a few examples of such events have been reported previously.Some skepticism of these results is appropriate, since such implosive events are rare within the literature and the method of faulting is likely to be unusual.Therefore, it is useful to consider other possible interpretations of the data presented.
One possible interpretation of the data is that the highly non-DC event is a result of multiple events occurring within a short space of time, and that the interference of waveforms creates the appearance of a single, non-DC event (Frohlich 1994).This has been suggested as a possible mechanism for generating highly non-DC (isotropic explosive) events around the Bár arbunga central volcano in NE-Iceland (Nettles & Ekström 1998).Visual comparison of waveforms from the non-DC event in Krafla with those of a DC event shows no obvious evidence of interference of waveforms (e.g.no significant or consistent secondary arrivals, and traces for the different events look similar in character, see Fig. S3 in the Supporting Information).So interference of multiple events is unlikely to be the cause of the apparent non-DC event.
A closing crack source mechanism generating Event A is favoured for several reasons.Melt and high-temperature geothermal fluid are present at depth, feeding the steam-dominated geothermal field at the surface (Gudmundsson & Arnórsson 2002).A region of low V p /V s ratio in vicinity of the cluster (Fig. 1C) has been hypothesized to be due to the presence of super-heated steam overlying rhyolitic melt (Schuler et al. 2015).The events are located at relatively shallow depths, in the vicinity of a large normal fault which was activated repeatedly by extension and normal faulting during the rifting episode.This means that it is possible for open cracks to exist (Frohlich 1994).They may then collapse due to a change in the hydrostatic pressure relative to the lithostatic pressure.From the moment tensor analysis conducted on the data, the orientation of these closing cracks is not well constrained, although this is to be expected from an event with only dilatational arrivals.Based on the volcanic-event classification of McNutt (2005), high-frequency events with clear onsets of P and S phases are classified as volcano tectonic (VT), whereas low-frequency events with emergent P and lack of S phases are classified as long period (LP).It has been previously hypothesized by Waite et al. (2008) that such implosive non-DC events may be generated by steam-filled cracks at depth, although the event presented by Waite et al. is an LP event (0.25-2 Hz), whereas the Krafla earthquakes are VT events with a frequency of 4-25 Hz and a dominant peak around 10 Hz at closest stations (Figs S3 and S4,Supporting Information).A CLVD solution may also fit Event A, although this would be a poor physical explanation for this earthquake, because a closing crack mechanism provides the most reasonable explanation given the present deflation of the Krafla magma chamber (Sturkell et al. 2008).
The moment tensor solutions clearly show the event plotting towards the solution for a closing crack.We are thus confident that Event A is a true non-DC event.Hence, unusual earthquakes in Krafla and other regions of microseismicity worldwide should not be discarded as spurious events: instead such events should be analysed as they may provide information about processes occurring at depth.
Microearthquakes have been used to track the propagation of dykes from volcanic systems in Iceland (e.g.White et al. 2011;Tarasewicz et al. 2012Tarasewicz et al. , 2014) ) and worldwide (e.g.Belachew et al. 2012;Falsaperla & Neri 2015).Reverse and normal focal mechanisms are documented around the tip of these propagating dykes.The intrusion of dykes will create a change in the rock volume, hence it is likely that opening (or closing) crack non-DC events may also be observed during propagation episodes.If other geodetic methods were available to monitor regional movements, such as inflation or deflation, observations of opening or closing crack microearthquakes could indicate volcanic activity and magma movement at depth.

C O N C L U S I O N S
A dense seismic network in Krafla has detected an example of a non-DC microearthquake.All first arrivals from this earthquake show dilatational polarities, with comprehensive coverage of the focal sphere.This event is located spatially close to a well-constrained DC event.Uncertainties in the velocity model, the V p /V s ratio and the location are tested for both events using a Bayesian approach.It is shown that both the DC and non-DC events are stable for the full range of uncertainties.The presence of a DC event indicates that the non-DC event is not due to a spurious seismic network or location method.The non-DC event is shown to have a large closing crack component.This is likely to arise from the current deflation of the magma chamber, occurring following the spreading episode in 1974-1985.

A C K N O W L E D G E M E N T S
The seismometers were borrowed from SEIS-UK under loans 891 and 914, with funding from the Natural Environment Research Council, grant NE/H025006/1.Juerg Schuler kindly provided data for the background V p /V s ratios shown in Fig. 1C.The Generic Mapping Tools (Wessel & Smith 1998) were used to produce the figures.

Figure 1 .
Figure 1.(A) Location map of the Krafla volcanic system.Caldera faults are marked in black, shear wave attenuations zones (Einarsson 1978) in purple, faults and fissures in brown and eruption sites in red (adapted from Hjartardóttir et al. 2015).Blue/white symbols show locations of seismometers: blue triangles are temporary stations operating from 2009 July, white triangles temporary seismometers operational from 2010 July to 2011 September and blue stars Iceland Meteorological Office permanent stations.Black box shows extent of (B).Inset shows location on a tectonic map of Iceland (Einarsson & Saemundsson 1987) with fissure swarms (yellow) and volcanic centres.(B) Map and orthogonal cross-sections of earthquakes discussed in the text.Events are colour-coded according to source mechanism.Red triangles show the location of the temporary seismometers.Event A is the closing crack event and Event B a reference double-couple event.Colour contours on the cross-section are modelled V p /V s ratios from Schuler et al. (2015), rotated by 8 • to align with orthogonal cross-sections.The low V p /V s ratio region is indicative of the presence of superheated steam above rhyolite melt.

Figure 2 .
Figure 2. (A) Lower hemisphere equal area projection of station locations on the focal sphere and associated vertical-component waveforms for the closing crack Event A. Waveforms are bandpass filtered 2-20 Hz and plotted with the same vertical scale.The horizontal scale is in units of seconds.All first arrivals that can be picked confidently have downward first motions which indicates that this event has a large non-DC component.(B) Moment tensor solutions of two well-constrained events.Event A is a closing crack event and Event B a double-couple solution, spatially close (<200 m) to Event A. The first column shows the variation in ray paths on the focal sphere due to the location uncertainty.Blue indicates negative polarity and red positive polarity, while green indicates no measured polarity due to unclear first arrivals.The second column shows the distribution of amplitudes and polarities for the best-fitting moment tensor solutions, with double-couple solution marked by black lines: blue shows negative arrivals shaded by amplitude, red shows positive arrivals shaded by amplitude, the positions of the stations on the focal sphere correspond to the maximum-likelihood NonLinLoc location.The third column shows the source PDF plotted on the fundamental eigenvalue lune, with blue corresponding to low probability and red to high probability.DC: double couple, CLVD: compensated linear vector dipole, TC+: tensile crack opening, TC−: tensile crack closing.No double-couple solution can fit the Event A data.