Ring Fault Slip Reversal at Bárðarbunga Volcano, Iceland: Seismicity During Caldera Collapse and Re‐Inflation 2014–2018

Microearthquakes reveal the kinematics of the Bárðarbunga caldera ring fault; during the 2014–2015 rifting event and gradual caldera collapse, and its subsequent, ongoing re‐inflation. Tightly constrained focal mechanisms during re‐inflation have reversed phase arrival polarities from events during the caldera collapse. Thus, the inner side of the steeply dipping northern caldera faults (averaging 81 ± 8°) has been moving upwards during the post‐eruptive period. Both precise relative relocations of the seismicity and fault plane solutions confirm that this is due to slip reversal on the same ring fault structure.

the caldera, seismicity is concentrated on the northern rim, supporting the inference from geodetic observations that the collapse was highly asymmetric. Ágústsdóttir et al. (2019) report normal faulting during the eruption on steeply inward-dipping faults, striking sub-parallel to the caldera rim, interpreted as representing a combination of piecemeal and trapdoor style collapse mechanisms (Acocella, 2007;Amelung et al., 2000).
There have been several reports that the polarity of fault motion changed after the eruption ended. Cross-correlation of earthquake waveforms recorded on a single seismic station shows a reversal in polarity just 2 months after the eruption ended (Jónsdóttir et al., 2017). Source parameter inversions derived by fitting regionally recorded earthquake waveforms show a similar result (Rodríguez-Cardozo et al., 2019). Between the end of the eruption and September 2020, the Global Centroid Moment Tensor (gCMT) project (Dziewonski et al., 1981;Ekström et al., 2012) reported six large-magnitude events (>M w 4.7) in the Bárðarbunga area, which all show reversed polarity compared to those during the eruption ( Figure S1 in Supporting Information S1). Geodetic data also provide evidence of caldera inflation initiating shortly after the eruption ceased. Li et al. (2021) tested different models to evaluate the roles of both renewed inflow of magma and viscoelastic relaxation following magma Focal mechanisms for six eruptive events, all normal faulting, in lower hemisphere projection in map view and rear hemisphere projection in cross sections. Triangles mark seismometers. Estimated trace of inner caldera fault (dashed) is from InSAR observations (Gudmundsson et al., 2016). Caldera rim (hatched) and bedrock surface (brown line) are from radio echo sounding (Björnsson & Einarsson, 1990). Blue line shows ice surface. Cross sections have no vertical exaggeration. withdrawal, concluding that both processes are likely important, with a combination of magma inflow to a sill at 10 km depth and slip on the caldera ring fault providing a good fit to the observations. However, these studies leave several open questions, most importantly whether this reversal in deformation represents a mirror image of that during the eruption. The limited hypocenter resolution in the aforementioned studies makes it impossible to answer this definitively. Furthermore, the geometry of the caldera ring fault or fault zone remains an open question (Ágústsdóttir et al., 2019), and this continued seismicity provides an opportunity to investigate this further. Throughout this paper, we refer to the caldera ring fault, but we discuss later the possibility that it could be a more distributed zone of faulting.
We here focus on two co-located large earthquakes on the northern caldera rim, and their associated fore-and aftershocks; one in 2014, during the eruption, and one in 2017, during the subsequent reinflation. The 2017 earthquake was chosen because the seismometer array was particularly dense at this time, including several temporary seismometers on the ice cap close to the caldera (Figure 1), providing more precise locations and tightly constrained fault plane solutions. The two mainshocks are in the same location within uncertainties, and apparently result from movement on the same fault, though the eruption (deflation) event has a normal faulting mechanism and the re-inflation event a reverse faulting mechanism. We also present hypocenter and fault plane solutions for fore-and aftershocks to each teleseismically-detected earthquake in Bárðarbunga between 2015 and 2018, to investigate the reactivation of caldera faulting after the eruption.

Seismic Data
Continuous seismic data were recorded on a dense local network of three-component broadband seismometers operated by Cambridge University, supplemented by seismometers from the Icelandic Meteorological Office (IMO), the British Geological Survey and University College Dublin ( Figure S2 in Supporting Information S1). Up to 72 seismometers were recording at once. They lie 2-100 km from the caldera, with 8 stations within 20 km and 44 within 50 km. This allows the recording of a wide range of takeoff angles for focal mechanism inversion. We studied syn-eruption events from the catalog of Ágústsdóttir et al. (2019), and used QuakeMigrate  to detect and locate earthquakes within 2 days either side of those reported in the gCMT catalog after the eruption ended. We use the same regional one-dimensional velocity model, from Ágústsdóttir et al. (2016). Earthquake magnitudes are from the gCMT catalog (M w ) where available, or from the IMO's Southern Iceland Lowlands network catalog (M lw ) where not (Jónasson et al., 2021;Rögnvaldsson & Slunga, 1993).
All events were analyzed manually, by picking P-and S-wave arrival times and assigning P-wave polarities using the Snuffler toolbox (Heimann et al., 2017). Waveform filtering was kept to a minimum; if the arrival remained unclear with a highpass filter above 4 Hz and lowpass below 20 Hz, no pick was made. Polarity was assigned with no filtering. Hypocenter locations were calculated using NonLinLoc (Lomax et al., 2000). Focal mechanisms were calculated using MTfit (Pugh & White, 2018). We found that good fits to the data were achieved for all events with inversions constrained to the double-couple (DC) moment tensor space, despite including observations from an unusually dense local seismic network.
Hypocenters were refined by relative relocation using hypoDD using the method of singular value decomposition (Waldhauser & Ellsworth, 2000). Absolute locations have a mean uncertainty of 0.59 km laterally and 1.16 km in depth, whereas relative relocations are an order of magnitude better, with a mean uncertainty of 38 m laterally and 129 m in depth. Depth is generally the least well constrained parameter, largely due to the lack of stations directly above the events. Improvement in the network density, particularly during 2017 when BARS and BARE were active on the glacier surface ( Figure 1; Figure S2 in Supporting Information S1), improved both absolute and relative location constraints for the 2017 earthquakes.
The Cambridge network yields an earthquake catalog with a magnitude of completeness M C ∼ 1.2 in this area . The gCMT catalog has an average M C ∼ 5.0 (Ekström et al., 2012).

Results
Thirty events were manually analyzed in this study, of which 24 show approximately east-west striking nodal planes and sub-vertical dip-slip faulting. The sense of slip is reversed between those during the eruptive period (normal faulting, n = 6), and post-eruptive period (reverse faulting, n = 18). Source parameters for each event are displayed in Table S1 in Supporting Information S1, and full focal spheres in Figures S3 and S4 in Supporting Information S1.
Six events show a large shift to shallower depth when relatively relocated (red boxed events in Figures S3 and S4 in Supporting Information S1). Their polarity picks are identical to those for events with better constrained locations, indicating that their different focal mechanisms are caused by differing take-off angles, uncertainties in which are not explicitly included in the MTfit inversion. These events are not discussed further in this paper.

Eruption Period (August 2014-February 2015)
Around 4,000 events, extending down to 6 km depth below sea level (bsl), were automatically detected and located within the caldera during the eruptive period by Ágústsdóttir et al. (2019). Six events from the eruptive period produced reliable fault plane solutions: a M w 5.1 event on 06 September 2014, with particularly good azimuthal station coverage, alongside the five next-largest events in the following two days ( Figure 1). The six relatively relocated hypocenters lie on or near the northern caldera rim, and between 2 and 4 km depth bsl (since the caldera rises to 2,009 m above sea level, they are 4-6 km below the surface).
The six well-constrained focal mechanisms each show a steeply dipping nodal plane striking between 078° and 097°, sub-parallel to the northern caldera rim, and dipping between 65° and 72°. They all indicate downwards motion on the inner side of the caldera ring fault. In cross-section, rear hemisphere projection focal mechanisms demonstrate this clearly, with one nodal plane consistently dipping steeply to the south, into the caldera ( Figure 1); this nodal plane is interpreted as the fault plane, fitting best with models of caldera collapse (Acocella, 2007). Section 4 provides further discussion on faulting and geometry.

Comparison of Co-Located Earthquakes During and After the Eruption
Focal mechanisms and first-motion waveforms for a M w 5.1 event in September 2014, during the eruption, and a M w 4.7 in September 2017, during the re-inflation period, highlight the most significant result of this study . Both focal mechanisms are tightly constrained, sharing similar fault planes striking approximately east-west. They demonstrate an exact polarity switch from normal to reverse faulting. First-motion polarities at all stations common to both events are reversed between the two earthquakes. The 2017 event has better nodal plane constraint due to better azimuthal station coverage, producing a smaller spread of possible fault plane solutions, particularly for the auxiliary plane (which defines the slip vector). Pressure (P) and Tension (T) axes on polar stereonets for the 6 eruptive events and 18 post-eruptive events (Figures 2e and 2f) demonstrate that the reversal in polarity between the eruptive and post-eruptive periods is consistent across all events we studied.

Post-Eruption Period (March 2015-December 2018)
Six post-eruptive earthquakes within the caldera were detected teleseismically and are included in the gCMT catalog. We investigated four of these that occurred between 2017 and 2018, and microearthquakes that occurred within two days before and after each (Figure 3).
Eighteen well-constrained focal mechanisms were obtained. As for the eruptive period, the post-eruptive focal mechanisms all show a steeply dipping nodal plane striking approximately east-west parallel to the caldera rim, again interpreted as the fault plane. Strikes range from 065° to 104°. In cross-section, rear hemisphere projection focal mechanisms clearly show the steep fault planes, with an average dip of 81 ± 8° southwards into the caldera.
All 18 well-constrained focal mechanisms show upwards motion on the inner side of the caldera ring fault, indicating reactivation of the caldera ring fault, which dropped down during the eruptive phase. Considering only the 16 coplanar events close to the inner caldera ring fault (box, Figure 3a), 14 of these show inward dipping fault planes between 70° and 88°. The remaining two show a steep outward dip, although slightly poorer constraint on their dips compared to others in this study means they could also be inward dipping within uncertainty ( Figure S4 in Supporting Information S1 shows full focal mechanism solutions).

Fore-and Aftershocks to the September 2017 M w 4.7 Event
The sequence of earthquakes around the M w 4.7 event on 07 September 2017 was examined in more detail, taking advantage of the excellent station coverage available at that time; ice stations BARE and BARS were active, as well as nearby KIS, DYN, and DJK, considerably improving coverage both in proximity and azimuth compared to the 2014 events. Most importantly, KIS is <2 km from the northern caldera rim, greatly enhancing location constraint, particularly in depth. All events with M lw > 1 between 06 September 2017 and 08 September 2017 were manually located, totaling 27 events ( Figure S6 in Supporting Information S1). This large earthquake and its associated events enabled us to calculate tightly constrained focal mechanisms for the smaller events, which is significant, as the point-source approximation on which our methods rely is more appropriate for smaller events. The focal mechanisms for the small events are consistent both with one another and with the larger events, allowing us to estimate the extent of the fault plane which slipped during the large event, discussed in Section 4.

Fault Plane Geometry From Relative Relocations
During the post-eruptive period, seismic activity remained concentrated in the north of the caldera, though the number and average magnitude of events are much smaller than during the eruption (Ágústsdóttir et al., 2019). Epicentral accuracies increased compared to the eruptive period, due to improved network geometry. Depth uncertainty remains similar throughout both periods. Most hypocenters relatively relocate within 0-4 km depth  Figure S2 in Supporting Information S1. Pressure (P, red) and Tension (T, blue) axes for: (e) All eruptive events; and (f) all post-eruptive events analyzed in this study. b.s.l. (2-6 km below the surface). Relative relocation collapses the hypocenter locations into a cluster just a few hundred meters across.
We use the relative relocations as a second method to investigate the caldera ring fault geometry. A range of best-fitting planes through the hypocenters was fitted using jack-knife testing (Tukey, 1958). All are sub-vertical (with mean dip of 88 ± 2° northwards), displayed as dashed lines on Figure 3b ( Figure S5 in Supporting Information S1 shows this without interpretation).

A Note on Double-Couple Versus Compensated Linear Vector Dipole Moment Tensor Solutions
All well-constrained focal mechanisms in this study could be fit by purely DC fault plane solutions, indicating slip on planar faults. This contrasts with the mechanisms for the larger events published in the gCMT catalog, which contain significant compensated linear vector dipole (CLVD) components, and possible volumetric contributions. Ágústsdóttir et al. (2019) also found most focal mechanisms to be purely DC, though a few required a CLVD component. However, the estimated fault rupture area for most earthquakes in this study (discussed further in Section 4) may be sufficiently small that it does not involve rupture of a significantly curved surface and so, in this case, quasi-planar fault slip would mean that the sources are likely to be DC. This is supported by the results of Rodríguez-Cardozo et al. (2021), who find that moment tensor solutions for earthquakes during the caldera collapse show larger CLVD components for larger magnitude events. Similarly, Rösler and Stein (2022) find, in a global review, that non-DC components, particularly for smaller earthquakes, are often artifacts of the inversion process.

Caldera Ring Fault Slip Reversal
The most significant result of this study is the mirror-image reversal of earthquake focal mechanisms between the eruptive and post-eruptive periods (Figure 2), indicating downward motion on the south side of the caldera ring fault during the eruption, and upwards motion during the post-eruptive period. This polarity reversal is also observed in moment tensor solutions from the gCMT catalog and suggested by results reported by Ágústsdóttir et al. (2019) (Figure S1 in Supporting Information S1), and corresponds with InSAR and GPS measurements from Li et al. (2021), indicating that the caldera is now re-inflating. The similarity of hypocenter locations and fault plane orientations between the two periods suggests that the mechanism of reversal is simple re-activation of the same fault or group of faults around the caldera rim.
This reversal is partly attributed to re-inflation of the magma storage region beneath the caldera (Li et al., 2021). The subsidence during the collapse period likely occurred as magma exited the crustal reservoir beneath the caldera, leading to a pressure drop and causing the roof to collapse (Gudmundsson et al., 2016). After the eruption ended, the caldera motion reversed to inflation, partly due to magma re-intruding into the crust beneath the caldera and exerting upwards pressure on the reservoir roof. Li et al. (2021) argue that viscoelastic rebound following magma withdrawal likely also contributes to the ongoing reinflation. It is believed that reversal began almost immediately post-eruption (Jónsdóttir et al., 2017;Li et al., 2021;Rodríguez-Cardozo et al., 2019), possibly as early as July 2015 (Ágústsdóttir et al., 2019). The delay until 2017 before the first teleseismically observed, large-magnitude earthquakes were detected likely results from the time taken to build critical stress before rupture after the sense of deformation reversed.
During the eruption, the caldera collapse was highly asymmetric, both in geodetic subsidence and in seismic moment release (Ágústsdóttir et al., 2019;Gudmundsson et al., 2016). Though these observations don't provide a complete picture of the ratio of seismicity and deformation, we tentatively interpret this as an indication that the asymmetric, trapdoor style mechanism continues, with the majority of fault slip occurring at the northern caldera rim.

Dip
The north-south cross-section in Figure 3b provides detailed insight into the geometry of the northern ring fault. Hypocenter locations in this study mostly do not lie on the topographic outer caldera rim, but are located around 1-2 km closer to the center of the caldera. Many of them align with the trace of an inner caldera ring fault, represented by the dashed line in Figure 3a and the arrow on the ice surface in Figure 3b, which was identified through InSAR imaging by Gudmundsson et al. (2016), capturing a M w 5.3 earthquake on 18 September 2014.
Relative relocations show hypocenters collapsed tightly into a band <2 km wide, perpendicular to the surface expression of the caldera rim, with sub-vertical best-fitting planes (dashed lines, Figure 3b). The most precise constraint on the fault planes for specific events, however, comes from individual focal mechanisms. All well-constrained focal mechanisms in this study show a steeply dipping fault plane sub-parallel to the caldera rim (rear hemisphere projections, Figure 3b), with a mean dip of 81° southwards into the caldera.
These individual fault plane dips differ from the findings of Gudmundsson et al. (2016), who report a steeply outward-dipping caldera fault on the north side. Rodríguez-Cardozo et al. (2021) find that where moment tensor solutions are close to DC (for smaller events), they indicate inward to sub-vertically dipping faulting on the northern caldera rim. Furthermore, when corrected for the geometry of the ring fault segment that ruptured, the CLVD moment tensor solutions for larger magnitude events also provide evidence for a sub-vertically or inward dipping ring fault.

Fault or Fault Zone
Comparison of the best-fitting planes through hypocenters with the individual fault plane solutions helps constrain the fault complexity. The improved hypocenter resolution of post-eruptive events in this study allows us to make this comparison more confidently than previous authors (Ágústsdóttir et al., 2019). The dip of the best-fitting planes through the hypocenters is sub-vertical, while almost all the individual fault planes from focal mechanism solutions dip more shallowly, between 70° and 88°, southwards into the caldera (Figure 3b). This discrepancy suggests that the zone of deformation at the northern caldera rim might best be described as a fault zone, comprising multiple imbricate faults.
The observations presented here support earlier interpretations that the northern boundary of the caldera is more complex than a single through-going fault plane. Holohan et al. (2011) andBranney (1995) suggest that faulting at calderas can be accommodated on a mixture of inward and outward dipping fault structures, which could explain the two outward dipping fault plane solutions identified in this study. However, notably, the dips of the best-fitting planes through earthquake hypocenters and of individual focal mechanisms are in closer agreement than previous estimates (Ágústsdóttir et al., 2019), and the locations and orientations of earthquake hypocenters and fault planes are strikingly consistent over the entire duration of the study (2014)(2015)(2016)(2017)(2018)). To answer definitively whether this zone of deformation consists of a single fault plane will require further improvement in the resolution of earthquake hypocenters and focal mechanisms, achievable by increasing station coverage, and by obtaining a better 3D velocity model.
In summary, we characterize the deformation at Bárðarbunga as dominantly comprising a mix of trapdoor and piston-type subsidence, with faulting most likely occurring within a narrow fault zone around the northern caldera rim.

Area of Fault Slip in an Individual Event
Detailed analysis of the M w 4.7 07 September 2017 earthquake and its aftershocks provides an opportunity to constrain the area of the fault plane that slipped in this event, and whether its size supports the reasoning discussed in Section 3.6 for the DC nature of our fault plane solutions. Aftershock distributions often outline the approximate area of the rupture plane for the larger magnitude "mainshock" earthquake. In Figure S6 in Supporting Information S1 we outline two possible, approximate rupture plane shapes based on the hypocenter distribution: one 3 × 1 km 2 rectangle and one 3 × 3 km 2 square. Calculations were then performed (Table S2 in Supporting Information S1) using the equation where 0 is seismic moment, is the shear modulus for an elastic solid, is fault rupture area and is average fault slip. We used combinations of possible values for each of the three variables to determine which estimates were most likely considering independent observations and additional information. Initial estimates for from the distribution of aftershocks in this study ( Figure S6 in Supporting Information S1) vary between 3 and 9 km 2 . Those for are from Bjarnason (2014), and references therein, who estimates 10 GPa at 1.5 km depth and 32.5 GPa at 5.0 km depth. Those for are from Gudmundsson et al. (2016), varying between 1.0 and 1.5 m.
We conclude that likely varies between 2 and 13 GPa. This result is consistent with estimates from seismicity and deformation during the 2018 Kīlauea caldera collapse (Anderson et al., 2019;Neal et al., 2019;Tepp et al., 2020). It is possible that the fault zone is either weaker with larger rupture area and/or slip, or the fault zone is stronger with smaller rupture area and/or slip. A GPS instrument in the center of the caldera recorded typically 0.3-0.5 m of subsidence associated with each earthquake of around M 5.2-5.4 occurring during the second half of September 2014 (Figure 3c of Gudmundsson et al., 2016). The low frequency of seismic signals from these earthquakes is also an indicator that the faults may be weak (Bjarnason, 2014). We conclude that the range of areas of fault zone are likely to be sufficiently small compared to the scale of the caldera not to require a significantly curved rupture plane, thus reinforcing the interpretation that the source of all but the largest caldera earthquakes is likely to be DC, as discussed in Section 3.6.

Conclusions
Motion on the caldera fault, or fault zone, reverses; the inner side of the caldera faults moves downwards during the eruptive period, and upwards during the post-eruptive period. This is reflected in waveforms from individual events and in the P and T axes of all events in this study. Comparing polarity observations from co-located events during and after the eruption, and precise relative relocations both indicate reversal of motion on the same ring fault structure.
Tightly constrained focal mechanisms show motion on steeply dipping faults on the northern caldera rim. Failure during individual earthquakes occurs on planes dipping slightly shallower than the best fitting plane through the event hypocenters, indicating that deformation occurs in a fault zone, though this result might also be influenced by un-modeled 3D velocity structure.