We obtained the event parameters, as well as travel times, from the SIL catalogue. This seismic event catalogue is produced by the Icelandic Meteorological Office using the South Icelandic Lowlands (SIL) network and (due to the unrest) some stations from the Reykjanet network63,64 as follows65:
Automated event detection and location.
Manual quality check and phase-pick adjustment for accepted events by analysts.
These picks are then used to obtain the final origin times and locations.
At the time there were 18 permanent seismic stations operated in the selected area, of which 5 came into operation in early 2020. We present these final locations directly using Voxler66,67.
Local earthquake tomography
The seismic traveltime tomography was done with the PStomo_eq algorithm16. From P- and S-wave traveltimes from local earthquakes the 3D P- and S-wave velocity structure is solved for simultaneously with new hypocentral parameters. The joint velocity-hypocentral problem is decoupled with the method of Pavlis and Booker68. Forward traveltimes are computed with finite difference solver time3D69,70. The conjugate gradients solver LSQR71 is used for solving for the velocity structure. The process is repeated over several iterations as the problem is non-linear. After each iteration new rays are computed in the updated velocity models from the updated earthquake locations.
Data selection for the tomography
The data were selected from the period 15 December 2019 to 27 February 2021, a few days into the swarm. A rectangular model area 90 km by 50 km was first selected based on the seismicity distribution and the location of the seismic stations. At the time there were 18 permanent seismic stations operating in the selected area, of which 5 came into operation in early 2020. For selection, events needed to have a minimum of seven picked phases (P or S), and a seismic station located within the distance of 1.5 times its focal depth to have good constraints on the hypocentral depths72. Some obvious outliers were also discarded. This resulted in the selection of almost 24 000 events (about 55% reduction). The reason we did not include later events for the tomography is that almost no further events occurred below 8 km depth (in the region of interest) and thus did not contribute to the imaging of the deep magma reservoir. Extending the selection further back in time did not contribute to the illumination of the area either, due to that few stations were operating in the area.
Model parameterisation
Seismicity is not evenly distributed within the model region. In the seismically active regions, the cut-off in seismicity is generally somewhere between 4 and 6 km depth. In an attempt to adjust the model discretisation accordingly, the model discretisation was adjusted with depth. Above 6 km depth we used rectangular model blocks 0.5 km thick and 1 km in horizontal directions. Below this depth the blocks were 1 km thick and 2 km in horizontal directions. Traveltimes were computed on a uniform grid of 0.25 km to maintain sufficient accuracy.
Model regularisationSmoothness constraints
To avoid ‘wild’ velocity variations, particularly in poorly resolved model regions, the model was regularised by applying smoothing constraints. This was implemented as soft constraints in which the Laplacian of the P- and S-wave fields was minimised along with the traveltime data73. The weight of the Laplacian function relative to the travel time data was controlled by a weighting parameter (c.f. Equation 3 in16). To find a suitable value for this parameter a series of tests were performed. An ’optimal’ value would be one that results in low data variance without the velocity models having unnecessary structure. Figure S2 shows the data variance versus the squared length of the velocity models for different values of the smoother. The model length was computed as the squared model deviation from the mean velocity at every layer. We selected a smoother with the value of 10 for both models, beyond which the data fit did not improve without simultaneously increasing the model length.
Cross-gradients constraints for structural similarity
The P- and S-wave velocity models are only coupled to each other in the earthquake locations. It is therefore common to employ a further constraint to limit the ratio between the models (Vp/Vs ratio damping). Here we applied a structural constraint between the P – and S-wave velocity models, by requiring that the cross product of the gradients in the respective models is everywhere zero. The resulting models are thus structurally similar74,75, without any constraints having been placed directly on the Vp/Vs ratios. It is a reasonable assumption that variation in state variables and composition affect both the P- and S-wave velocities, though it is well known that rigidity and compressibility are differently affected. The cross-gradient constraints do not force the P-and S-wave velocity models to be similarly affected (as would a Vp/Vs ratio damping). The assumption is rather that variations in the P-and S-wave velocity models occur at the same locations. A change in one model may thus be opposite to the change in another model. Synthetic tests show that this results in more realistic Vp/Vs ratios than if e.g. damping towards a pre-defined Vp/Vs ratio is implemented75.
Final model
The final velocity models are discussed in the main text. A Vp/Vs ratio high is shown between 9 and 12 km depth at about 20 km along the cross section beneath Fagradalsfjall (Fig. 3) that is interpreted to represent the upper part of a magma reservoir located in the lower crust close to the crust-upper mantle boundary and feeding the 2021–2023 Fagradalsfjall eruptions.
Model appraisalCheckerboard reconstructions
Demonstrating synthetic checkerboard reconstruction tests is a common way to show which parts of a seismic model is robust. In those, checkers of alternating velocity perturbations are superimposed on a starting model and synthetic travel times are computed. Sometimes noise is added to this data, which is then inverted. Regions in which the checkers are well reconstructed are then deemed reliable in the real models. This may provide a general overview of which parts of the models are well resolved. To make the test ‘difficult’ and demonstrate that the cross-gradients constraints do not damp the Vp/Vs ratios per se, the P- and S-wave velocity checkers have reversed polarities. Below 6 km depth the checkers are larger, a consequence of the inversion cells being larger below this depth. Figure S3 shows the result of Vp/Vs ratio of the ‘true’ (synthetic) and the reconstructed models in the right and left panels, respectively. The Vp/Vs ratio is well reconstructed between 1 and 7 km depth, in general where seismicity occurs. In the region of the high Vp/Vs ratio anomaly in the depth range of 9 to 12 km depth we note that the checker is reconstructed, though not to its full strength. The coefficient of correlation between the reconstructed and original checkers is computed and the outline of the 0.6 isoline is shown in Figs 3 and S3 to outline were the checkers are fairly well reconstructed.
Hypothesis tests
To make targeted tests for certain features in a model, a specific anomaly may be removed from the final model. A few more iterations are then performed with the real data in order to investigate if the anomaly is required by the data, i.e. if the removed features are ‘put back’ into the models by the data. In Fig. S4 we have replaced the high Vp/Vs anomaly below 9 km depth with ‘normal’ velocities (panels a and c), resulting in a worse data fit. The result after a final few iterations demonstrate that the anomaly is put back by the data, returning to the previous data fit (panels b and d).
Characteristic model test
Whereas a checkerboard test may provide an overview of which parts of a model is well resolved, other synthetic models may be investigated to answer specific questions about the model. Fig. S5 is a synthetic model showing a region of high Vp/Vs ratios, about 6 km across and extending to great depth (right panel). The high Vp/Vs body also features a core of even higher Vp/Vs ratios (2.04). In the synthetic model a low Vp/Vs ratio layer resides above the high Vp/Vs ratio region. The model has some resemblance with the final model. The reconstruction shows that the region of high Vp/Vs ratios is clearly observed. It is, however, broader than the true model, and the core is not recovered. The depth extent of the Vp/Vs anomaly is also not possible to determine. The reconstruction also shows that the low Vp/Vs ratio layer is well reconstructed in the center of the model, but smeared out towards the south west in the model. Fig. S6 shows a much broader high Vp/Vs ratio anomaly (14 km across). In the reconstruction this model is also smeared out, but less so than the smaller anomaly in Fig. S5.
Lessons learned from the model tests
The hypothesis test shows that high Vp/Vs ratio anomaly beneath Fagradalsfjall at depths between 9 and 12 km (Fig. 3) is required by the data. When the anomaly is replaced by ‘normal’ velocities in the hypothesis test, the inversion puts it back into the model (Fig. S4). We also tested many different starting models (not shown), and the feature is persistent in the models regardless of starting models. The checkerboard reconstruction test (Fig. S3) shows that in the centre of the model the checkers are well reconstructed down to 8 km depth. Below this depth, the checkers are smeared out. Two characteristic model tests (Figs. S5 and S6) show that the dimensions of the high Vp/Vs ratio anomaly are difficult to determine. A 6 km wide anomaly appears wider than it is. Testing for a much wider feature (14 km) shows that even though it is smeared out, the anomaly would certainly be wider than in our results. From this we infer that the top part of the magma reservoir is likely 10 km or smaller in the horizontal direction. These tests indicate that the shallowest part of the reservoir at 9 km depth is robustly reconstructed, but also that it is impossible for us to determine its depth extent.
Relocated seismicity
The entire seismic catalogue between 24 February and 15 March 2021, and similarly for the 2022 and 2023 eruptions, was relocated in the final velocity models. The location errors (1 SD) are shown in Fig. S7. The results indicate that the events are most accurately located in the east-west direction, probably due to the east-west elongation of the seismic network. The depth is also the parameter with the largest error. Some 90% of the relocated events have smaller depth uncertainty than 0.6 km, and 60% of the events have an uncertainty in depth of 0.3 km or less.
Further support for the high Vp/Vs anomaly indicating a magma reservoir.
Though Vp/Vs interpretation is non-unique, as is common in geophysics, the following observations support our interpretation as to the high Vp/Vs anomaly coinciding with the top of the source magma reservoir at about 9 km depth below the eventual eruption site. (1) The location of the anomaly below the eruption site. (2) The events, particularly the earthquake swarms and their migration, leading up to the eruption. (3) The eruption itself and particularly the location of the volcanic fissure. (4) The long-period seismicity coinciding with the anomaly. (5) The geochemical and petrological results indicating that the primitive magma came from a deep-seated reservoir and not from a shallow magma chamber11,18,19,29,36.
That the high Vp/Vs would be due to pressurised liquid water76 is not tenable for the present anomaly, partly because of the great depth (and thus pressure) and high temperatures – far above 580 °C (the brittle-ductile transition beneath the Reykjanes Peninsula16 occurs at temperatures between 580 °C and 750 °C). Any available volatiles and/or meteoric fluids (unlikely to be found at these depths), must be in a compressible state, which would produce not higher but lower Vp/Vs ratios. The propagation of seismicity, which starts at the top of the high Vp/Vs anomaly, is also a very specific signal favouring our interpretation of the anomaly as the top of a magma reservoir.
The location of the anomaly also coincides with where one would expect the source reservoir of the eruption. This follows, first, because the path of a regional dike, such as the present one, is normally close to vertical5,57. Second, as indicated by Eq. (1) and discussed earlier, a reservoir would be expected to be directly below the boundary zone (Fig. 1b), where the crust is comparatively thin and highly fractured (Fig. 2), namely at the site of a local potential energy minimum. Third, the roof of the reservoir should be located where the dike-propagation induced seismicity begins, as is the case here. Fourth, as the source region deepens with time, we would expect the seismicity to deepen as well (as is seen in syn-eruptive seismicity). The deep long-period (DLP) seismicity at the location of the reservoir77 could also be expected to occur inside a reservoir preceding and during an eruption. The reservoir is a poroelastic body and thus only partially molten – most of the reservoir is solid5. Several mechanisms have been proposed for earthquakes inside magma reservoirs58,78,79,80, but the most likely one is increase in pore-fluid pressure due to melt/magma accumulation. Pore-fluid pressure increase triggers earthquakes on existing fractures, as well as on new fractures, as is well known from rock-physics experiments for developing enhanced geothermal (EGS) and hydrocarbon reservoirs58,67,81,82. Melt migration and pore-fluid pressure changes follow from the anomaly being the top of the source reservoir, the roof rupture, and dike-segment injections. Thus, in our interpretation, DLP events at the observed depth (10–12 km) are exactly as expected.
Most of the petrological and geochemical results of Halldorson et al.18 (their Fig. 4) and Troll et al.19 agree with the top of the reservoir being at a depth of about 9 km. The fountain tephra results, which Halldorson et al.18 base most of their interpretation on, suggest a greater depth for the origin of the magma, but that magma was not issued from the fissure until 38 days into the eruption (see their Fig. 3a). By that time, some 20 million m3 of material had already erupted13. The source depth of 15 km depth for magma erupted after 38 days of eruption11 is in perfect harmony with our reservoir model. Although there were not enough seismic signals to image the deeper parts of the reservoir, volcanotectonic and petrological considerations indicate that the reservoirs beneath the Reykjanes Peninsula, including that of Fagradalsfjall, are many kilometres deep1,5,19 (Fig. 2). As the eruption progresses, the magma is drawn from deeper parts of the reservoir, in accordance with Eq. (1) below. More specifically, Halldorsson et al.18 (see also Bindeman et al.11) show that the magma source deepened during the 2021 eruption: the magma erupted early in the eruption was more depleted and came from a source at lower pressure than later in the eruption, when the magma was more enriched and from a greater pressure. Thus a stratified reservoir is a likely source – as has been suggested for reservoirs in Iceland in general and for those beneath the Reykjanes Peninsula in particular1,5. In such a reservoir, the depleted magma in the top part of the reservoir would be the first to erupt, whereas the enriched would erupt as the eruption progressed – exactly as observed in the present eruption. We conclude that, although the resolution at depth is not good for small features, the top part of the reservoir is properly resolved.
Reservoir modellingDarcy’s law
Darcy’s law for flow in porous media may be expressed as follows5,17:
$$\vec{q} = \frac{Q}{A} = – \frac{{k\rho _{m} g}}{{\mu _{m} }}\nabla \emptyset$$
(1)
Here the vector \(\vec{q}\) is the discharge or Darcy velocity and denotes the volumetric flow rate Q per unit area A normal to the direction of flow, k denotes the intrinsic permeability, ρm the density of the magma or melt, g the acceleration due to gravity, and μm the dynamic (absolute) viscosity of the magma/melt. The symbol \(\emptyset\) denotes the potential energy or total head of the magma/melt, and ∇∅ is thus the gradient of the potential energy. The negative sign is to indicate that magma/melt migration is driven from regions of higher potential energy (head) to regions of lower potential energy. Thus, magma/melt in a partially molten upper mantle is driven to regions of local minimum potential energy. While part of the magma/melt flow in the reservoir may be through fractures, in which case the cubic law would apply, the fractures are likely to be mostly of grain-size dimensions and their fluid transport may be approximated through equivalent porous flow5,18.
Conditions for reservoir rupture
The conditions for magma-reservoir rupture is as follows5:
$$\:{p}_{t}={p}_{l}+{p}_{e}={\sigma\:}_{3}+{T}_{0}$$
(2)
Here \(p_{t} = p_{l} + p_{e} ~\) denotes the total magmatic pressure in the reservoir, pl the lithostatic stress at the site of the rupture of the roof of the reservoir, pe the reservoir excess magmatic pressure, that is, the pressure in excess of σ3, the minimum compressive (maximum tensile) principal stress, and T0 the local in situ tensile strength of the roof at the rupture site. Eq. (2) implies that when the total fluid pressure in the reservoir reaches the minimum principal compressive stress plus the in-situ tensile strength, the chamber roof ruptures and a magma-filled fracture, here a dike-segment, is injected.
Magmatic overpressure and dike dimensionsMagmatic overpressure
For a magma-filled fracture, the overpressure po (assumed constant at a given depth in the fracture) is obtained from the following Eqs. 5 ref24:
$$p_{0} = \frac{{\Delta u_{I} E}}{{2L\left( {1 – v^{2} } \right)}}$$
(3)
where ΔuI is the opening of the dike fracture, E is Young’s modulus and ν is Poisson’s ratio of the host rock, and L is the strike-dimension or length of the dike fracture.
When a dike-segment begins to propagate up into the reservoir roof towards the surface, the magmatic overpressure po in the dike-segment (magma-filled fracture) can be estimated thus5
$$\:{p}_{0}={p}_{e}+\left({\rho\:}_{r}-{\rho\:}_{m}\right)gh+{\sigma\:}_{d}$$
(4)
where pe is the excess magmatic pressure in the reservoir at the time of roof rupture (equal to the in-situ tensile strength of the roof rock), ρr is the average host-rock density, ρm is the average fluid (magma) density, g is the acceleration due to gravity, h is the dip-dimension of the dike (the height of the dike-segment above its source), and σd is the differential stress (the difference between the maximum and the minimum principal stress) in the host rock at the location where the dike overpressure is of interest. In Eq. (4) the term \(\:\left({\rho\:}_{r}-{\rho\:}_{m}\right)gh\) represents the buoyancy.
We focus on the dike thickness at a crustal depth of 2 km where much of each of the dike-injections of 2021, 2022, and 2023 became arrested (Figs. 5, 6 and 8, S8 and S9). Here we consider the 2021 injection, when much of the multiple dike was emplaced. If we use typical generalised crustal density values for the units that constitute the crust of Iceland28, the reservoir depth of 9 km, the excess pressure/ in-situ tensile strength of the roof as 3 MPa5, the above magma density of 2800 kg m− 3, and σd as 1 MPa5, Eq. (4) gives the dike overpressure at 2 km depth as about 7 MPa. To obtain the dike thickness, we rewrite Eq. (3) and solve for the dike opening or thickness thus:
$$\:{\Delta\:}{u}_{I}=\frac{2{p}_{0}\left(1-{v}^{2}\right)L}{E}\:\:\:$$
(5)
which yields about 4 m thickness at 2 km depth.
Comparison with other dike-dimension and dike-overpressure estimates
Other estimates of the dimensions of the dike are generally similar to our results. Thus, the strike-dimension is estimated at about 9 km32,48. A dip-dimension of 8 km has also been estimated, but in some cases as the depth from the volcanic fissure at the surface to the base of the dike48. In our estimate, however, the 7 km dip-dimension is that of the main (arrested part of the) dike, from a crustal depth of 2 km to the depth of 9 km (the top of the magma reservoir). The maximum thickness (dike-fracture opening) estimate using the model of elastic dislocations in an elastic half-space48 (that is, a homogeneous and isotropic crustal segment, so with no layering) yields about 3 m. This is not greatly different from our result of 4 m which, though, is based on an entirely different approach using layering and principles from fluid mechanics and fracture mechanics.
Elastic-dislocation models are not suitable for explaining dike-propagation paths, however. For these, appropriate principles of fluid and solid (including fracture) mechanics are needed. The mechanical details as to how dikes select their propagation paths are given elsewhere5,57, but in the present paper we infer the paths of the main dike-segments from 2021, 2022, and 2023, from rupture and segment injection in the roof of the reservoir, to arrest of the main part of the dike, and finally to the eventual eruption and volcanic fissure formation of a tiny part of the main dike (Fig. 5).
One main parameter for dike propagation is the magmatic driving pressure or overpressure. This is the pressure that ruptures the crust, drives the dike propagation, and contributes largely to the eventual dike dimensions. One model of the 2021 eruption has a 13 km tall vertical pipe, with a lateral cross-sectional area of a few square metres, transporting magma from a depth of 19 km in the mantle to the base of the dike48. We know of no plausible mechanisms by which a pipe of such dimensions could be generated, particularly given the short time window for its formation. In this model a constant density difference of 300 kg m− 3 is assumed between the host rock and the magma up to the base of the dike48. For a 13 km tall pipe, the overpressure would then be at least 38 MPa at the base of the dike. Since the model also assumes no density difference between the rock and magma in the dike itself, the 38 MPa overpressure should be largely maintained inside the dike.
With an overpressure of 38 MPa then, from Eq. (5) and using the same elastic constants as above, the thickness of the dike would be about 14 m for a dike dip-dimension of 6 km and about 19 m for a dip-dimension of 8 km – thus including the feeder-dike part. (The authors use a dike dip-dimension of 8 km in part of their paper and 6 km in other parts.) These openings or dike thicknesses are 5–6 times larger than the maximum thickness of 3 m that the authors estimate using their dislocation model. If the same dike overpressure would be assumed to extend to the surface (8 km dip-dimension), then, from Eq. (5), the opening of the 180-m-long volcanic fissure at the beginning of the eruption would have been about 2.6 m, or more than 12-times the estimated opening24.
Using Eq. (4) for a layered crust, we estimate the maximum magmatic overpressure in the dike at 2 km depth as about 7 MPa. Using this overpressure, the maximum dike thickness at 2 km depth becomes about 4 m. Similarly, from Eq. (5) and the length/aperture ratio of the volcanic fissure, we obtain the overpressure of the feeder-dike at the surface as about 3 MPa.
The results also indicate that all the dike-injections during the events of 2021, 2022, and 2023 used largely the same paths. More specifically, the injections of 2022 and 2023 used parts of the already established dike from the main 2021 injection as their paths, resulting in a multiple dike. Such dikes are very common and are thought to make eruptions mechanically easier and thus more likely83.