One year of cyclic unrest in a hydrothermal eld as a harbinger of a volcanic eruption

One year of deformation and seismicity prior to a volcanic eruption in March 2021 at an oblique plate boundary in Iceland created a unique opportunity to study the interaction between upwelling magma and geothermal processes. We apply poroelastic modelling to explain satellite geodetic data showing three uplift and subsidence cycles at the Svartsengi geothermal eld and use gravity data to constrain the density of intruded material. We use recordings on optical cable to generate a high-resolution earthquake catalogue and developed new waveform stacking and migration methods to detect and locate 39,500 earthquakes. The resulting model explains the geodetic, gravity, and seismic data by magmatic derived gas intruded into a horizontal sealed aquifer at 4 km depth in the roots of the geothermal eld at the top of up-doming brittle-ductile boundary. The total injected volume is estimated 9.5·107m3 with optimal density of 840 kg/m3. Our results suggest upward migration of three packages of volcanic gas along the brittle-ductile boundary from a subcrustal magmatic source 8–10 km east of the geothermal eld, with important implications for the dynamics leading to the eruption.


Background
High temperature (HT) geothermal elds are associated worldwide with volcanism and uid 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 eld (Fig. 1) reveals a previously unexplored cyclic interaction between tectonic spreading, magmatic reservoirs, and the detachment of supercritical magmatic uids interacting with deep hydrothermal waters.
On January 22 nd , 2020, an earthquake swarm started 3 km east of Svartsengi. Simultaneously, continuous GNSS stations showed uplift at the geothermal eld close to the re-injection site 1 of the power plant, followed by subsidence. Three cycles of such in ations and de ations occurred until July 2020, always followed by continuous de ation and diminishing seismicity. A similar in ation started in August 2020 in the centre of the Krýsuvík HT eld, 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 elds, followed by a ssure 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 years 2 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 AD 3 . Each volcanic episode might last for 1-3 centuries with occasional basaltic lava ows from N45°E trending ssures extending into the adjacent plates 4 . HT geothermal elds with reservoir temperature of 240-330°C at 1-3 km depth 5 have formed at the intersection of the seismic zone and the main volcanic ssure swarms (Fig.1). One well, IDDP-2 at the Reykjanes HT eld, was drilled to 4.6 km depth 6 where bottom hole temperature is estimated to be about 600°C 7 just beneath an underpressurized aquifer (<35 MPa) revealed by circulation loss in the water lled well during drilling.
The upper crust of the RP is close to 4.5 km thick 8, 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 elds 13, 14,46 . The estimated temperature at the BDB is 600-700°C 1,15) The 76 MW e power plant at the Svartsengi HT eld 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 industries 16 . 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 elds at the RP plate boundary and the time lag until a distant ssure eruption broke out is unusual and previously not understood. Using a comprehensive modelling approach, we show for the rst time that the detachment and ascent of magma-derived volatiles into a sealed aquifer at the BDB beneath Svartsengi HT eld 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 eld
Analysis of twelve months of satellite InSAR time-series data (see methods) processed in both ascending and descending con guration 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 ssures 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 rst 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 nal episode with an uplift rate of 0.5 mm/day produced a cumulative uplift of 32 mm within 60 days. On July 18 th , an earthquake of magnitude M4.1 occurred near the in ation centre 17 . 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 in ation-de ation cycles
The deformation pattern suggests volume increase at depth. We rst model the uplift with a simple point source in an elastic half-space 18 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 ·10 6 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 models 19 . 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 uid 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 in ow history over the full 12-month period of unrest, with an in ow of uids during periods of in ation and zero in ow during periods of de ation (Fig. 2e). The resulting constant intrusion rates are 13.4, 8.1 and 5.1 m 3 /s, respectively, causing a total intrusion volume of about 9.5 ·10 7 m 3 . The induced deformation and stress changes related to the in ow 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 uids, 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 pro le extending north and west from the centre of uplift (Fig. S6).
Gravity changes between the rst 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 uids. 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 uid at 4 km depth is close to 840 kg/m 3 or smaller if uid compressibility is large. This disquali es magma as a direct cause of the pressure pulses at 4 km depth, and is supported by the fast subsidence, and points to the in ow of lowdensity material like water and/or CO 2 .

Seismicity
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 O ce (IMO) catalogue 20 . Speci cally, we estimated the rate changes during the deformation changes in space and time. For an estimated completeness magnitude, M c , of 1.5, we rst determined the average background rate between 2000 and 2020 on a 3D grid surrounding the deformation centre. The resulting rate changes were nally 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 uid pressure increase in the roots of the Svartsengi HT eld. This nding 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 rst time, integrated distributed acoustic sensing (DAS) data, from a 21 km long bre optic telecommunication cable buried 80-90 cm below the surface from the southern tip of Reykjanes to Grindavík crossing the Svartsengi HT eld 21 (Fig. S8). The new catalogue covers the period from February 1 st to August 30 th . 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 ts the deformation, gravity and seismic data (Figs. 2 & 3). It consists of a pressure and volume increase due to uid intrusion into an aquifer at 4 km depth that diffuses with time and has density up to 840 kg/m 3 .
These results disqualify magma as the intrusion material. Although re-injection uid 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 uids, mainly CO 2 . The temperature in the aquifer at the bottom of the Svartsengi HT eld 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/m 3 for CO 2 as intruded mass at 330°C, an absolute pressure of 180 MPa is needed 22 . 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 uid is temporarily connected to parts of the deeper crust. Considering the uid compressibility might reduce the density estimate and the required pressure 23,24 . Contribution of sulphuric gasses will have a similar effect. An elevated concentration of volatiles in the lower crust beneath the Svartsengi HT eld is also independently con rmed by high V p /V s ratios in seismic tomography 25, 26 .
Our earthquake analysis shows that the BDB domes up from 6-7 km depth to 4 km depth beneath the centre of in ation, where the deformation data suggest the volume increase to occur. Measurements in the deep exploration wells IDDP-1 at the Kra a geothermal eld and IDDP-2 at the Reykjanes geothermal eld 6 show that underpressurized aquifers exist near the BDB. The crust is in both cases fully elastic and the uid pressure in the rock is hydrostatic down to the BDB. We propose a model that ts 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 Moho 27 which is at 15 km depth 10 . 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 ew upward along the BDB toward the underpressurized aquifer (< 35 MPa) at 4 km depth at the bottom of the convective HT eld 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 in ow with exponentially decreasing intensity and increasing recurrence times. The exponential decay is expected for pressure-driven draining from a deep reservoir with nite volume 28 . Cyclic transport of uid batches through the channel is expected when a valve mechanism is present and injects accumulated CO 2 volume into the channel once the valve is opened. The rst three cases the uid batches were directed to Svartsengi HT eld but the fourth one in August was directed to the Krýsuvík HT eld where conditions as in Svartsengi are expected (Figs. S10 & S11).
Our model is consistent with the eruption that started March 19 th , 2021. In contrast to violent initial phase of eruptions derived from a magma chamber, the eruption rate was gentle during the rst 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 CO 2 injected into the roots of the HT elds, we can estimate that the minimum volume of degassed magma beneath the present eruption site is of the order of 10 km 3 (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 rst time the role and importance of high temperature geothermal elds in the complicated mechanism leading to volcanic eruptions at oblique divergent plate boundaries of the oceanic crust.

InSAR data processing
We exploited the Copernicus Sentinel-1A and 1B satellite SAR data, which are available over Iceland in ascending (08/01/ -15/12/2020) and descending (06/01 -13/12/2020) orbits with a revisit time of 6 days for each track. By calculating a single interferogram for the cumulative period between 20/01/-17/07/2020, we found four fringes in both geometries, representing a shift of approximately 11 cm in the line of sight (LOS) direction (Fig. S4). A more detailed InSAR time series analysis shows three in ation and de ation episodes (Fig. S2): rst uplift 19/01 -18/02/2020, second uplift 07/03 -24/04/2020, and third uplift 18/05 -17/07/2020. Due to the 6-day revisit time of Sentinel1A/1B acquisitions the exact days of in ation/de ation episodes may be shifted by a few days. The time-series analysis was carried out using the Small Baseline Subset (SBAS) algorithm 29 , implemented in SARscape, which is based on the combination of interferograms characterized by small normal and temporal baselines, allowing to maximize spatial and temporal coherence and therefore the quality of the interferograms (maximum temporal baseline 18 days, multilooking of 7:2 for range and azimuth direction, Goldstein lter window size of 32 pixels, coherence threshold of 0.2, a ground resolution of 30 m). The topographic phase contribution was subtracted using the ArcticDEM digital terrain model of 2 m spatial resolution. InSAR time-series were compared to data from a GNSS station provided by the IMO and available through the University of Iceland 30 , showing good agreement for both ascending and descending dataset (Fig. S2 a and b). Ascending and descending LOS displacement maps were combined to derive the vertical and the EW components of the ground motion (Fig. S5).

Gravity data processing
Each of the four gravity campaigns (Fig. S6) lasted for a few days. Elevation and gravity at the second northernmost station (HS22) were used as a reference, but free-air corrected of up to 2.6 μGal to account for small elevation changes at this site between campaigns. The gravity value at the end stations (HS16 and SNH25) were less constrained or had di cult local conditions due to snow and ice and were therefore regarded as unreliable. Finally, the gravity data were corrected for tidal and latitude effects and free air gravity computed using elevation changes obtained from InSAR measurements. The data were not corrected for ocean tides, as these effects can be considered negligible.
The uncertainty in similar gravity campaigns is generally of the order of 10-15 μGal for individual data points. These error limits are, however, highly dependent on external conditions like weather prior to and during the measurements, the local site condition, the nearby anthropogenic activity, or seismic noise. We selected quiet and calm days to minimize the measurement uncertainty. Our dataset provides consistent maximum values close to the centre of uplift which decrease with distance. However, a few data points are outliers, both with respect to the general trend with time and the nearby stations. It is our estimate that, apart from obvious outliers and with respect to internal consistency of the data, the error limits are close to ± 5 μGal at the centre of uplift.
We corrected for long term background gravity reduction caused by the net production of geothermal uid at Svartsengi 31 . The ratio of the total gravity decrease from 1976 to 2014 to the corresponding net mass production, gives a gravity decrease of 0,67·10 -9 μGal/kg. By applying this to the estimated production of 2020 we get a value of 5.1 μGal/year at the centre of uplift.
Finally, we estimated the gravity increase during the rst week of uplift in January 2020 prior to the rst gravity campaign. At that time, the central uplift was already 4 cm. Computing the ratio between the free-air corrected gravity increase and total uplift between the rst and the second gravity campaigns reveals 170 µGal/m resulting in 7 µGal gravity increase during the rst week of uplift (Tab. S1).

Poroelastic modelling
To explain the observed surface uplift, we rst tested the widely used Mogi model 18,32 which describes a pressure increase within a spherical chamber in a homogeneous elastic half-space. The best-t source depths are found to be 4.0, 4.4 and 4.9 km for the 3 uplift episodes, respectively. Though this simple elastic model can reproduce the cumulative surface uplift in space (Fig. S3), it, has several issues: 1) The subsidence motion after each uplift episode seems too fast and large to be explained by the cooling process of the intruded magma (Fig. 2d). 2) The induced change in Coulomb stress is concentrated directly around the assumed magma chamber and therefore cannot explain the triggered seismicity, which has been observed in a widespread area (Fig. 3). 3) Using a realistic estimate of magma density (~2700 kg/m 3 ), it is di cult to explain all campaign gravity satisfactorily with the episodic magma intrusions based on the deformation data. 4) Most importantly, there is no indication for the existence of a magma chamber under the Svartsengi uplift area. Actually, anomalous low Vp/Vs ratios are observed in the lower crust beneath Svartsengi 25,26 indicating gas-saturated porous rocks rather than partial melting. Therefore, we present a more realistic poroelastic diffusion-deformation model using the semi-analytical tool POEL 33 , which simulates strongly coupled diffusion-deformation processes induced by injection tests. Based on Biot's linear poroelasticity theory, a poroelastic medium is de ned usually with 5 parameters, i.e., the shear rigidity µ, the drained Poisson ratio ν, the undrained Poisson ratio ν u , the Skempton coe cient B and the hydraulic diffusivity D. In the rst step, we adopt the 1D seismic SIL reference model (Tab. S3) 34 . As the diffusion effect might be negligible in the seismic frequency band, µ and ν u can be derived simply from the density and seismic velocities given in the SIL model. Except for a thin permeable aquifer, all other layers are assumed to be nearly impermeable (D=0.001 m 2 /s). Additionally, we assume that the whole half-space medium is fully saturated (B=0.99) and has a standard drained Poisson ratio (0.25). The depth z, and poroelastic parameters ν u and D of the aquifer layer were optimized to best t the InSAR data. Note that the thickness of the aquifer can be xed (here 10 m), as the surface deformation is only sensitive to the aquifer transmissivity, which is proportional to the product of thickness and hydraulic conductivity.
The source is represented by episodic uid injection into this thin and con ned aquifer. The central location of the injection coincides with the maximum surface uplift at 63.870°N/22.465°W. Each of the 3 uplift episodes is attributed to a period with constant injection rate. The timing of the 3 injection periods is estimated directly by the upliftsubsidence turning points of the InSAR time series with a resolution of 6 days, i.e., Episode 1 from day 12 to 42, Episode 2 from day 60 to 108, and Episode 3 from day 132 to 192 of year 2020. Note the surface uplift exhibits a slightly elliptic pattern, which can be interpreted reasonably by a dominant diffusion along the main strike of the Svartsengi HT eld (Fig S1). Therefore, we use a Gaussian line source instead of a point source. The orientation and the characteristic length of the source was optimized, too. Note that the parameters of the aquifer and the source geometry have a nonlinear control on the surface deformation, while the 3 constant intrusion rates are linearly related with the surface uplift. We optimize the nonlinear parameters by a Monte Carlo search. For each set of given nonlinear parameters, we estimate the 3 intrusion rates by the least-squares method. The nal best-t model is de ned by the minimum mis t. For the quasi-static induced poroelastic deformation, only the volume, rather than the density, of the intruded uid is relevant. The gravity measurements provided an independent constraint on the density of the intruded uid. The internal deformation has no effect on the free-air gravity change for spherical pressure sources in a homogeneous halfspace 35 . Based on this, and on predictions from superposition of spherical pores 36 , we suppose the internal deformation is a secondary effect and only consider the primary effect of the intruded mass on the free-air gravity change. For this purpose, we use the POEL output for Darcy ux to calculate the spatio-temporal redistribution of the intruded mass within the aquifer layer and its contribution to the gravity variations on the surface. A simple comparison between the predicted and observed free-air gravity anomalies yields a uid density estimate of about 840 kg/m 3 . The mis ts of the model to the data are all within the measurement uncertainties (Fig. S7).

Joint analysis of DAS and a local seismic network by waveform attribute stacking
The DAS strain-rate data have been recorded using an Silixa iDAS (version 2) interrogator with a sampling frequency of 1 kHz, a gauge length of 10 m and spatial channel offset of 4 m along the full length of the bre. The DAS data has been included as 20 virtual single-component stations by extracting the DAS signal at regular 500 m intervals along the bre. For each virtual station, the signal was integrated over 36 m (9 channels) along the cable to improve the signal to noise ratio.
The employed automatic detection method 37 is based on an image function (IF) computed for a grid of potential source positions and times. The IF is computed from the stacking of time-back-shifted waveform attributes of P-and Sphases using the local velocity model (Tab. S3). Waveform attributes of P-phases employ a smoothed STA/LTA function 38 calculated from ltered seismograms. S-phase attributes are calculated from a smoothed squared signal which is more sensitive to phases of highest amplitudes. Smoothing causes the detector to be more robust against errors in the assumed seismic velocities and reduces the computational cost. Waveform attributes are normalized by a moving average with a duration longer than the duration of a typical transient event. Seismic events appear as spatiotemporal peaks in the 4D image function. A detection is registered when the IF exceeds a certain threshold value. The position of the spatio-temporal peak is used to provide the origin time and location. The detector is implemented in the Lassie software package, distributed under the Pyrocko framework 39 (Tab. S2).

Seismicity model
Earthquake triggering is expected for positive changes of the induced CFS, while no triggering is expected for negative changes 40,41 . For simplicity, we assume that the number of earthquakes is proportional to positive CFS-changes and consider stress shadowing (Kaiser effect) by only allowing triggering if the absolute stress exceeds all precursory values 42,43 . Thus the predicted cumulative number of induced events becomes N(x, t) ~ max (CFS(x, time<=t)).
The CFS-change is given by S(x, t) = τ max (t) + f p(x, t), where f is the friction coe cient, τ max and p are the induced maximum shear and pore pressure, respectively, based on the poroelastic model described above. Here we use the maximum shear value to account for variable receiver mechanisms. However, we nd that the results are robust if we calculate CFS for a xed right-lateral NS strike-slip mechanism. We calculated the CFS in a spatial box with a dimension of 30 km in NS-and 60 km in EW-direction, centered at the geodetically inverted deformation source, with a grid spacing of 0.5 km and at ten depth levels from 0.5 km to 9.5 km.
The susceptibility to stress changes is assumed to be proportional to the fault density estimate, based on preceding seismicity, so-called background activity. For that purpose, we analyzed the IMO-catalog using M≥1.5 events that occurred in this region between 2000 and 2019. Note that the catalog is complete for this magnitude cutoff. The rate of M≥1.5 background events in the selected region is found to be r=0.76 days -1 , which translates to r = 2.41 M≥1 events days -1 assuming a Gutenberg-Richter magnitude distribution with b=1.
The epicenters are laterally smoothed using a Gaussian lter with a standard deviation of 1 km, and the depth distribution is approximated by a Gaussian distribution with a mean of 5.7 km and a standard deviation of 1.8 km. After normalization, this resulted in a probability p r (x i ) for the occurrence of an event in a given 3D-grid cell x i .
The spatiotemporal rate of expected events is simply equal to R(x i , t) = r p r (x i ) + c p r (x i ) d/dt CFS(x i , t) for positive stress changes and CFS exceeding its preceding maximum value and R(x i , t) = r p r (x i ) else. The model has two parameters, (i) the friction coe cient f, which has typical values in the range from 0.1 to 0.8, and (ii) the susceptibility parameter c. We nd that f only slightly affects the results and set it to f=0.7 44 . The parameter c is set by the condition that the predicted total number of events in the space-time volume is equal to the number of observed earthquakes.

Cyclic injection model
We explain the unrest at the Svartsengi HT eld by at least two interacting processes: (1) continuous, slow in ow of magma from the mantle into the base of the lower crust and (2) cyclic rise of low-density pressurized uids to a depth of 4 km, along the updoming BDB. The low-density uids could be magmatic volatiles such as CO 2 . The volatiles either rise from melting source at greater depths or exsolve directly from the magma that has accumulated at the BDB.
The drainage of a deep reservoir by viscous channel ow is similar to the gravity-driven out ow from a tank through a pipe. If viscosity is dominant, the ux follows an exponential decay with q(t) = q 0 e -t/t0 , where q 0 is the initial ux at time t=0 with t 0 = r 2 h 2 (t=0)/q 0 (Fig. 5A & information sheet IS1). The ascent height h or the length L of the inclined channel do not enter the problem. The modeling of the total in ux after three uplift episodes leads to a draining decay time of t 0 = 115 days and an initial ux of q 0 = 3.2 m 3 /s (Fig. 5b).
The cyclic ascent of nite-volume uid batches can be explained by the fracture-mechanical concept of buoyancydriven uid ascent. The process is controlled by the lling of uid batches with ux rate q(t), their detachment from the feeding reservoir, and the ascent through a channel or fault to the out ow point at 4 km depth (injection into the sealed aquifer). The problem can be illustrated with the transport of continuously arriving packages with an "elevator system".
Each elevator is lled to its maximum load before it departs. The volume of intrusion batches at the upper aquifer has been estimated at ~34 10 6 m 3 , ~33 10 6 m 3 and ~26 10 6 m 3 for the rst, second, and third batch, respectively (Fig 5b). A decrease of the batch volume can be understood from the continuous reduction of the fracture roughness of rock the more often a batch rises. If the ux of arriving packages is high, the time to load the elevator is short. If the ux continuously decreases, the loading time increases. The ascent velocity of the uid batch is almost independent of the ux at the '' lling station'' and only determined by the properties of the channel and the enclosed uid volume. A constant ascent velocity is reasonable from a fracture mechanical approach and e.g. supported by laboratory experiments of nite volume uid ascent in gelatine 45 .
The model explains the observed time function of the accumulated, injected volumes, the onset and inter-event time of the three uplift cycles (Fig. 5 b & c). The t has been obtained with the requirement that the accumulated volume at the end of phase 3, ΣV = V 1 +V 2 +V 3 , and the onset of the rst injection phase are explained. The model is also consistent with the onsets of the uplift phases 2 and 3, and maybe a fourth phase at Krýsuvík in September 2020 (Fig. 5c & Figs. S10 & S11). What the model cannot resolve is the ascent times of the uid batch after the detachment from the reservoir to the injection point at the aquifer at 4 km depth. Therefore, all times measuring the beginning of a new uplift phase in Fig. 5b and c are arrival times of the uid batch at 4 km depth. For instance, the start of the lling of the rst uid batch is estimated for the 18 th of December 2019 (Fig. 5c), which correlates with the occurrence of the strong earthquake swarm from 15 th -20 th of December 2019 at 5-6 km depth and about 7 km east of the centre of uplift. Therefore, we suggest that the in ow in the channel occurred at this location. This would also link the unrest beneath Svartsengi to the eruption that started in March 2021 just above the December 2019 earthquake swarms (Fig. S9). The increasing duration of the injection into the aquifer indicates that the pressure difference between the head of the uid batch and the sealed aquifer is decreasing with every arrival of new over-pressurized batch. This is expected if the vertical height of the batch (its volume) decreases but may also indicate that the aquifer is more and more pressurized by the in ow of uids. Consequently, the layer above the aquifer was uplifting and bending. Parts of the overpressure in the aquifer diffused away, but not all before the arrival of the next batch.
Declarations Figure 1 Overview