# Deformation and seismicity decline before the 2021 Fagradalsfjall eruption

### Routine earthquake places and seismicity fee

Seismicity within the space is monitored by the Icelandic Meteorological Office (IMO), the place earthquakes are recorded on the SIL nationwide seismic community. Arrival occasions of P and S waves from the earthquakes are mechanically detected by the SIL evaluation system, which locates the occasions and assigns a preliminary magnitude. The occasions are then manually reviewed, revised and up to date as required. These routine analyses use one frequent one-dimensional velocity mannequin (SIL mannequin) for all occasions in Iceland42,43. The whole variety of earthquakes mechanically detected by the SIL system in the course of the interval 24 February to 14 March 2021 was about 45,000, with 15% manually reviewed. In the next interval, 15–19 March, 19% of the roughly 7,500 mechanically situated earthquakes have been reviewed. To enhance seismic monitoring of the unrest on Reykjanes Peninsula, eight extra stations on the peninsula, owned by the Czech Academy of Sciences (CAS) and operated by Iceland GeoSurvey (ÍSOR), have been related to IMO’s monitoring system and are additionally used within the computerized and reviewed customary places. Three of those stations (lsf, iss and lag) have been added to the system in 2020 and an extra 4 (lat, ash, moh, faf) have been added in early March 2021 (the primary three on 5 March, the fourth on 11 March). On 5 March one station (odf) from University of Cambridge was additionally included within the routine computerized and guide analyses. The lately elevated station density on Reykjanes Peninsula has improved the magnitude of completeness (Mc) for seismicity within the space, from its earlier 0.5 or higher estimate for the interval 2002 to 2013 (ref. 44). However, throughout intense seismic swarms with massive earthquakes the automated location system saturates and background noise vastly will increase, leading to vital lower in sensitivity. For the extraordinary interval analysed right here, Mc for Reykjanes Peninsula is estimated to be round M1, and thus the seismicity fee evaluation is proscribed to M > 1 earthquakes. Seismicity fee (Fig. 1c) is estimated individually for the dyke and the plate boundary areas, by summing up all occasions in a 2-km-wide strip across the dyke and a 2.5-km-wide strip across the plate boundary axis (earthquakes within the overlap space are divided between the dyke and the plate boundary areas). Figure 1c additionally experiences variety of earthquakes in different areas. With solely a fraction of the seismicity having been manually reviewed, the seismicity evaluation relies on computerized places.

### Earthquake focal mechanism and magnitude willpower

After guide evaluation of the earthquake places, focal mechanisms and magnitudes are calculated, and the occasions added to the SIL catalogue. The calculation relies on a grid search via all potential strike, dip and rake angles for which first-motion P-wave polarities match and for which P- and S-wave amplitudes are inside acceptable limits of these noticed45. Magnitude calculations within the SIL evaluation system are likely to underestimate magnitudes for occasions of MW > 3.5 (ref. 46). Better magnitude estimates for these occasions are mechanically calculated by a separate real-time evaluation, which repeatedly evaluates noticed peak floor velocity (PGV) in several frequency bands at every seismic station47. With these values and computerized places from the SIL system for occasions of M > 3.5 the second magnitude, MW, is calculated for every station utilizing a ground-motion prediction equation for velocity48. The last MW estimate for every occasion is a strong imply of a 3rd of the calculated values across the median, from stations outdoors the clipping vary (round 15 km) to a distance of 230 km from the occasion49. This worth is powerful and compares properly with internationally decided second magnitudes. Supplementary Table 1 lists all the decided magnitudes for MW ≥ 4.0 occasions on Reykjanes Peninsula from February to April 2021. These magnitudes are used to calculate the cumulative second launched within the largest occasions in the course of the unrest interval.

### Relative relocation of earthquakes

The seismic exercise in the course of the interval 24 February to twenty March 2021 is situated on and round each the dyke intrusion and the numerous roughly north–south strike-slip faults triggered by the intrusion. Figure 1 and Supplementary Table 1 present that 15 strike-slip occasions of M ≥ 4.5 have been triggered close to the dyke, leading to many aftershocks on their fault planes within the magnitude ranges above round M1.5. Therefore, in an effort to map the general options and dimensions of the dyke intrusion and to attenuate the contamination by aftershocks from the roughly north–south strike-slip faults, an information set of smaller magnitude occasions was chosen for high-precision relative relocation. The information set contains all 1,321 earthquakes to date routinely situated and reviewed, within the native magnitude vary 0.5 < M < 1.4.

The relocation methodology makes use of double variations of each ‘picked’ absolute arrival occasions and improved relative arrival-time measurements obtained via cross-correlation of P and S waveforms from completely different earthquakes50. The process iteratively inverts the weighted sq. sums of absolute P and S arrival-time variations, in addition to the double variations of (1) absolute arrival occasions of P and S waves, (2) relative arrival occasions of P and S waves and (3) relative S–P arrival occasions. The ensuing distribution of occasion places has a excessive inside/relative location accuracy and is due to this fact properly suited to map particulars of energetic subsurface faults and dyke intrusions51,52. The velocity construction on the Reykjanes Peninsula is considerably completely different from the usual SIL-velocity mannequin used within the routine evaluation, influencing particularly the hypocentre depths, which are usually roughly 1 km too deep within the routine evaluation. Therefore, a better-fitting one-dimensional velocity mannequin is used, derived from travel-time and phase-velocity evaluation of physique waves propagating alongside the Reykjanes Peninsula53 and the Reykjanes-Iceland Seismic Experiment seismic refraction profile54. Fifteen seismic stations inside a distance of fifty km from the earthquakes are used within the relocation together with all IMO stations and the three CAS/ÍSOR stations related to the system in 2020. Exclusion of stations closest to the dyke is intentional to fulfil the necessities for the relocation methodology to correctly work, because it assumes that waves from carefully situated occasions journey via the identical velocity construction. The ensuing earthquake hypocentre depths are referenced to the imply elevation of the closest stations, at roughly 80 m above sea degree.

The ensuing relocated occasion distribution (Fig. 3 and Extended Data Fig. 3) exhibits two primary segments within the roughly 9-km-long dyke intrusion, differing in strike and depth vary. A northern section, at the least 4.5 km lengthy, has a extra easterly strike, of round 45°, and reaches depths of 8 km. A 4.5-km-long southern section strikes 25° and reaches depths of 6 km. The change in depth is slightly sharp on the intersection of the 2 segments (B on Extended Data Fig. 3). Extended Data Fig. 3 additionally exhibits a clustering of occasions round a depth of 0.5 km within the eruption area (the star and white line on the determine). The color coding of occasions in Extended Data Fig. 3 reveals the general time evolution of the exercise in the course of the interval. Seismicity begins at a depth close to the centre of the dyke (B on Extended Data Fig. 3) and prompts the northern section between 24 February and three March, earlier than propagating southwards activating the southern section from round 4 March to 19 March. The space across the intersection (B on Extended Data Fig. 3) stays energetic for the entire time interval.

### Geodetic modelling

We inverted the geodetic information utilizing a Bayesian strategy with a modified model of GBIS software program27 to acquire the vary of supply parameters that may clarify the noticed deformation. An preliminary inversion was carried out utilizing GNSS observations from 18 stations (a mix of marketing campaign and steady measurements) overlaying the whole pre-eruptive interval from 23 February to 19 March 2021, and two Sentinel-1 interferograms (from descending monitor T155 overlaying the interval 23 February to 19 March 2021 and ascending monitor T16 overlaying the interval 19 February to 21 March 2021). This Bayesian inversion software program makes use of a Markov-chain Monte Carlo strategy and the Metropolis Hastings algorithm55,56. The joint likelihood distribution operate for the varied supply parameters was obtained by operating two million iterations. In our preliminary inversion we modelled the dyke with an oblong dislocation57, permitting for opening and shear, and the noticed motion on the central axis of the plate boundary with a second dislocation, additionally permitting for opening and shear. We modelled the floor displacements ensuing from the 2 largest tectonic earthquakes as separate rectangular dislocations, permitting for shear movement solely. A subsidence sign near the eruption web site was modelled as a deflating level supply58 for simplicity, though this deflation could also be the results of regular faulting associated to the dyke intrusion. We assumed a Poisson’s ratio of 0.27. All of the supply parameters have been allowed to range, with broad bounds outlined a priori.

The output parameters from this inversion confirmed {that a} dyke was intruding alongside vertical planes outlined by the seismicity. In a second inversion, we substituted the only dyke dislocation for 2 related dislocations with strikes mounted to that indicated by the relocated seismicity (N45°E within the northern section and N23.5°E within the southern section). We allowed the place to range, recognizing that absolutely the places for the relocated seismicity aren’t correct, and located that the popular location was round 400 m west of the relocated seismicity and handed via the eruption fissure.

For our last inversion we mounted the placement of the dyke to the popular location from the second inversion with most lengths (4.5 km for every section) and depths (7.5 km within the north and 6 km within the south) derived from the seismicity. We additionally mounted the placement, dip (vertical) and strike (N65.4°E) of the plate boundary dislocation from the second inversion. Very little opening (9 cm) was assigned to this section within the second inversion, so we allowed for shear movement solely. Given the vertical nature of each the dyke and plate boundary segments, we allowed for strike-slip shearing in a strike-slip sense solely, with no dip-slip movement. We divided the dyke and plate boundary section into patches, 750 × 750 m2 within the higher 3 km and 1.5 × 1.5 km2 beneath a depth of three km, and solved for the opening and/or slip independently on every patch. We allowed the placement of the deflating level supply to range, and our mannequin positioned it at a depth of 1,480–1,670 m with a quantity change within the vary of −2.6 to −3.4 × 106 m3.

To constrain the each day incremental opening alongside the dyke, we used a constrained linear least squares strategy. We used 100 samples from the posterior distribution of our last inversion for the whole pre-eruptive interval to constrain the full opening for every dyke patch and solved for incremental opening utilizing GNSS solely. To scale back the noise within the GNSS time collection, we first filtered them utilizing a 2-day triangular transferring window, earlier than calculating the incremental displacements. We estimated the each day opening for every patch, with a constraint that solely constructive opening was allowed. To account for systematic errors affecting the entire GNSS community, we concurrently solved for a each day reference body adjustment in east, north and up. We then summed the each day quantity change of each patch to offer a each day quantity change fee and calculated the imply depth (centre of gravity) of this each day quantity change.

The trade-off between each day quantity change and imply depth of magma emplacement was studied with a simulation check. An incremental dyke opening inversion was carried out on the simulated information to check whether or not or not the noticed trade-off between depth and quantity change for the each day options (Fig. 3) influences the general relationship between depth and quantity change. We divided the full opening of the utmost a posteriori likelihood answer from our distributed opening inversion into 28 days of equal quantity change, with shallowing imply depth. We used these each day dyke opening fashions to simulate the displacements on the GNSS stations and added consultant noise. We then inverted the simulated GNSS information in the identical manner as earlier than. The outcomes are proven in Extended Data Fig. 8. Despite the trade-off between depth and quantity for particular person days, the simulated fixed quantity change fee is retrieved, with no dependence on depth.

Ground deformation attributable to all 64 recorded earthquakes Mw ≥ 4, apart from the 2 largest occasions on 24 February and 14 March, was evaluated by the next process. Each earthquake was assumed to be on an oblong north–south putting vertical strike-slip fault in accordance with the dominating fault mechanisms. The central location and depth of every fault was acquired from the IMO hypocentre catalogue. For all earthquakes smaller than M5, an oblong fault of two × 2 km2 was assumed; for M ≥ 5 a size of 4 km and a peak of two km was used for the earthquake. For the mapping of very shallow hypocentres (with depths of lower than 1.1 km), 1 km was added to the central depth to keep away from artefacts from a man-made fault protruding of the bottom. The second magnitude was transformed to seismic second M0 utilizing the relation M0 = 101.5M+9.05 (ref. 59) and a shear modulus of μ = 30 GPa to acquire the imply slip s from the seismic second, M0 = μAs, the place A is the fault space. The cumulative second of all 62 faults (that’s, excluding the 2 largest earthquakes) was equal to a single MW 5.9 earthquake. The software program Coulomb60 was used to calculate the floor deformation from the faults in a dense (0.1 × 0.1 km2) grid. The ensuing displacement subject (Extended Data Fig. 9) has the identical traits because the geodetic mannequin produced by the primary inversion utilizing distributed slip alongside the dyke along with the plate boundary section and the 2 primary faults, however is of a lot smaller magnitude. Considering the contribution of smaller faults within the inversion course of as described above didn’t enhance the match to the geodetic information, so these aren’t thought of right here. The results of the smaller faults might be anticipated to be already accounted for in the primary inversion.

### Conduit stream modelling

We modeled the stream to the bottom of the dyke via a cylindrical conduit. In actuality, the conduit could also be a unique form and extra fracture-like when it’s forming. However, it’s fashioned within the ductile decrease crust beneath the locking depth of the crust, the place plate tectonic stress doesn’t construct up in the identical method as within the topmost elastic crust. Furthermore, magma stream in a fracture of restricted dimensions might quickly focus in direction of the widest a part of such a fracture, and the efficient stream path might develop into extra cylindrical in form32,61,62. In any case, a non-circular cross part wouldn’t considerably have an effect on our conclusions (see beneath). Laminar stream of an incompressible fluid within the conduit is assumed, pushed by the distinction between the stress on the magma supply area and the magmastatic head, and by treating the dyke as a reservoir that presents negligible resistance to stream,

$$P-{rho }_{{rm{m}}},gleft(L+hright)=frac{8nu L}{{r}^Digital}u+frac{{rho }_{{rm{m}}}}Digital{u}^Digital,$$

(1)

the place P is the stress on the magma supply, ρm is the density of the magma, (g) is gravitational acceleration, L and r are the size and radius of the conduit, respectively, h is the peak above the bottom of the dyke at which magma is intruding, ν is the dynamic viscosity and u is the imply magma velocity within the conduit. The time period on the left is the driving stress, the primary time period on the appropriate is the viscous loss attributable to laminar (Hagen–Poiseuille) stream and the second time period on the appropriate is the dynamic stress loss. Using the connection between volumetric stream fee and velocity, (Q={rm{pi }}{r}^Digitalu),

$$P-{rho }_{{rm{m}}},g(L+h)=frac{8nu L}{{pi r}^Media}Q+frac{{rho }_{{rm{m}}}}{2{{pi }^Digitalr}^Media}{Q}^Digitalapprox frac{8nu L}{{pi r}^Media}Q,$$

(2)

within the case the place (Lgg {rho }_{{rm{m}}}{rm{Q}}/16pi nu ). Note that within the case when h = 0 and (P={rho }_{{rm{c}}}gh), the place ρc is the common density of the crust above the magma supply area, this reduces additional to the extra acquainted relationship (Q={{rm{pi }}r}^Medialeft({rho }_{{rm{c}}}-{rho }_{{rm{m}}}proper)g/8nu ) (ref. 63). Differentiating (2) with respect to h offers

$$-{rho }_{{rm{m}}},g=frac{8nu L}{{{rm{pi }}r}^Media}frac{{rm{d}}Q}{{rm{d}}h}.$$

(3)

In our conceptual mannequin, we assume that the stress to carry the dyke open is supplied by the deviatoric stress within the higher crust and the dyke is successfully a tank being crammed by a conduit coming into its base. Dividing (P) into the stress as a result of weight of the crust above and beneath the bottom of the dyke, and contemplating the stream fee into the bottom of the dyke, QD, from (2)

$${rho }_{{rm{u}}},gD+bigtriangleup rho gL=frac{8nu L}{{{rm{pi }}r}^Media}{Q}_{D}$$

(4)

the place ρu is the common density of the higher crust (above the bottom of the dyke), D is the depth to the bottom of the dyke and Δρ is the common distinction between the densities of the decrease crust and the magma.

Combining (3) and (4) offers

$$frac{{rm{d}}Q}{{rm{d}}h}=frac{-{rho }_{{rm{m}}}{Q}_{D}}{{rho }_{u}D+bigtriangleup rho L}.$$

(5)

Thus, if we are able to measure the influx fee on the base of the dyke and the change of influx fee with peak within the dyke, we are able to decide the size of the conduit:

$$L=frac{-{rho }_{{rm{m}}}{Q}_{D}}{bigtriangleup rho frac{{rm{d}}Q}{{rm{d}}h}}-frac{{rho }_{{rm{u}}}D}{bigtriangleup rho }.$$

(6)

Assuming values of ρm= 2,700 kg m−3, ρu = 2,700 kg m−3, Δρ = 300 kg m−3 and D = 6 km, and values estimated from our information of  QD = 32 m3 s−1 and (frac{{rm{d}}Q}{{rm{d}}h}) = −0.0043 m2 s−1 offers L = 13 km, which places the supply depth at an estimated 19 km. Note that this result’s unbiased of the viscous stress loss time period and doesn’t due to this fact depend upon the conduit being cylindrical in form (each Q and dQ/dh depend upon the viscous time period, which cancels when taking the ratio). Also be aware that the mannequin doesn’t require that ρm is the same as ρu. Assuming g = 9.8 m s−2 and ν = 100 Pa s then offers an estimate for (r) of 0.9 m. We thus infer that the cross-sectional space of the conduit is of the order of some sq. metres. If the efficient stream path is just not corresponding to that of a cylindrical conduit, then its cross-sectional space might be anticipated to be of comparable magnitude, though viscous drag is determined by the conduit form58.