Background

The Reykjanes Peninsula in southwest Iceland is an oblique rift, and part of the North-American-Eurasian plate boundary (Fig. 1). It has volcanic systems with NE-SW trending fissure swarms, arranged en echelon (Sæmundsson et al. 2020). A series of N-S trending strike-slip faults, in particular in areas between the volcanic systems, release shear along the plate boundary (Pagli et al. 2003; Clifton and Kattenhorn 2006; Einarsson 2008; Einarsson et al. 2020; Sigmundsson et al. 2020).

Fig. 1
figure 1

Map of Reykjanes Peninsula displaying checked earthquakes ≥ M2 from 15 December 2019 to 4 August 2022. The location of the modeled dikes (dike 1 from Sigmundsson et al. 2022, and dike 2 and dike 3 from this study) is displayed by yellow (dike 1), green (dike 2), and blue (dike 3) colored lines at Fagradalsfjall. The location of the eruptive fissures from the 2021 (Pedersen et al. 2022; Hjartardóttir et al. 2023) and 2022 eruptions are marked as orange lines and the 2021 and 2022 lava flows as dark red–filled regions (from Pedersen et al. 2022 and Gunnarson et al. 2023 respectively). Lava flows from 900 to 1240 AD are displayed by the light red–filled regions (Jóhannesson and Sæmundsson 1998). Pre-existing fractures are displayed as dark gray lines (Clifton and Kattenhorn 2006). Eruptive fissures from 900 to 1240 AD as red dashed lines (modified from Clifton and Kattenhorn (2006) and Sæmundsson et al. (2010)). The yellow fissure swarms in the inset map are from Einarsson and Sæmundsson, 1987. The background digital elevation model is the IslandsDEM (LMI 2022)

Full size image

Since February 2021, diking events and eruptions have occurred at Mt. Fagradalsfjall (Fig. 1 and Fig. 2), which are a part of a major period of unrest on the Reykjanes Peninsula. The unrest began with an earthquake swarm in the Fagradalsfjall region between 15 and 20 December 2019 (Barsotti et al. 2023). Following a short quiescence, activity re-commenced on 21 January 2020, with a small swarm of earthquakes near Grindavík, about 10 km southwest of Fagradalsfjall (Flóvenz et al. 2022). Deformation in this area, indicative of inflation, was confirmed by movement on two GNSS (Global Navitation and Satellite System) stations in this region and on Sentinel-1 interferograms. Geodetic inversions of these observations suggested the inflation resulted from a pressure increase at a depth of about 3–5 km, in a sill-shaped area that can be modeled with opening of a near-horizontal rectangular dislocation, centered directly west of Mt. Thorbjörn (Geirsson et al. 2021; Flóvenz et al. 2022). This inflation was followed by two additional sill-type intrusions in a similar location, between 6 March–17 April and 15 May–22 July 2020, respectively (Barsotti et al. 2023; Cubuk-Sabuncu et al. 2021). The three intrusions near Mt Thorbjörn comprised a total volume change of about (9–11)×106 m3. In mid-July 2020, inflation was detected in the Krýsuvík volcanic system, approximately 10 km east of Fagradalsfjall. This episode of inflation lasted approximately 6 months. Geodetic inversions indicated the observed signal could be fit by an increase in pressure within a Mogi source (Mogi 1958) located at a depth of ~6 km (Parks et al. 2022). From December 2019 to January 2021, seismicity shifted along different regions across the Peninsula, in relation to continued unrest triggered by magma migration and release of tectonic stress (Geirsson et al. 2021).

Fig. 2
figure 2

Map of Fagradalsfjall region, Reykjanes Peninsula, displaying checked earthquakes ≥ M1 from 30 July to 3 August 2022. These represent 10% of all detected earthquakes ≥ M1 that occurred during this time period. Earthquakes ≥ M4 are displayed as colored stars. Horizontal GNSS displacements are shown by black arrows. The location of the modeled dike (final Dike 3 uniform opening model spanning entire intrusive period) is displayed by the red line. The location of the eruption which started on the 3 August 2022 is displayed by the red triangle. The approximate location of the plate boundary is marked by the magenta dashed line (Michalczewska et al. 2014). The diking region (used to select data for Fig. 3) is displayed by the green rectangle

Full size image

Intense earthquake swarms commenced on 24 February 2021, concentrated at both Fagradalsfjall and also extending across a 20-km segment along the plate boundary, with earthquakes up to Mw5.64. Deformation was detected at continuous GNSS geodetic stations, and subsequent Interferometric Sythethic Aperture Radar Analysis (InSAR) of Sentinel-1 data confirmed the observed deformation was primarily the result of dike emplacement (opening and shear slip), as well as shear slip along the plate boundary segment. Geodetic inversions indicated a ~9-km-long dike with a total intruded volume of around 3.4×107 m3 (Sigmundsson et al. 2022). The shear slip on the dike patches multiplied by their areas is 4.4×107 m3 (right-lateral) and shear slip on the plate boundary (PB) patches multiplied by their areas is 4.0×107 m3 (left-lateral) (geodetic moment divided by shear modulus, Table 1). These numbers show shearing contributed importantly to the deformation. When multiplied by a shear modulus of 15 GPa, the geodetic moment on the PB is 0.6×1018 Nm. The combined numerical values of the shearing moments for both the PB and the dike correspond to a geodetic moment of 1.3×1018 Nm. Seismic activity occurred both in the diking area and in areas to the west and east, due to triggered fault slip, releasing previously accumulated stress along the plate boundary. The estimated seismic moment release in earthquakes is ~1×1018 Nm (Sigmundsson et al. 2022). If the shear on the dike is considered to be mostly aseismic, then the earthquake moment should be equal to the moment of the slip on the plate boundary alone. However, some of the modeled slip on the dike (thus the moment assigned to the dike) may actually be occurring on PB-crossing N-S faults that also cross or reach the dike (off-dike fault slip). Thus, the moment on the PB alone (0.6×1018 Nm) is a lower bound and the moment of the PB+dike (1.3×1018 Nm) is an upper bound of what we expect to slip seismically. The actual seismic moment falls in between these two values. The stored tectonic stress was initially released most rapidly, followed by a decline in deformation and seismicity preceding the eruption onset, on 19 March 2021 in Geldingadalir at Mt. Fagradalsfjall. The eruption continued until 18 September 2021 and produced a lava field covering an area of 4.8 km2 with an extruded bulk volume of (150 ± 3)×106 m3 (Pedersen et al. 2022).

Table 1 Parameters of diking events at Fagradalsfjall between 2021 and 2022
Full size table

A second dike intrusion occurred beneath Fagradalsfjall between 21 and 27 December 2021. The intruded median volume from geodetic modeling was estimated at 19×106 m3, the median top depth of the dike at 2 km, and the initial medium magma inflow rate was ~22 m3/s. This intrusion did not culminate in an eruption. A fourth sill intrusion west of Mt. Thorbjörn, began forming at the end of April 2022, in a similar location to previous intrusions there. The intrusion lasted until the end of May 2022 and the estimated median volume change was 4×106 m3. Following several months of quiescence, a third dike began propagating beneath Fagradalsfjall on 30 July 2022, culminating in an effusive eruption on 3 August 2022. This intrusion and corresponding deformation and seismicity are the focus of this study.

Data and methods

Seismicity

Seismicity is recorded by the national seismic network of the Icelandic Meteorological Office (also referred to as the South Iceland Lowland (SIL), Veðurstofa Íslands (VI) network). Arrival times of P- and S-waves from earthquakes are automatically detected by the analysis system, which locates the events and calculates a preliminary magnitude (Böðvarsson et al. 1999; Jakobsdóttir 2008). The events are then manually reviewed and revised as required. However, during periods of high seismic activity, such as that analyzed here, only a fraction of the seismicity has been manually reviewed. The seismicity analysis (Fig. 3) is therefore based on automatic locations. During intense seismic swarms with large earthquakes, the automatic location system may saturate and background noise greatly increases, resulting in a significant decrease in sensitivity. For the intense period on Reykjanes peninsula analyzed here, the magnitude of completeness, which normally is 0.5 or better (Panzera et al. 2017), is estimated to be around M1. The seismicity rate analysis is therefore limited to M≥1 earthquakes (Fig. 3). The automatically estimated magnitudes are based on a local magnitude scale, ML, applied to decay of peak amplitudes with distance (Gudmundsson et al. 2006). These, however, tend to be increasingly underestimated for larger magnitudes. To improve the ML estimate, a corrective relationship between ML and MW was generated from a dataset of 206 events from Reykjanes peninsula, of magnitudes 2.8 ≤ ML ≤ 5.3. The correction, MW = ML + 0.2, was applied to all the events. The mechanism of the MW5.3 event was calculated by the SIL analysis software (Rögnvaldsson and Slunga 1993).

Fig. 3
figure 3

Timeline of earthquakes ≥ M1, within the diking area (see Fig. 2 for location) from a 21 to 27 December 2021 (second diking event) and b 30 July to 4 August 2022 (third diking event). Panels show (upper): number of earthquakes per hour (blue) and cumulative number of earthquakes (red) (middle) magnitude of earthquakes (black) and cumulative seismic moment versus time (red) (lower): 8-h solutions for GNSS stations KRIV, LISK, and FEFC, plotted along the direction of maximum motion during each dike intrusion. The LISK timeseries in a is incomplete due to power outages during this time period. The thin red vertical line in b shows the time of eruption onset on 3 August 2022

Full size image

Interferometric Synthetic Aperture Radar processing

Interferograms were generated using the Interferometric Synthetic Aperture Radar Scientific Computing Environment software (Rosen et al. 2012). The images were multilooked by 2×6 (azimuth look × range look) for Sentinel-1, 3×3 for COSMO-SkyMed, and 5×5 for TerraSAR-X, to increase the signal-to-noise ratio. Topography was corrected for using the IslandsDEMv1, a high-resolution DEM of Iceland made by the National Land Survey of Iceland, resampled to a 10-m grid (LMI 2022). Unwrapping was undertaken using Snaphu (Chen and Zebker 2002). Water/ocean and the 2021–2022 Fagradalsfjall lava flows were masked out to ensure only stable ground motion was used for modeling. Errors in the interferogams were estimated using the fitVariogram function in the GBIS software (Webster and Oliver 2007; Bagnardi and Hooper 2018). A part of the interferogram in a non-deforming area is selected to build a variance-covariance matrix that is then fit using an exponential function. Atmospheric perturbations are thus included in the uncertainty estimates, but not corrected for. Atmospheric perturbations in the interferograms are expected to be of different scale and shape than the signals related to the dike intrusions studied here. However, smaller deformation signals may be obscured, both by the larger events and by atmospheric effects.

Continuous Global Navigation Satellite System processing

The Reykjanes Peninsula GNSS data were analyzed using GAMIT 10.7 (Herring et al. 2018) combined with data from over 60 IGS (International GNSS Service, https://igs.org/) stations and IGS satellite orbit information. Four- and eight-hour GNSS solutions (“ultra-rapid”) were analyzed using a subset of the IGS stations, for which hourly data are available within 20 min of the hour. Each 8-h session contains the previous 8 h of Reykjanes Peninsula GNSS monitoring data and 24 h of Global Reference Station data. For key stations (FAFC and GONH), we estimate site positions splitting the 8-h data into two 4-h solutions. During the processing of both daily and sub-daily solutions, we solved for satellite orbits and earth rotation parameters, estimating atmospheric zenith delay every 2 h. We corrected for ocean loading using the FES2004 model (Lyard et al. 2006) and apply the IGS14 azimuth and elevation-dependent phase center variation model for all antennas (Schmid et al. 2016; Rebischung and Schmid 2016; Montenbruck et al. 2015). We estimate site positions for the daily solutions shortly after 00 UTC as well as three times a day for the sub-daily solutions, starting the processing shortly after 00, 08, and 16 UTC. We used GLOBK to combine solutions and estimated daily site positions in the ITRF2014/IGB14 (Altamimi et al. 2012) reference frame. For Continuous Global Navigation Satellite System (cGNSS) stations with sufficient length of observations, we estimated and removed linear, annual, and semi-annual terms for each GNSS station based on interseismic data collected prior to the beginning of unrest. For continuous stations installed during the unrest, the values estimated for the nearest continuous GNSS station were used. The site positions are published about 1.5 h after the hour.

Geodetic modeling

Geodetic inversions are undertaken utilizing both InSAR and GNSS observations to constrain the range of source parameters that can explain the observed deformation during the second and third diking events. Where possible, a combinatinon of ascending and descending interferograms were utilized. This improves the inversion result by reducing the non-uniqueness of the source parameters, particularly in terms of the latitude and longitude of the source. The geodetic inversions were carried out using a modified version of GBIS software (Bagnardi and Hooper 2018). This Bayesian inversion software generates final a posteriori probability distribution functions for all model parameters, using a Markov-chain Monte Carlo approach incorporating the Metropolis Hastings algorithm (Hastings 1970; Mosegaard and Tarantola 1995). In our initial inversions, we modeled the second and third dikes with a rectangular dislocation (Okada 1985) allowing for opening (D in Fig. 5b, Fig. 9b, and S10b), and the observed movement at the central axis of the plate boundary (PB) with a second dislocation, allowing for shear (PB in Fig. 9b and S10b). For the third diking event, we also modeled the surface displacements resulting from the MW5.3 tectonic earthquake as a separate rectangular dislocation, allowing for strike-slip and dip-slip motion (F in Fig. 5b). A subsidence signal to the west of Mt. Thorbjörn was initially modeled as a deflating point source (Mogi 1958) for simplicity (M in Fig. 5b). However, this subsidence may have been the result of normal faulting, so in the final dike 3 model we included an additional rectangular dislocation to account for this deformation. We assumed a Poisson’s ratio of 0.27. Wide bounds were defined a priori for each of the source parameters to enable them to vary. For these models, the joint probability distribution function for the various source parameters was obtained by running one million iterations. The uncertainties for individual parameters are reflected in their one-dimensional marginal posterior probability density functions, and trade-offs between parameters are reflected in the respective two-dimensional marginal probability density functions.

For our final inversion, we fixed the location of the dike and the PB to the preferred location from these inversions detailed above. We allowed for only opening on the dike and shear motion on the PB. We divided the dike and plate boundary segment into patches. For the second diking event, the dike patches are 500 × 500 m2 from 0 to 7 km. For the PB, the patches are 500 × 500 m2 in the upper 2 km and 1.0 × 1.0 km2 below a depth of 2 km. For the third diking event, the dike patches are 500 × 500 m2 in the upper 3 km and 1.0 × 1.0 km2 below a depth of 3 km. For the PB in this model, the patches are 750 × 750 m2 in the upper 3 km and 1.5 × 1.5 km2 below a depth of 3 km. We solved for the opening and slip on each patch of the dike and PB respectively. For these final models, the joint probability distribution function for the various source parameters was obtained after running two million iterations. Smoothing was not applied, as this is not required when using a Bayesian approach. The non-uniqueness is captured in the posterior probability distribution, and the median value of slip/opening reflects this, giving a smoother result than any individual slip/opening solution in the distribution.

Initial observations of dike intrusion July–August 2022 and eruption warning

A third dike intrusion at Fagradalsfjall began on 30 July 2022, with seismicity providing the first indications. The earthquake activity commenced at 08:29 UTC approximately 2 km northeast of the forthcoming eruption site. Over the next 3 h, a continuous swarm of small earthquakes (less than M1.5) migrated ~800 m westward at a depth of ~8 km. Around noon, the seismicity moved into the upper crust, migrating in the following few hours to ~3-km depth and propagating laterally ~1.5 km southwest towards the eruption site and ~1.5 km northeast towards Mt. Keilir (location shown on Fig. 6). As the seismicity moved into the upper crust, earthquake magnitudes also increased, up to greater than MW4. In the 24 h prior to the onset of seismicity, on the morning of 30 July, only 2 small earthquakes (less than M1.0) had been detected within the diking area, whereas in the following 24 h over 1700 earthquakes were detected and nearly 600 earthquakes with magnitudes greater than M1.0 (Fig. 3). The activity continued and strong earthquakes were felt in the capital area over several days, the largest of which was MW5.3 that occurred on 31 July outside the diking area, approximately 2 km to the NE of Grindavík (Fig. 2). Based on studies of earlier intrusive activity on the Reykjanes Peninsula (Sigmundsson et al. 2022), at the Krafla volcano (Einarsson and Brandsdottir 1980, Einarsson and Brandsdóttir 2021), and at the Bárðarbunga volcano (Sigmundsson et al. 2015), the timeline and spatial distribution of seismicity suggested it was resulting from a dike intrusion together with triggered seismic activity in nearby areas releasing stored tectonic stress. Deformation signals due to both a dike intrusion and tectonic earthquakes were anticipated.

To improve temporal resolution, the GNSS data from Reykjanes Peninsula were analyzed every 8 h, providing three datapoints per day with a 1.5-h delay. Nearby, GNSS data show a delayed response to the seismic swarm, with onset of deformation starting after 16:00 UTC on 30 July when the seismicity had already moved to within a few kilometers of the surface. The ground displacements were suspected to be a result of diking because of the observed horizontal displacements outward from Fagradalsfjall, at stations located at distances up to 15 km to the northwest and southeast (Fig. 2, Fig. 3, and Fig. S1).

Initial constraints on the deformation from InSAR were based on a Sentinel-1 image acquired on ascending Track 16, acquired at 18:59 UTC on the evening of 1 August. Due to download delays, the data was received at 07:25 UTC on 2 August and an interferogram rapidly formed with the previous image available in the same orbit on 20 July (12-day time span) (Fig. 4). Line-of-sight (LOS) changes were up to 16 cm, with the main observed signal consistent with a dike intrusion, confirming the prior interpretation based on GNSS observations and seismicity. Additional features in the interferogram relate to deformation associated with the MW5.3 earthquake (blue ellipse, Fig. 4(a)) and a localized signal west of Mt. Thorbjörn (blue arrow, Fig. 4(a)).

Fig. 4
figure 4

Sentinel-1 interferogram spanning 20 July to 1 August 2022, wrapped (a) and unwrapped (b). Blue ellipse in (a) shows deformation associated with the MW5.3 and MW4 earthquakes on 31 July 2022 and blue arrow, a localized signal, likely associated with earthquakes on 1 August 2022

Full size image

Once the initial interferogram and GNSS displacements were available, geodetic modeling was carried out the same day utilizing the data spanning this new intrusive period. A Bayesian inversion was undertaken utilizing both cGNSS and InSAR data (Sentinel-1 interferogram spanning 20/07/22 to 01/08/22) as the input. This first deformation model considered three separate sources embedded within an elastic halfspace: (i) opening of a dike, modeled as uniform opening on a rectangular plane (dislocation) with no slip (Okada 1985), (ii) a second dislocation allowing for strike-slip and dip-slip motion (with optimal parameter values: length 3.4 km, width 1.0 km, depth 1.2 km, dip −89.9°, strike 21°) to model the Mw5.3 earthquake near Grindavík, and (iii) a point pressure source (Mogi 1958) to model a small subsidence signal identified to the west of Mt. Thorbjörn. The subsidence was likely related to earthquakes that occurred in this area on 1 August, but a Mogi source was used to simplify the model as much as possible and optimize the run time. A best fit model (Fig. 5) indicated that the top depth of the dike was very shallow (median depth of 1.1 km with 95% confidence interval (CI) 1.0–1.3 km) with a volume change of 10.3 million m3 (CI 6.3–14 million m3) (over a 58-h period) and a high magma inflow rate (~49 m3/s, CI 30–67 m3/s). Median values of model parameters indicate a length of 2.7 km (CI 2.3–3.0 km), width of 2.9 km (CI 2.0–3.9 km), and strike of N36°E, extending from the northern edge of the 2021 lava field to the northwestern side of Mt. Kistufell (Fig. 6).

Fig. 5
figure 5

Preliminary geodetic model for Dike 3, using data covering the onset of diking until 1 August 2022. Panels show a input data: observed horizontal GNSS displacements (black arrows) and line-of sight change in the Sentinel ascending interferogram (T16 20/07/2022-01/08/2022) (color), b predicted ground displacements and LOS change, c residuals, and d distribution of the dike volume change in the series of models produced in the inversion. Red star in a shows location of eruption site on 3 August 2022. In b, red line shows modeled dike (D), black line is modeled fault (F), and white star is Mogi source (M)

Full size image
Fig. 6
figure 6

Projected surface location of the modeled dikes. The yellow line is dike 1 (intruded February to March 2021) from Sigmundsson et al. (2022), green line is modeled dike 2 (December 2021), initial dike 3 model (July–August 2022, using data up until 1 August 2022) prior to eruption onset is the dashed blue line, and final dike 3 model (July–August 2022) is the solid blue line. The location of the 2021 eruptive fissures are displayed as orange stars. The location of the eruptive fissure that opened on 3 August 2022 (~13:30 UTC) is shown by the orange line. The 2021 lava outline is displayed as the red-filled region with yellow outline (from Pedersen et al. 2022) and 2022 lava flow with blue outline (from Gunnarson et al. 2023). GNSS stations are marked as filled black circles

Full size image

In addition, the seismic activity had started to display indications that it was declining (Fig. 3), in a similar manner to that which occurred prior to the 2021 eruption, and which was precursory to that eruption (Sigmundsson et al. 2022). A meeting was held with Civil Protection, staff from the Icelandic Meteorological Office (IMO) and University of Iceland at 15:00 UTC on 2 August, to discuss the new intrusion and latest data/modeling results. That evening (at 17:49 UTC, 2 August), IMO issued a warning that the likelihood of a new eruption in the coming days was high. This was reported on the IMO website, in local media and on the evening news (IMO 2022).

An effusive eruption started the next day (3 August 2022) at approximately 13:30 UTC when a 375-m-long fissure opened to the NE of the 2021 eruptive craters. The projected surface location of the dike (from the optimal model) was within 49 m to the north and 110 m to the south of the fissure opening (Fig. 6). On 18 August 2022, the volcanic tremor began to steadily decrease until it stopped completely in the early hours of 21 August 2022. The activity at the crater ceased at the same time, with no visible extrusion of fresh lava thereafter. The eruption was declared ended on 21 August 2022.

Effects over a wider area and improved deformation models

A COSMO-SkyMed (CSK) SAR image was acquired in ascending orbit at 06:41 AM on 4 August 2022 and an 8-day interferogram was generated utilizing the previous image acquired in the same track (interferogram spanning 27 July to 4 August, Supplementary Information Fig. S2a and b). This showed a substantial amount of deformation related to both the dike intrusion and large earthquakes (one>M5 and 11>M4), which occurred during this period. Up to 29 cm of LOS motion was observed in relation to the dike intrusion. The interferogram indicated that the dike had not lengthened (after 1 August 2022). However, it was evident in this interferogram that movement had also occurred along the central axis of the plate boundary, similar to that observed during the February–March 2021 dike intrusion (red arrows on Fig. S2a).

On 8 August, reports were received of large fractures observed on Mount Hagafell (Fig. 2), located northeast of Grindavík (Fig. 7). Although it was suspected that these were the result of the MW5.3 earthquake that occurred on 31 July 2022, there was significant concern from the general public that a new eruptive fissure might open close to Grindavík. A meeting was called with Civil Protection on 9 August 2022 and additional geodetic modeling focused on this region was presented. This second deformation model utilized the ascending CSK interferogram spanning 27 July to 4 August 2022 as the input (Fig. S2a and b) and the surface displacements resulting from the tectonic earthquakes in this area were modeled as rectangular dislocations, allowing for strike-slip and dip-slip motion. We tested incorporating one or two dislocations in the modeling. The best-fit model indicated that movements occurred on two north-south striking fault planes to the east of Mount Hagafell (Fig. S3), with predominantly right-lateral strike-slip motion on the eastern fault and a combination of strike-slip and dip-slip on the western fault. However, the requirement of the second fault may in fact be unnecessary—possibly just compensating for additional movement along the PB in this region (not included in this particular model). The use of only one satellite viewing geometry here (one ascending interferogram) also affects the quality of the derived model parameters. Several fault parameters were not well constrained including the longitude location of fault 1 and the width of fault 2. There was also a residual signal to the east related to the dike intrusion at Fagradalsfjall.

Fig. 7
figure 7

Slump deposits mapped on the summit of Mount Hagafell. a Aerial photograph (from Loftmyndir ehf) of summit of Mount Hagafell draped on top of IslandsDEM (LMI 2022), showing the location of mapped slump deposits (red stars). b IslandsDEM (LMI 2022) showing location of mapped slump deposits on Mount Hagafell (red star) and mapped fractures north-east of Mount Hagafell (blue star), location of Mw5.3 earthquake on 31 July 2022 (yellow star) and additional Mw4 earthquake that occurred within 1.5 s (magenta star). Red line in b is the modeled dike 3 which was intruded during July–August 2022. c Example photo of slumping at Mount Hagafell

Full size image

The MW5.3 earthquake was preceded by a smaller foreshock 1.5 s earlier, about 600 m away. This closeness in time and space caused their waveforms to overlap, contaminating the mechanism and magnitude calculation of the foreshock. In the routine reviewing in the immediate aftermath of the events, the magnitude of the foreshock was overestimated by approximately a whole magnitude unit. Later, more careful inspection of the first event and comparison with the second event resulted in a MW4 estimate for the foreshock, therefore significantly decreasing its relevance in the deformation modeling. The MW5.3 earthquake is dominated by right-lateral strike-slip. A component of right-lateral strike-slip on north-south oriented faults as inferred is in line with what is expected from the tectonic setting; an array of north-south-oriented strike slip faults has been mapped on the Reykjanes Peninsula (Clifton and Kattenhorn 2006; Einarsson et al. 2023). The dip-slip component on the second fault is also consistent with extension across the fissure swarms, which is accommodated by normal faulting oblique to the plate boundary. There was no evidence from GNSS or seismicity for intrusive activity in this region.

A team from IMO visited the summit of Mount Hagafell on 10 August 2022, where slump features were observed. Soil had sunk up to 60 cm, in zones up to 4 m wide (Fig. 7), likely the result of collapse into older pre-existing fractures (Fig. 7) due to amplified shaking in this region during the MW5.3 earthquake. The team also visited the site of the modeled faults northeast of Mount Hagafell, where new fractures were observed (Fig. S4). The main fracture at this location was approximately 50 m long, with measured openings up to 40 cm and correlates with a NE-SW trending fracture system, observed clearly on the drone image. This is representative of en echelon arranged surface fracturing, associated with the north-south oriented right-lateral strike-slip faults, as inferred from geological studies of surface features of such faults on the Reykjanes Peninsula (Clifton and Kattenhorn 2006; Einarsson et al. 2023), and in agreement with the modeled faults (Fig. S3).

The Sentinel-1 acquisition on 11 August 2022 from Track 155 did not include data in the region of interest. An additional CSK interferogram was generated on 12 August 2022 (spanning 4–11 August) (Fig. 8). This interferogram displayed no deformation in the vicinity of Grindavík, confirming no intrusive activity in this area. However, the interferogram did display a small NE-SW trending subsidence signal (~−2 cm) in close proximity to the new dike intrusion at Fagradalsfjall. To determine the cause of this signal—whether related to atmospheric delay from the eruption plume or continued evolution of the dike—decomposition of ascending and descending interferograms into approximate east and up displacements (near-east and near-vertical) was undertaken (with same procedure as described by Drouin and Sigmundsson (2019)), covering two separate time periods (prior to and after the eruption onset, Fig. S5 and Fig. S6). The difference in the near-east and near-vertical displacements during the two time periods (Fig. S7) indicates a small amount of continued inflation of the dike. This is in agreement with the horizontal movement on GNSS station FAFC directly to the east of the dike which displayed continued easterly motion until mid August (Fig. S8).

Fig. 8
figure 8

COSMO-SkyMed (CSK) ascending interferogram spanning 4 to 11 August 2022, wrapped (a) and unwrapped (b)

Full size image

In order to provide stronger constraints on the parameters of the dike intrusion (immediately prior to eruption onset) and for an improved understanding of other processes that occurred during this diking episode, a series of additional geodetic models were undertaken with varying complexity. A third deformation model, undertaken using two Sentinel-1 interferograms (ascending T16 20/07/2022-13/08/2022 and descending T155 30/072022-23/08/20222, Fig. S2) and GNSS data up until 4 August 2022, was run, incorporating the dike (with uniform opening only), a fault to the N-E of Grindavík to account for the Mw5.3 earthquake on 31 July 2022 and a second fault to the west of Mt. Thorbjörn. Both faults were allowed only to slip, without opening. Although two earthquakes occurred on 31 July 2022, in the region to the N-E of Grindavík, we only model the larger Mw5.3 earthquake here, as the smaller Mw4 earthquake is not considered to contribute significantly to the broader deformation field modeled here. As shearing was also evident on the new interferograms (both in the vicinity of the dike and the plate boundary axis, Fig. S2), we also included a vertical dislocation (Okada 1985) allowing for slip (strike-slip motion only) along the central axis of the plate boundary. The optimal dike model was very similar to that obtained from the earlier model, with a dike median length of 2.7 km, a top depth of 1 km, and a strike of ~N35°E, but with significantly more opening (median value of 2.8 m). The median length of the PB segment was 13.9 km with a width of 1.6 km, strike of N63°E and shear of 0.19 m. The parameters for each source in this model were well resolved after one million iterations. This model also suggests that the dike did not lengthen after 1 August but did continue to inflate. This may be representative of the bend in the eastern displacements from GNSS observations at FAFC (Fig. S1 and Fig. S8). Based on the progression of the seismicity, at the onset of this intrusion at ~08:30 UTC on 30 July 2022, the magma was initially at a depth of ~7–8 km in the area directly west of FAFC. During the first 4 h, the dike propagated both vertically and laterally to the SW and in the following 3 h propagated up to ~3-km depth, at which time FAFC started moving to the SW (Fig. S8, East component). Between 21:00 UTC on 30 July and 01:00 UTC on 31 July 2021, the dike started propagating laterally to the NW and seismic activity was mainly concentrated there through to 1 August. This NW propagation and continued inflation of the dike, may be the cause of the rotation of the signal at FAFC from SW to SE commencing in the evening of 31 July (Fig. S8, East component), and continuing there through to 3 August. Triggering of a MW4.7 earthquake, ~1 km NW of the station on 1 August may be indicative of dike widening. The effect of the earthquake itself, however, on the movement of FAFC is insignificant compared to that of the dike. It should also be noted, there appears to be a delay between the migration of the seismicity to the NW and the rotation observed at FAFC. Analyzing these signals at a site very close to the intrusive activity can be complicated; the propagating dike tip may not always produce seismicity (Smittarello et al. 2019) and when it does, it may not relate to the exact location of the dike tip (Dieterich et al. 2000). In addition, the earthquakes used here are not relative relocations, so there is an associated error in their latitude/longitude/depth locations. Since the various models and seismicity do not indicate any significant lateral extension of the dike after 1 August, this continued motion to the SE, observed on FAFC, is likely associated with widening of the dike. A small motion to the SE continues in the days following the eruption onset, likely associated with the initial high influx rate.

To gain an improved understanding of the distribution of opening along the dike, a multipatch model was undertaken using the same input data as detailed above. This involved setting the location of the PB and dike based on the optimal model parameters derived from the uniform model. The dike and PB extended from the surface to a depth of 7 km, comprising a denser patch sampling in the upper 3 km. Two million iterations were run to ensure the model had fully converged. This model can reproduce the observable deformation with minimal residuals (Fig. 9).

Fig. 9
figure 9

Final geodetic model for dike 3. Panels show a input data: observed horizontal GNSS displacements (black arrows) and line-of-sight change in the Sentinel ascending interferogram (T16 20/07/2022-13/08/2022) (color), b predicted ground displacements and LOS change, c residuals, and d distribution of the dike volume change in the series of models produced in the inversion. Red star in a shows location of eruption site on 3 August 2022. In b, red line shows modeled dike (D), magenta line is plate boundary segment (PB), black lines are modeled faults (F)

Full size image

In order to undertake a comparison of the 2nd and 3rd diking events, similar models were undertaken for dike 2 that formed in December 2021. The first model included a simple dike (with uniform opening only) and a vertical dislocation allowing for strike-slip motion only, to model slip along the PB. The model incorporated two Sentinel-1 interferograms as the input (ascending T16 16/12/2021-28/12/202 and descending T155 20/12/2021-07/01/2022, Fig. S9) and total GNSS displacements from 20 to 29 December 2021. Following this, the multipatch model was run in a similar manner to that detailed above, allowing for only opening on the dike and slip on the PB (Fig. S10). A comparison of distributed opening along both dikes is displayed in Fig. 10. The 95% CI volume change calculated from these models for dike 2 is 12–25 million m3 and for dike 3, 10–19 million m3.

Fig. 10
figure 10

Median distributed opening along dikes 2 and 3. a Dike 2 (intruded between 21 and 27 December 2021 and b dike 3 (intruded between 30 July and 3 August 2022). Red star shows the location of the eruptive fissure that opened on 3 August 2022, following the emplacement of dike 3

Full size image

Discussion

Considering the long timescales of historic eruptive periods on the Reykjanes Peninsula (Sæmundsson et al. 2020), repeated magmatic intrusions are likely to continue intermittently—potentially for several hundred years—some of which may solidify in the upper crust while others will migrate to the surface and culminate in an eruption. During the Krafla fires (1975–1984), 45% of the emplaced dike intrusions resulted in an eruption (20 dike intrusions and 9 eruptions, Einarsson and Brandsdóttir (2021), Buck et al. (2006)), whereas at Afar, Ethiopia, 29% of dike intrusions formed in a rifting episode between 2005 and 2010 resulted in an eruption (14 dike intrusions and 4 eruptions, Wright et al. (2012)). In these rifting episodes, eruptions tended to occur more frequently toward the later stages of the rifting episodes, whereas the initial dike intrusion events did lead to minor eruptions (Ebinger et al. 2010). If activity continues on the Reykjanes Peninsula and this analogy holds, additional and more voluminous eruptions would be expected as a continuation of the 2021–2022 activity. The size of the lava fields produced in 2021 and 2022 is indeed small compared to the amount of lava erupted in the previous period of eruptive activity on the peninsula (older historical lava fields; see Fig. 1). When comparing the 2021–2022 activity on the Reykjanes Peninsula to other rifting episodes like at Krafla and in Afar, it is important to consider that there is a difference in tectonic setting. Highly oblique spreading occurs on the Reykjanes Peninsula (Fig. 2), with a large component of shearing, unlike in Krafla or Afar. As a result, the magmatic segments of the Reykjanes Peninsula are arranged en echelon and overlap. Furthermore, they do not host a well-developed central volcano. The Reykjanes Peninsula events are eventually analogs of events taking place in oceanic spreading contexts where shearing is a major component of deformation, as in incipient, magma-rich transform zones.

There is significant uncertainty in forecasting the likelihood of a future eruption, which is essential to mitigate risk and minimize the impact of volcanic hazards—especially when an intrusion is within close proximity to towns and important infrastructure. There are both notable similarities and differences between the three recent diking events at Fagradalsfjall—although only dikes 1 and 3 culminated in an eruption. All dikes were intruded in a similar region at Fagradalsfjall, close to where the central axis of the Reykjanes Peninsula oblique plate boundary passes through. The central axis of the plate boundary is where the rate of strain accumulation due to plate movements is highest during quiescent periods, as their gradients in surface displacement rates are the highest (Sigmundsson et al. 2022). Dike 1, emplaced during February–March 2021, was significantly longer (~9 km) and also comprised 2 segments (with strikes of N45°E and N23.5°E). The median top depth of the dike was at 1.5 km (from initial uniform model) and initial magma inflow rate was 30–35 m3/s—although this declined to <10 m3/s prior to eruption onset on 19 March 2021. Dike 2, emplaced during December 2021, was ~4 km long, striking N42°E. The median top depth was ~2 km, and initial magma inflow rate was 22 m3/s. Dike 3, emplaced during July to August 2022, was ~3-km-long striking at N35°E. The median top depth was 1 km and initial magma inflow rate was 49 m3/s (Table 1).

During the dike propagations, significant shearing was observed along segments up to 20 km in length, along the plate-boundary (Sigmundsson et al. 2022). To compare the amount of cumulative geodetic moment released, we multiplied the median amount of left-lateral strike-slip (computed from each patch within the PB from the distributed shear models) by the length and width of the PB patches. Our results indicate that the shearing that occurred during the 2nd and 3rd diking events was comparable; however, this was only about one-quarter of the amount that occurred during the first diking event in February–March 2021. The first dike intrusion in February–March 2021 lasted for a period of approximately 3 weeks, whereas the second lasted 7 days and the third dike intrusion lasted only 5 days. There was a clear decline in both seismicity and deformation prior to the onset of the first eruption (Sigmundsson et al. 2022), whereas only a decline in seismicity was observed prior to the second eruption. Continued movement was observed on several GNSS stations (FAFC, GONH, LISK, and STAN, Fig. S1) after the 3 August 2022. This may be a result of the higher magma inflow rate into dike 3 and continued inflation of this dike following the eruption onset. The main similarities between dikes 1 and 3 are the shallow depth of the intruded magma forecast from the initial models (< 1.5 km top depth) and the higher initial magma inflow rates when compared to dike 2. It should be noted though that these values are dependent on the time span of the data that was utilized in the geodetic models—the estimated inflow rate is an average for the time span of the geodetic data.

The decline in seismicity was used in forecasting the likelihood of an eruption during the third diking event, in conjuction with additional geodetic observations and modeling results. During this decline in seismicity (initially identified on 1 August 2022), horizontal motion at GNSS station FAFC changed direction and began moving SE indicating continued inflation of the dike, the modeled magma inflow rate was high, and top depth of the dike was shallow. A decay in seismicity would, however, also be observed if the dike stopped propagating without reaching the surface and triggering an eruption (e.g., towards the end of the 2nd diking event, which did not culminate in an eruption (Fig. 3)). Nevertheless, since declining seismicity may be indicative of an impending eruption, a high level of alert should be maintained for several days, even after most of the seismicity stops. The character of seismicity and additional datasets should be interpreted in parallel for improved understanding of the processes and potential outcomes. A decline in seismicity was observed prior to the 2021 Fagradalsfjall eruption (Sigmundsson et al. 2022), and prior to the Bárðarbunga-Holuhraun eruption in 2014. In the days following a 14-day long dike propagation within the Bárðarbunga volcanic system (prior to eruption onset), the rate of seismic moment realease was considerably lower than that observed during the peak of seismic activity associated with the dike propagation (Ágústsdóttir et al. 2019; Ruch et al. 2016; Sigmundsson et al. 2015).

The August 2022 eruptive fissure occurred where dikes 2 and 3 overlap (Fig. 6). The seismicity preceding both dike intrusions started at depths >8 km in the area NW of GNSS station FAFC (blue-filled circles on supplementary Fig. S11). On 21 December 2021, the seismicity began at ~18:30 UTC and remained near FAFC for about 4 h, gradually shallowing to 6 km over the next 5 h, then moved to NE for 2 h. The direction of dike propagation then changed, and the dike began propagating rapidly to the SW at a depth of ~6 km and migrated roughly 5 km SW along the strike of the dike 2 intrusion during the following 5 h but did not propagate further upwards. The seismicity accompanying dike 3, on the other hand, propagates 4 km upwards in the first 5 h and about 1.5 km laterally towards SW then NE over the first 15 h. The multipatch opening models displayed in Fig. 10 show that during the second dike intrusion, opening was distributed over a wide depth range, with the bulk of the opening at depths >6 km, but shallowing towards the NE end of the modeled dike. The initial deeper propagation (during the first 5 h) does not appear to significantly contribute to the deformation, although the dike possibly extends up to 1 km further to the NE at a deeper level (based on seismicity), which is incorporated within the 95% CI bounds of the modeled length. Dike 3 displays a much narrower band of opening which is concentrated at a much shallower level, between 1 and 3 km. It may be possible that on 30 July 2022, an additional pulse of melt from depth was injected into preferential migration pathways, developed within the northeastern segment of dike 2, which in turn triggered the upward propagation of dike 3.

Conclusions

During extended episodes of volcanic unrest characterized by repeated magmatic intrusions, it is important to forecast the likelihood of an eruption (for example, which dike intrusions will culminate in an eruption and which may not) and likely location of fissure openings, to minimize risk and potential exposure to hazards at the time of eruption onset. A warning of high probability of an eruption at Fagradalsfjall was issued on 2 August 2022 based on joint interpretation of ground deformation, geodetic modeling results, and seismicity patterns. Observations were interpreted in terms of a dike intrusion characterized by a high magma inflow rate (about 49 m3/s) to a shallow depth (dike top at about 1 km depth). This was associated with an observable decline in seismicity, interpreted as a marker of declining release of tectonic stress in the area. An eruption commenced at ~13:30 UTC on 3 August 2022, when a 375-m-long eruptive fissure opened up within 49–110 m of the projected surface location of the modeled dike. Interpretation of ground deformation patterns over a wider area assisted by interpretation of seismicity suggested no upflow of magma in other areas of the Reykjanes Peninsula during this time, but significant fault slip, surface effects, and fractures, as well as shearing across an ~14-km-long segment of the central axis of the oblique spreading plate boundary on the Reykjanes Peninsula. The dike intrusion that preceded the eruption was characterized by a high magma inflow rate to a shallow depth with an observable decline in seismicity detected within the 38 h preceding to eruption onset. These parameters may be key to determining which future dike intrusions on Reykjanes Peninsula, and eventually elsewhere, will culminate in an eruption. The joint interpretation of multiple datasets (e.g., seismicity and deformation) in conjuction with geodetic modeling results is key to forecasting further activity.