High temperature (HT) geothermal fields are associated worldwide with volcanism and fluid convection at plate boundaries and hot spots. Recent volcanic and seismic unrest at the Reykjanes Peninsula (RP) plate boundary in SW Iceland and the Svartsengi HT field (Fig. 1) reveals a previously unexplored cyclic interaction between tectonic spreading, magmatic reservoirs, and the detachment of supercritical magmatic fluids interacting with deep hydrothermal waters.
On January 22nd, 2020, an earthquake swarm started 3 km east of Svartsengi. Simultaneously, continuous GNSS stations showed uplift at the geothermal field close to the re-injection site1 of the power plant, followed by subsidence. Three cycles of such inflations and deflations occurred until July 2020, always followed by continuous deflation and diminishing seismicity. A similar inflation started in August 2020 in the centre of the Krýsuvík HT field, 20 km east of Svartsengi (Fig. 1). In February 2021, crustal extension and an intense earthquake swarm revealed the formation of a NNE striking magmatic dyke in between the two HT fields, followed by a fissure eruption about 9 km east of Svartsengi.
The spreading axis of the Mid-Atlantic Ridge comes on land at the SW corner of the RP. There it bends into a 60 km long N70°E striking oblique plate boundary expressed by a 5-10 km wide seismic and volcanic zone, with large episodic earthquake swarms every 20-40 years2 with magnitudes up to M 6, mostly on N-S trending strike-slip faults. Volcanic episodes occurred at intervals of 800-1000 years during the past 4000 years, the last one ended in 1240 AD3. Each volcanic episode might last for 1-3 centuries with occasional basaltic lava flows from N45°E trending fissures extending into the adjacent plates4. HT geothermal fields with reservoir temperature of 240-330°C at 1-3 km depth5 have formed at the intersection of the seismic zone and the main volcanic fissure swarms (Fig.1). One well, IDDP-2 at the Reykjanes HT field, was drilled to 4.6 km depth6 where bottom hole temperature is estimated to be about 600°C7 just beneath an underpressurized aquifer (<35 MPa) revealed by circulation loss in the water filled well during drilling.
The upper crust of the RP is close to 4.5 km thick8,9,10 and composed of basaltic extrusives mixed with a downward increasing occurrence of intrusives. The lower crust down to Moho at ~15 km depth is thought to be made of intrusives with no evidence of melt 10,11,12. The brittle-ductile boundary (BDB) below the seismic zone on the RP is generally at 6-7 km depth, doming up to 4-5 km depth below the HT fields13,14,46. The estimated temperature at the BDB is 600-700°C 1,15)
The 76 MWe power plant at the Svartsengi HT field and the Blue Lagoon Spa (Fig. S1) are the heart of the Geothermal Resource Park, providing electricity and hot and cold water to 25,000 residents and local industries16. The annual production is about 450 kg/s, of which about 300 kg/s are re-injected into the reservoir.
The relation of cyclic deformation and earthquake activity at two distinct HT fields at the RP plate boundary and the time lag until a distant fissure eruption broke out is unusual and previously not understood. Using a comprehensive modelling approach, we show for the first time that the detachment and ascent of magma-derived volatiles into a sealed aquifer at the BDB beneath Svartsengi HT field can explain the strong uplift and subsidence cycles and induced seismicity and can be interpreted as precursor to a coming eruption through the instability of a deep magma reservoir in the lower crust or upper mantle.
Transient deformation in the geothermal field
Analysis of twelve months of satellite InSAR time-series data (see methods) processed in both ascending and descending configuration reveals an elliptical area exceeding 80 km² affected by episodic uplift and subsidence, along with minor horizontal displacements (Fig. 2a). The elliptical area has a major axis of ~12 km striking N60°E and minor axis of 8.5 km. The major axis follows the orientation and main strike of the geothermal reservoir (Fig. S1) but deviates both from the N45°E strike of volcanic fissures and the N to NNE strike-slip faults mapped at the surface on the RP.
Three events of sharp uplift episodes were followed by gradual subsidence. The duration of the three successive uplift episodes increased while the uplift rate decreased correspondingly (Fig. S2). The first episode had the fastest uplift rate of 2.2 mm/day and climaxed at a maximum uplift of 66 mm after 30 days. It was followed by 18 days of subsidence, a total of less than 10 mm. The second episode had an uplift rate of 1.1 mm/day and produced an uplift of 55 mm in 48 days. It was followed by a faster subsidence episode lasting approximately 24 days with a subsidence of 16 mm. The final episode with an uplift rate of 0.5 mm/day produced a cumulative uplift of 32 mm within 60 days. On July 18th, an earthquake of magnitude M4.1 occurred near the inflation centre17. It was followed by a subsidence period visible on SAR until mid-December with an average subsidence rate of 0.2 mm/day and total subsidence of 38 mm. Overall, the cumulative total uplift observed at the centre of displacement was close to 150 mm and by December 2020 the actual uplift was reduced to 90 mm.
Modelling of inflation-deflation cycles
The deformation pattern suggests volume increase at depth. We first model the uplift with a simple point source in an elastic half-space18 to retrieve information about the location, depth and strength of the source. Inversion modelling suggests a source depth of 4.4 km. The horizontal position is stationary and varies only around 330 m between the three uplift episodes. The inferred volume change for the three uplift episodes is 4.70, 3.55 and 2.75 ·106 m³, respectively, and cumulatively explains over 90% of the observed deformation. (Fig. S3).
To better match the elliptical shape of the uplift pattern, we also used rectangular dislocation models19. The uplift episodes can be modelled by three 10 m thick, 7-9 km long and 30-50 m narrow nearly horizontal intrusions at a depth between 3.7 – 4.4 km with an average strike of N60°E.
To explain the uplift and subsidence cycles, the time-dependent gravity anomalies and seismicity in the region, we tested a poroelastic model considering a strongly coupled diffusion and deformation process. This model suggests a thin, permeable porous aquifer layer in the depth of 4 km, embedded in a multi-layered poroelastic halfspace. The elastic parameters are based on the velocity model used for earthquake locations (Tab. S3). During the three uplift episodes, the aquifer was pressurized through fluid intrusion along a N60°E trending line source where the magnitude of pressurisation was attenuated from the centre along the line by a Gauss-taper with a standard deviation of 2 km. We inverted for the inflow history over the full 12-month period of unrest, with an inflow of fluids during periods of inflation and zero inflow during periods of deflation (Fig. 2e). The resulting constant intrusion rates are 13.4, 8.1 and 5.1 m3/s, respectively, causing a total intrusion volume of about 9.5 ·107 m3. The induced deformation and stress changes related to the inflow and pore-pressure diffusion can explain the observed deformation cycles, particularly the fast subsidence following each cyclic uplift, but also the observed free-air gravity anomalies and seismicity patterns in time and space.
Implications from free-air gravity anomalies
Gravity combined with deformation data provides additional insights, as changes in gravity depend on the mass of the intruded fluids, while deformation depends on volume change. We measured gravity in four consecutive campaigns (Tab. S1) at 10-12 existing permanent stations along an L-shaped profile extending north and west from the centre of uplift (Fig. S6).
Gravity changes between the first and second campaign show a consistent free-air corrected gravity increase of 10-14 μGal around the centre of uplift (Fig. S7, Tab. S1). This value is used in the poroelastic modelling to distinguish the type of intruded fluids. A consistent free-air corrected gravity decrease of 4-8 μGal is observed at the uplift centre from April to October. Between October 2020 and February 2021, the free-air corrected gravity decrease continues at most stations. The maximum decrease averaged over the three closest stations at the centre was 5 μGal.
The poroelastic model explains the free air gravity changes (Fig. 2b & Fig. S7) provided that the density of the injected fluid at 4 km depth is close to 840 kg/m3 or smaller if fluid compressibility is large. This disqualifies magma as a direct cause of the pressure pulses at 4 km depth, and is supported by the fast subsidence, and points to the inflow of low-density material like water and/or CO2.
The observed seismicity provides constraints on the unrest processes at depth. In general, earthquake rate increases where stresses are elevated, and earthquake locations constrain the geometry of the deformation source. Migrating seismicity fronts provide evidence of propagating intrusions or migration of pore pressure in aquifer systems. We followed two approaches to incorporate seismicity: 1) analysing the seismicity rate changes relative to the background activity, and 2) compiling a highly accurate earthquake catalogue during the uplift episodes.
To evaluate the background seismicity, we used the national Iceland Met Office (IMO) catalogue20. Specifically, we estimated the rate changes during the deformation changes in space and time. For an estimated completeness magnitude, Mc, of 1.5, we first determined the average background rate between 2000 and 2020 on a 3D grid surrounding the deformation centre. The resulting rate changes were finally compared to model predictions based on Coulomb failure stress (CFS) changes induced by the deformation source derived from the InSAR data. Our model (Fig.3) explains the spatial and temporal distribution of the earthquakes as seismicity induced by fluid pressure increase in the roots of the Svartsengi HT field. This finding provides independent evidence for the geodetically constrained poroelastic aquifer model.
For a more detailed analysis of the spatiotemporal patterns during the unrest period, we created a new catalogue for 2020, using 26 seismic stations spaced up to 30 km around Svartsengi and, for the first time, integrated distributed acoustic sensing (DAS) data, from a 21 km long fibre optic telecommunication cable buried 80-90 cm below the surface from the southern tip of Reykjanes to Grindavík crossing the Svartsengi HT field21 (Fig. S8). The new catalogue covers the period from February 1st to August 30th. We developed a new waveform stacking and migration method applied to continuous data streams for the detection and localisation of the smallest earthquakes. We detected 39,500 earthquakes with magnitudes M >-1 and localise the majority automatically. The locations have high quality, both laterally and vertically, since the sensors and DAS cable were located directly above the uplift zone, and the azimuthal distribution of all stations was extraordinarily good.
The local seismicity during the Svartsengi uplift is very shallow. The lack of deeper earthquakes below the uplift centre supports the updoming of the BDB to a depth of 4 km (Fig. 4b). Interestingly, the long and narrow axis of updoming of the crustal seismicity resembles well the elliptical uplift pattern (Fig. S9). The new catalogue also shows that the very shallow seismicity mostly occurs at the centre of uplift and dips to lower levels outside the uplift-affected region, indicating that these earthquakes are probably induced by stresses generated by the bending of the crustal layer above the deformation source.
A new and wholistic model of the pre-volcanic processes
Our poroelastic model fits the deformation, gravity and seismic data (Figs. 2 & 3). It consists of a pressure and volume increase due to fluid intrusion into an aquifer at 4 km depth that diffuses with time and has density up to 840 kg/m3. These results disqualify magma as the intrusion material. Although re-injection fluid might comply with the density value, the re-injection rates are far too small to explain the volume and gravity change. Therefore, we propose the intrusion of magma-derived supercritical fluids, mainly CO2. The temperature in the aquifer at the bottom of the Svartsengi HT field is likely to be close to 330°C prior to the intrusion but not higher than the 600°C at the BDB. To get a maximal density of 840 kg/m3 for CO2 as intruded mass at 330°C, an absolute pressure of 180 MPa is needed22. This value exceeds the estimated 100-110 MPa lithostatic pressure at 4 km depth and can explain the uplift. Such a high overpressure may only develop if the fluid is temporarily connected to parts of the deeper crust. Considering the fluid compressibility might reduce the density estimate and the required pressure23,24. Contribution of sulphuric gasses will have a similar effect. An elevated concentration of volatiles in the lower crust beneath the Svartsengi HT field is also independently confirmed by high Vp/Vs ratios in seismic tomography25, 26.
Our earthquake analysis shows that the BDB domes up from 6-7 km depth to 4 km depth beneath the centre of inflation, where the deformation data suggest the volume increase to occur. Measurements in the deep exploration wells IDDP-1 at the Krafla geothermal field and IDDP-2 at the Reykjanes geothermal field6 show that underpressurized aquifers exist near the BDB. The crust is in both cases fully elastic and the fluid pressure in the rock is hydrostatic down to the BDB.
We propose a model that fits available observations of gravity, seismicity and surface deformation, including the uplift, subsidence, occurrence and time scales (Fig. 4 & Fig. S11). Decompression melting of slowly rising magma occurs below the Moho27 which is at 15 km depth10. Volatiles exsolved from the melting process, migrate upwards and become trapped at the BDB at ~7 km depth, generating strong overpressure, but not high enough to lift the overburden (~220 MPa). This might have triggered the earthquakes swarms late 2019. After reaching a certain limiting volume, magmatic volatiles flew upward along the BDB toward the underpressurized aquifer (< 35 MPa) at 4 km depth at the bottom of the convective HT field where they intrude the aquifer with pressure high enough to cause the uplift (>110 MPa).
The feeder channel model for the lower crust explains episodic inflow with exponentially decreasing intensity and increasing recurrence times. The exponential decay is expected for pressure-driven draining from a deep reservoir with finite volume28. Cyclic transport of fluid batches through the channel is expected when a valve mechanism is present and injects accumulated CO2 volume into the channel once the valve is opened. The first three cases the fluid batches were directed to Svartsengi HT field but the fourth one in August was directed to the Krýsuvík HT field where conditions as in Svartsengi are expected (Figs. S10 & S11).
Our model is consistent with the eruption that started March 19th, 2021. In contrast to violent initial phase of eruptions derived from a magma chamber, the eruption rate was gentle during the first month when it substantially increased. In our model the volatiles from the degassing magma had already detached from the initial magma batch prior to the eruption rising, from the upper mantle and trapped at the BDB where it forced a different path through the crust to the underpressurized geothermal reservoirs. From the total amount of CO2 injected into the roots of the HT fields, we can estimate that the minimum volume of degassed magma beneath the present eruption site is of the order of 10 km3 (information sheet IS2). Consequently, the available volume of magma is neither a limiting factor for the longevity of the eruption nor the erupted volume.
In conclusion, our model explains for the first time the role and importance of high temperature geothermal fields in the complicated mechanism leading to volcanic eruptions at oblique divergent plate boundaries of the oceanic crust.