Introduction

The Reykjanes Peninsula in south-west Iceland is an oblique rift zone forming part of the divergent plate boundary between the North American and Eurasian plates and the onshore continuation of the Mid-Atlantic Ridge. Rifting is accommodated within six volcanic systems, arranged in an en-echelon pattern along the plate boundary (Fig. 1a) (Sæmundsson and Sigurgeirsson 2013; Sæmundsson et al. 2020). These mark the areas with the highest density of tectonic fractures and eruption fissures. The intrusion of magma in NE-SW oriented dykes accommodate an extensional component of the rifting, with the remaining strike-slip component accommodated by N-S oriented strike-slip faults, thought to be capable of producing magnitude 6 earthquakes (Einarsson 1991; Björnsson et al. 2020). Dating of Holocene lava flows has shown that volcanism and rifting along the Reykjanes Peninsula is episodic, with most of the volcanic systems experiencing events every 800–1000 years for the last 4000 years, except Fagradalsfjall, which has not experienced volcanism for 6000 years (Sæmundsson et al. 2020). The volcanic activity seems to affect each system sequentially, jumping successively from east to west, with the jumps spaced 100–200 years apart. Seismicity along the Reykjanes Peninsula is also episodic, with earthquake swarms observed at intervals of around 30 years during the last century of instrumentation (Einarsson 2008).

Fig. 1
figure 1

Location of the study area and seismic stations used in this study. (a) Shows the Reykjanes Peninsula with volcanic systems shaded orange and mapped fractures delineated by thin black lines (Clifton and Kattenhorn 2006; Sæmundsson and Sigurgeirsson 2013). The location of seismic stations used in this study are indicated by the triangles and coloured by the dates the instruments were available. The volcanic systems are named: R—Reykjanes, S—Svartsengi, F—Fagradalsfjall, K—Krýsuvík, B—Brennisteinsfjöll and H—Hengill. The region outlined by the purple dashed line is shown in Fig. 5. The inset shows the location of the Reykjanes Peninsula (red box) in relation to Iceland. (b) Shows the region within the black dashed rectangle in the upper panel. The new lava erupted during 2021 is indicated by the orange polygon and the vent locations by the orange circles (Pedersen et al. 2022). The names of stations mentioned in the text are labelled

Full size image

The 2021 Fagradalsfjall eruption was the first on the Reykjanes Peninsula for about 800 years. It was the culmination of a sequence of seismic swarms and intrusion events along the western part of the peninsula that started in 2020 (Cubuk-Sabuncu et al. 2021; Geirsson et al. 2021; Flóvenz et al. 2022). In this study, we focus on the seismicity around Fagradalsfjall, recorded by a seismic network operated by multiple research groups. Due to the density of stations, we have obtained an exceptionally high-resolution image of the seismicity before and during the eruption in 2021.

This study investigates a cluster of seismicity located below the brittle-ductile transition (6–9 km). Seismicity at these depths has been observed before in Iceland (Key et al. 2011b; Greenfield and White 2015; Hudson et al. 2017) as well as elsewhere in the world (e.g. Reyners et al. 2007; Shelly and Hill 2011; Wech et al. 2020) and is thought to be due to fluids (either CO2 or magmatic) moving through the lower crust. Using automated detection algorithms and detailed relocation of the seismicity, we present a timeseries of deep seismicity preceding the 3-week-long dyke intrusion before the eruption started. During the Fagradalsfjall eruption itself, we compare the occurrence of deep earthquake swarms to changes in the nature of the eruption (i.e. discharge rate or fountaining style). In doing so, we show that the deep seismicity is likely linked to the magmatic plumbing system stretching from near-Moho magma storage bodies to the surface.

Tectonic, crustal and seismic background

Crustal structure

Using seismic refraction techniques (RISE survey, Weir et al. 2001), the crust beneath the Reykjanes Peninsula can be split into two layers: an upper layer with a steep velocity gradient and a lower layer where the velocity gradient is low. The Moho is located beneath this lower layer, at a consistent depth of 15 km (Weir et al. 2001). The upper layer is thought to be equivalent to oceanic layer 2 and corresponds to sheeted dykes, sills and lava flows. Cumulate gabbros form the bulk of the lower crust (oceanic layer 3) (Brandsdóttir and Menke 2008). Both layers are thickened relative to typical oceanic crust, but are formed through similar processes.

Little is known about the 3D velocity structure of the Reykjanes Peninsula. Local earthquake tomography of the eastern Reykjanes Peninsula and the South Iceland Seismic Zone have suggested that the plate boundary is characterised by reduced velocities between 6 and 7 km depth relative to the surrounding area (Tryggvason et al. 2002). More recent surface wave tomography results using ambient noise observations indicate that the upper crust is highly heterogeneous, with low-velocity bodies either correlating with known high-temperature regions or sites of past volcanism (Martins et al. 2020). Electromagnetic imaging methods and seismic observations have revealed that the shallow crust beneath the Reykjanes Peninsula is similar to that observed in high-temperature geothermal fields elsewhere in Iceland except for a lack of deeper conductive bodies (Hersir et al. 2020; Karlsdóttir et al. 2020) and shallow bodies which attenuate shear-waves (Brandsdóttir and Menke 2008). From lithological studies, it has been suggested that the heat source of the Reykjanes high-temperature geothermal fields are repeated and cooled magmatic intrusions (Friðleifsson et al. 2020).

Magmatic plumbing systems

Magmatic plumbing systems beneath Icelandic volcanoes are generally thought to be complex, with multiple sills and intrusions distributed through the lower crust (Maclennan et al. 2001; Cashman et al. 2017; Maclennan 2019). Magma batches are thought to accumulate in storage reservoirs and can reside at depth for thousands of years before eruption or intrusion (Mutch et al. 2019a). During active periods, magmas can ascend rapidly through the crust (Mutch et al. 2019b) and reactivate previously intruded bodies (Tarasewicz et al. 2012). Magma ascent, especially to shallow depths, is focussed beneath central volcanoes (Greenfield et al. 2016), although a significant proportion can occur away from these areas (Key et al. 2011a, 2011b; Hartley and Thordarson 2013).

Volcanic deep long period events

It is estimated that between 5 and 20% of the melt generated beneath mid-ocean ridges is erupted at the surface (White et al. 2019). Thus, insight into magma intrusions can help reveal the process by which the crust beneath Iceland is built. In this study, we investigate a cluster of low-frequency earthquakes below the brittle-ductile transition beneath the Reykjanes Peninsula. Earthquakes such as these are termed deep long period events (DLPs) and have been observed beneath volcanic regions around the world (Nakamichi et al. 2003; Nichols et al. 2011; Shapiro et al. 2017; Hensch et al. 2019; Wech et al. 2020; Kurihara and Obara 2021) including other volcanoes in Iceland (Greenfield and White 2015; Hudson et al. 2017). DLPs share some similarity with long-period (LP) earthquakes observed in the shallow crust (Chouet and Matoza 2013), although their greater depth (10–45 km) precludes many of the mechanisms invoked for shallow LP events. Low-frequency earthquakes (LFEs) and tectonic tremor beneath subduction zones also share some similarities with DLPs (Shelly et al. 2007; Peng and Gomberg 2010; Shapiro et al. 2018) and the mechanisms invoking the transport of fluids through the lower crust/upper mantle are often used to explain the occurrence of both DLPs and LFEs. The diversity of DLPs observed around the world (Chouet and Matoza 2013) or even beneath a single volcano (Nakamichi et al. 2003; Kurihara and Obara 2021) hints at a range of mechanisms which can cause DLP seismicity. In addition to the transport of either magmatic or volatile phases through the crust (Nakamichi et al. 2003; Shapiro et al. 2017), recent observations of long duration swarms of DLPs beneath both dormant and active volcanoes suggest mechanisms related to the build-up of fluid pressure from secondary boiling (Wech et al. 2020) or rapid exsolution of gas bubbles (Melnik et al. 2020).

Methods

Dataset

Our core seismic network consisted of a network of 3-component, broadband seismic stations deployed near the town of Grindavík (Fig. 1a) between June 2020 and August 2021. The network was originally deployed as 9 instruments extending as far east as the stations 8F.FAGN and 8F.FAGS on the western side of Fagradalsfjall (Fig. 1b). This was originally intended to monitor the ongoing seismic activity related to inflation within the Svartsengi volcanic system, triggering strike-slip fault movements (Geirsson et al. 2021; Flóvenz et al. 2022). The network footprint was extended eastwards in August 2020 and the station density around the Fagradalsfjall eruption site was increased in April 2021. The data from our network were supplemented by 20 instruments (18 broadband, two short-period) operated by Iceland GeoSurvey (ÍSOR), the Icelandic Meteorological Office, the Czech Academy of Sciences, the University of Washington and the University of Potsdam (Fig. 1). Sixteen of these instruments were active throughout the time period discussed here, one was active from March 2021 and three were active from June 2021. All the deployed stations were not permanent installations and some data was lost due to power issues. Station availability information is documented in the Supplementary Information (10.5281/zenodo.6136954).

One-dimensional velocity model

We calculate a bespoke one-dimensional (1D) velocity model for the region beneath our seismic network using VELEST (Kissling et al. 1994). We first manually picked P- and S-wave phase arrival times for a selection of well-located events extracted from a preliminary earthquake catalogue—located using a regional velocity model (Weir et al. 2001). Events fulfilling the criteria: maximum azimuthal gap < 180°; root-mean-squared residual time < 0.15 s; > 10 arrival time picks of which > 5 were S-wave arrival time picks; and minimum epicentral distance < hypocentral depth, were retained. This resulted in a catalogue of 253 events with raypaths covering the entire study area and 4447 arrival time picks (2359 P-wave, 2088 S-wave). VELEST inversion results can be sensitive to the starting model, so we first generated a suite of starting models defined by three variables: the P-wave velocity at the surface (1–4 km/s), the P-wave velocity (Vs) at 20 km depth (6–8 km/s) and the depth where the velocity gradient changes (5–10 km). S-wave velocities (Vs) were calculated using a constant Vp to Vs ratio (Vp/Vs) of 1.78. The model parameterisation was guided by typical 1D velocity profiles and Vp/Vs ratios observed both on the Reykjanes Peninsula and elsewhere in Iceland (Weir et al. 2001; Mitchell et al. 2013; Greenfield et al. 2020b). At this stage, the model was parameterised using 1-km-thick layers in the upper 8 km, and 2-km–thick layers from 8 to 20 km depth. We then selected the 50% best-fitting models as defined by the final data variance and calculated the mean model to use as a starting model for the next stage.

We used the mean model as a seed to generate 5000 new starting models with random depth parameterisation (layers ranging in thickness from 0.2 to 5.0 km) and velocities (± 2 km/s). We again selected the 50% ‘best’ models and calculated P- and S-wave models from the ensemble (Fig. 2). By using the ensemble to calculate velocity models, we allowed the data to guide the choice of where layers are required rather than prescribing them a priori. The final model was then calculated by using the ensemble model as input to a final run of VELEST. This ensured that our output model was at a minimum in misfit and allowed us to manually vary the damping parameters for individual layers (e.g. if they were poorly constrained by the data) and to invert for static station corrections. The final model is similar to that from RISE (Weir et al. 2001), but with slightly slower upper crustal velocities (Fig. 3). The relatively small number of earthquakes at depth (Fig. 2a) means that this velocity model is only constrained to ~ 11 km depth. At greater depths we use the RISE model. Interestingly, our data indicate a low-velocity-zone in the Vs model at the depth of the deeper earthquakes (> 8 km) (Fig. 3).

Fig. 2
figure 2

Results from velocity model inversion using VELEST. The ensemble of models from the range of starting models is indicated by the background shading with bright colours showing a high number of similar models. The mean model for each of the P-wave (a), S-wave (b) and Vp/Vs ratio (c) is delineated by the magenta line. The preferred model used in the next stage of the modelling and the model from RISE (Weir et al. 2001) with a Vp/Vs of 1.78 are indicated by the blue and red lines respectively. A depth histogram of the earthquakes used in the inversion is shown on the left panel and the number of rays which hit each layer in the final model is indicated by the white line on the middle panel

Full size image
Fig. 3
figure 3

Final inversion for the velocity model. The starting and final models are shown by the black dotted and solid lines respectively. The model from RISE (Weir et al. 2001) is shown by the red line

Full size image

Earthquake catalogue creation

We constructed our initial earthquake catalogue using QuakeMigrate (QM), an open-source Python package for automatic earthquake detection and location using waveform migration and stacking (Smith et al. 2020; Winder et al. 2021b, 2021a). This program uses continuously recorded waveforms to calculate a STA/LTA (short-term-average/long-term-average) onset function. Peaks in this function correspond to the arrival times of seismic phases and the width of each peak correlates with the pick uncertainty (Drew et al. 2013). To generate the earthquake catalogue for this study, we filter the waveform data between 2 and 16 Hz (to avoid the strong oceanic microseism on Iceland) and use 0.02 and 0.1 s long windows to calculate the STA and LTA respectively. The resulting STA/LTA functions for each station and component are migrated through a pre-calculated travel-time lookup table and peaks in the resulting coalescence function indicate hypocentres and origin times of events. A static detection threshold of 1.65 is used throughout the entire dataset. While this may miss events during quieter periods, we found that with the constantly changing number of stations and noise levels due to storms or the eruption, a relatively high and constant threshold gave more consistent results with fewer artefacts.

Local magnitudes for the events were obtained by first removing the instrument response from each waveform and then applying a standard Wood-Anderson response. Maximum zero-to-peak amplitudes were then measured for each horizontal component. Magnitudes were calculated by taking the logarithm of the amplitude and correcting for attenuation using a logA0 operator. We found that using standard logA0 operators (for Southern California, Hutton and Boore 1987) produced large residuals when applied to this data. Even logA0 curves estimated for other regions in Iceland (central Iceland, Greenfield et al. 2020b) did not produce satisfactory results. We therefore proceeded by calculating a new logA0 curve for our network by directly inverting for the logA0 fitting constants \(n\) and \(K\), the local magnitude, \({ML}_{i}\) and the station-component static corrections, \({C}_{jk}\), in the following equation (Keir et al. 2006; Illsley-Kemp et al. 2017) of the form \(Gx=y\),

$${\mathrm{log}}_{10}{A}_{ijk}+2={ML}_{i}-n{\mathrm{log}}_{10}\frac{{r}_{ij}}{17}-K\left({r}_{ij}-17\right)-{C}_{jk},$$
(1)

where \({A}_{ijk}\) is the measured zero-to-peak amplitude for event \(i\), station \(j\) and component \(k\) and \({r}_{ij}\) is the hypocentral distance between earthquake and station. We applied the constraint that all the station corrections should sum to zero and removed clipped observations. The resulting logA0 operator was calculated to be,

$$-\mathrm{log}{A}_{0}=0.6902{\mathrm{log}}_{10}\frac{{r}_{ij}}{17}+0.0318\left({r}_{ij}-17\right)+2.$$
(2)

Station-component corrections are tabulated in the Supplementary Information and vary between ± 0.5 magnitude units. Importantly, due to the small aperture of the network, Eq. (2) is only valid to a hypocentral distance of 40 km.

An additional effect due to the small aperture of our network is that many of the large magnitude earthquakes were clipped on a large proportion of stations. Thus, calculation of accurate magnitudes of large events is difficult. We tried to account for this by using the accelerometer channels recorded on some of the instruments to calculate magnitudes. These have the advantage of not saturating, but are recorded at a low sample rate (5 samples-per-second) and are less sensitive than the velocity channels. We processed these data by removing the gain and integrating twice to obtain the displacement. Due to uncertainties in accounting for the response of the acceleration channels, we opt to scale the magnitudes estimated from the acceleration channels to those from the velocity channels. We selected only those observations with signals above the noise and isolate those observations simultaneously made on both the velocity and acceleration channels. Then, assuming we can account for attenuation linearly (i.e. following \(-\mathrm{log}{A}_{0}=mr+c\)), we inverted for the parameters \(m\), \(c\) and individual station corrections. The magnitude can then be calculated from the acceleration recordings using the equation.

$${ML}_{i}^{acc}={\mathrm{log}}_{10}{A}_{ij}+0.0959{r}_{ij}-0.6520+{C}_{jk}^{acc}$$
(3)

Cross-correlation analysis

Differential travel-times between earthquake pairs recorded on common stations were measured using cross-correlation (Poupinet et al. 1984). We first pre-processed the waveforms by detrending, tapering and bandpass filtering between 2 and 16 Hz. Snippets of data were cut from the traces between 0.26 s before the arrival time and 1.5 s after. The traces were then cross-correlated with each other to calculate the maximum normalised cross-correlation coefficient and associated lagtime. Interpolation was used to increase the precision of the lagtime observations. We inputted these differential travel-time measurements, correlation coefficients and the QM catalogue locations into GrowClust (Trugman and Shearer 2017), a program to cluster and then relocate seismicity using relative travel times and initial reference locations, using the default clustering and weighting parameters.

Cross-correlating earthquake waveforms recorded on common stations also allows us to extract information on waveform similarity between all pairs of events. Waveforms are more similar for a given earthquake pair if the earthquakes occurred close to one another (i.e. have similar raypaths) and have a similar mechanism. Similar earthquakes can then be grouped into families using hierarchical clustering in which the normalised correlation coefficient (\(CC\)) (or more accurately \(1-CC\)) is used as the ‘distance’ between earthquake pairs recorded on a given station. Thus, a family of events is defined where the \(CC\) was higher than a user-defined threshold (0.75). By following the logic above, all the events within a family must be located close to each other (e.g. ~ 600 m for an S-wave velocity of 3.37 km/s and a frequency of 2.5 Hz) and have a similar mechanism.

To calculate the earthquake families, we use the station 7E.LAT because it had data recorded throughout the study period. Unfortunately, closer stations such as 8F.MERA or 7E.FAF had gaps where no data was recorded and other stations which recorded continuously did not detect many of these small earthquakes. We used the horizontal components aligned on the S-wave arrival because the signal-to-noise ratio of the S-wave is higher than the P-wave for these events. Similar observations have been made at other swarms analysed in Iceland (Greenfield and White 2015).

Real-time seismic amplitude measurement

Real-time seismic amplitude measurement (RSAM) is calculated at the station 8F.MERA (see Fig. 1b for location) because it is the closest station to the eruption site with a continuous seismic record. We first decimated the recordings on the two horizontal components from 250 to 50 samples-per-second (sps) before detrending and applying highpass filter with a corner at 1.5 Hz. The magnitude of the two horizontal components was calculated and the data was smoothed using a rolling mean (10 min long window). Finally, for storage and plotting the RSAM trace was resampled at 3-h intervals.

Results

The automatic earthquake catalogue generated using QM is presented in Fig. 4 and consists of 127,899 earthquakes ranging in magnitude from − 1 to 5.8. More than 70% of the earthquakes occur during four swarms where the earthquake rate exceeds 1000 earthquakes per day (Fig. 4e). These swarms occur in July 2020, August 2020, October 2020 and February–March 2021. The first three of these swarms are associated with the activation of N-S and ENE-WSW strike-slip faults across Fagradalsfjall. The large and prolonged swarm in February–March 2021 is due to the intrusion and propagation of a shallow dyke with associated triggering of seismicity in the surrounding regions due to the change in stress. The dyke first propagated north-east before propagating south-west from its initiation point (yellow star, Fig. 5), eventually reaching 10 km in length. The eruption then occurred 2 km southwest of this initiation point (orange circles, Fig. 5). The observed seismicity dropped shortly before the eruption, presumably due to the release in stress as magma accumulated near the surface. In this study, we focus on seismicity deeper than the prominent brittle-ductile transition imaged at 6.5–7.0 km depth beneath Fagradalsfjall (Fig. 4a). Upon visual inspection, all the deeper earthquakes had a lower frequency content than the shallow volcano-tectonic (VT) events (Fig. 6). Therefore, we denote these earthquakes as DLPs due their similarity to other events seen in Iceland and around the world.

Fig. 4
figure 4

Earthquake catalogue from QuakeMigrate. Earthquake hypocentres are indicated by the red dots with events identified as being located beneath the brittle-ductile transition indicated in blue. (a) Shows a cross-section from west to east with a histogram of event depth located on the left of the panel. The brittle-ductile transition is indicated by the black dashed line. A map view is shown in (c). Histograms showing the magnitude and reported horizontal location uncertainty are shown in (b) and (d) respectively. (e) Shows a histogram of the number of events per day (note the change in the y-scale). The black line indicates the cumulative histogram of all events while the vertical red line indicates the start of the eruption in March 2021. The blue line is a cumulative histogram of earthquakes deeper than the brittle-ductile transition. No vertical exaggeration

Full size image
Fig. 5
figure 5

Deep seismicity located beneath Fagradalsfjall. (a) and (b) show the latitude and depth variation through time respectively. Earthquake locations and time from QM (grey circles) and GrowClust (GC) (coloured circles) are shown together with a histogram of the number of earthquakes per day indicated by the black bars. Real-time seismic amplitude measurement (RSAM) calculated at station 8F.MERA is indicated by the red line while the onset of the eruption is indicated by the vertical orange line. Phases of the eruption (Pedersen et al. 2022) are delineated by the blue dotted lines and labelled with roman numerals. Panels (e) and (f) show cross-sections along the east–west extent of the area shown in the map views (c) and (d) respectively. The background grey circles in (c), (d), (e) and (f) show the original QM locations. The coloured circles in (c) and (e) indicate earthquakes relocated after manual picking and processing with NonLinLoc. In (d) and (f), the coloured circles indicate the GC relocations after waveform cross-correlation. The orange circles indicate the Fagradalsfjall eruption vents and the yellow star indicates the approximate location the dyke originates from. Representative horizontal and vertical uncertainties for the plotted earthquakes are indicated by the errorbars in (e) and (f)

Full size image
Fig. 6
figure 6

Frequency content of a deep (left column) and shallow (middle column) earthquake. The magnitude and hypocentral distance of the two earthquakes are similar. The north–south, east–west and vertical components recorded at station 8F.MERA and filtered with a 2 Hz highpass filter are shown in the top, middle and bottom panels respectively. Waveforms are aligned on the P-wave arrival time. Pick times and windows used to calculate the frequency are indicated by the vertical red lines and grey boxes respectively. The frequency spectra for the deep earthquake (blue) and the shallow earthquake (black) are shown on the right panels

Full size image

Location of DLPs

The DLPs are located approximately 2 km north-east of the Fagradalsfjall eruption site (orange circles, Fig. 5). Initial locations using QM have a location uncertainty of 1.5–3 km horizontally and 2–4 km vertically and indicate a sphere of seismicity centred at a depth of 10 km and with a diameter of 4 km (blue circles, Fig. 4c). Manual repicking and relocation of a subset of the higher signal-to-noise ratio events using NonLinLoc (Lomax et al. 2009) revealed a tighter cluster of seismicity (coloured circles, Fig. 5c), centred in the same location but with reported uncertainties of between 0.5 and 1.5 km horizontally and between 1 and 2 km in depth.

The double-difference relative relocation method (GrowClust, Trugman and Shearer 2017) reduces the reported relative location uncertainty to < 500 m both horizontally and vertically, revealing finer scale structure (coloured circles, Fig. 5d, f) than that visible from the automatic (grey circles, Fig. 5) or manually refined locations (coloured circles, Fig. 5c, e). Six seismogenic clusters are identified from regions with a higher density of earthquakes (Fig. 7). Clusters a–c are primarily active before the dyke intrusion, while clusters e and f are active after the dyke intrusion and during the eruptive period (Fig. 7). The single southern cluster (cluster d, Fig. 7), with a tight spatial distribution of earthquakes, was active only on December 19, 2020, during a particularly intense swarm of deep seismicity. This 'December 2020 Swarm' is further analysed in a section below.

Fig. 7
figure 7

Detailed image of the QM seismicity relocated using GC (red circles). (a) and (c) Show located seismicity through time with respect to longitude and latitude respectively. Onset of the eruption and the phases from Pederson et al. (2022) are indicated by the orange and blue dashed lines respectively. On (b), the relocated earthquakes are plotted on top of the number of earthquakes in 0.2 km bins, emphasising the clustering of earthquakes in some regions. Clusters identified by high-density regions are labelled a–f. Eruptive vent locations are indicated by the orange circles and the approximate location from which the dyke originates is marked by the yellow star

Full size image

Noticeably, there is a change in the centroid location of the DLPs with time, as events occurring early in our observation period (June–Dec 2020) located ~ 1 km to the NE and ~ 1 km deeper than later events (March–August 2021). The final reported errors of the relatively located seismicity are less than 500 m, so we can be confident that the change in the DLP centroid location is a robust observation.

DLP event magnitudes and distribution through time

We detect 747 earthquakes in the deep group of seismicity after removing false positives. The magnitude range of the deep seismicity is more restricted than the shallower activity with magnitudes ranging from − 0.4 to 1.5. Analysis of the Gutenberg-Richter distribution (Fig. 8) shows that we obtain a magnitude of completeness of 0.5 and a b value of 1.99 ± 0.13, similar to that observed in other volcanic areas (McNutt 2005). Due to the emergent and low amplitude P-wave arrivals, we could not calculate focal mechanisms for the deep earthquakes.

Fig. 8
figure 8

Gutenberg-Richter distribution for the deep earthquakes. Red triangles show the number of earthquakes in each magnitude bin (0.05 magnitude units) and the black crosses show the number of earthquakes with a magnitude greater than M (following the equation \({\mathrm{log}}_{10}N=a-bM\) where \(a\) and \(b\) are constants). The fit to the distribution is indicated by the blue line with the magnitude of completeness (Mc) indicated by the blue inverted triangle

Full size image

Through time, the deeper events occur in a series of swarms, each lasting between 1 and 4 months. Within each swarm, the event rate can reach a maximum of 18 events per day (during July 2021), but it is more usual for only 1–2 events to occur in a day (Fig. 5a, b). A noticeable feature of the data is the cessation of the deep seismicity when a dyke intrusion starts on 24 February, 2021 (Fig. 5) in the upper crust. This is coincident with a large increase in the seismicity rate, as indicated by the RSAM (Fig. 5b). Thus, it is possible that deep seismicity continues throughout the dyke intrusion phase (24 February–19 March) but that the high VT event rate masks any of the small-magnitude deep earthquakes. From our detection of two distinct swarms after the onset of the eruption phase (19 March onwards), it is clear that the noise from the tremor does not obscure the deeper earthquakes.

We identify 28 families with more than five events from the cross-correlation analysis. These families can be associated with, but not directly linked to, the clusters shown in Fig. 7. Both datasets come from the cross-correlation of events detected by QM. However, due to the location uncertainty of the GC relocations, the low-frequency content of the DLPs and because we use a single-station approach to identify the families, clusters and families do not directly overlap with one another (except for one case, see the 'December 2020 swarm' section below). The largest two families (Families 1 and 2, Fig. 9), representing ~ 22% of the catalogue, had events in them which occur throughout the observation period. Other families were active for shorter periods of time and tend to either be active before the dyke intrusion (e.g. Families 4 and 6, Fig. 9) or after (e.g. Family 5, Fig. 9). Similar results are observed from the relative relocation of the seismicity (Fig. 7). Here, clusters of seismicity tend to be active either before or after the dyke intrusion. However, some clusters and families contain earthquakes throughout the time-period. Such observations indicate that the seismogenic region and driving mechanism causing DLPs must be non-destructive and changing little throughout the observation interval.

Fig. 9
figure 9

Earthquake families detected at station 7E.LAT. The left column shows the occurrence of each family (row) through time and the right column shows the aligned waveforms (grey) and the stacked waveform (thick, black line)

Full size image

DLP event frequency content

The frequency content of the deeper seismicity is noticeably lower than that of the shallower seismicity (Fig. 6). The shallow seismicity has a broad range of frequencies present with significant energy from 2.5 to 15 Hz. However, the deeper seismicity has a more restricted range of frequencies present with energy from 1 to 5 Hz. Such a change in the frequency content is unlikely to be due to distance-dependent attenuation because despite the different depths of the two events plotted in Fig. 6, hypocentral distances are similar. This must mean that either (i) attenuation is very high in the lower crust or immediately around the deep seismicity, (ii) that the deeper seismicity has a lower stress drop or (iii) that the rupture velocity of the deep seismicity is lower.

December 2020 swarm

On 19 December 2020, a swarm occurs in which all the earthquakes locate to the same position (within error) (cluster d, Fig. 7) and share similar waveforms (Family 3, Fig. 9). Close inspection of the output of QM and the observed waveforms reveal that many small events are present below the detection threshold we used. To detect these small and low signal-to-noise ratio events we used a match-filter approach (Shelly and Hill 2011; Greenfield and White 2015; Chamberlain et al. 2017). We isolated a random event from those detected by QM (origin time, 2020–12-19T13:25:13) and generated the template by cutting a 3-s-long window around the S-wave arrival pick on the north–south component of station 8F.MERA. This station was chosen because it is the closest station to the deep earthquakes and has the highest signal-to-noise ratio. To reduce computational cost, we decimated the originally recorded 250 sps data to 50 sps and pre-processed the template by detrending and bandpass filtering between 2 and 8 Hz. The same pre-processing steps are applied to the continuous data recorded on the north–south component of station 8F.MERA between 14 and 23 December 2020. The template is then cross correlated with the continuous data.

We chose to detect an earthquake when a peak in the cross-correlation function had a value greater than 0.8. By using such a high threshold, we can assume that because the event waveforms are so similar the earthquakes must be at the same location and have the same mechanism. We can then calculate the magnitude of each newly detected event by comparing the amplitude of the recorded waveform with that of the template waveform.

The result of this procedure is the detection of 2344 extra events (Fig. 10), two orders of magnitude more than in our original catalogue over the same time interval. Earthquakes within this swarm first appeared at 2020–12-17T02:00 and continue till 2020–12-21T22:00. At its most intense, 54 events were detected in an hour and events occurred every 2 s at their peak. The magnitude of the earthquakes in the swarm range between − 1.2 and − 0.1 and peak near the middle of the swarm (Fig. 10a). Using the distribution of interevent times through the swarm, the swarm can be split into a number of sub-swarms during which the magnitude and the interevent time are correlated (Fig. 10c, d). Sub-swarms were active for 20–35 min and started with an intense period of seismicity where the interevent time is less than 100 s and the magnitudes tend to be lower. The interevent time and magnitude then increase systematically reaching a plateau between 250 and 400 s and a magnitude of − 0.2 to − 0.1. The next swarm often started immediately afterwards.

Fig. 10
figure 10

Evolution of a deep earthquake swarm beneath Fagradalsfjall in December 2020. The variation in magnitude and interevent time through the duration of the swarm are plotted in panels (a) and (b) respectively. Event times are coloured by the normalised correlation coefficient to the template event. The lower panels show a zoom in of the swarm between 2020–12-19T10:00 and 2020–12-19T13:00. The magnitude is plotted in (c), the interevent time in (d), the waveform recorded on the north channel of station 8F.MERA in (e) and the correlation coefficient in (f). In (e) and (f), the time of the triggered events (those above a normalised correlation coefficient of 0.8, red dashed line (f)) are indicated by the red circles

Full size image

Discussion

The deep seismicity we observe beneath Fagradalsfjall is located in the mid-crust (10–12 km depth), beneath the site from which magma was intruded laterally NE and SW within the upper crust. Individual seismic events in the mid-crust lack the high-frequency content found in the upper crustal VT seismicity and have frequencies (1–5 Hz) similar to those of DLPs (Chouet and Matoza 2013). In common with other DLP events, P-waves have low amplitudes and are often not observed at all, especially for the smallest events. Relative relocation of the cluster of seismicity reveals it to be located in a sphere with a diameter of 2 km. Except for the intense DLP swarm in December 2020, before the dyke intrusion on 24 February 2021, DLP seismicity is located to the north-east in at least three different clusters (clusters a–c, Fig. 7). During the eruption, the location of the DLPs shifted approximately 1 km south-west and 500 m shallower in two, more productive clusters (Fig. 5d, f; clusters e and f, Fig. 7). The December 2020 swarm of DLPs was located 1.5 km south-east of the seismically active volume (cluster d, Fig. 7).

Deep long-periods earthquakes are usually interpreted to be caused by the transport of magma through the lower crust. Beneath the Klyuchevskoy Volcanic Group on the Kamchatka Peninsula, DLPs and tremor are observed to be linked to eruptive episodes (Shapiro et al. 2017; Journeau et al. 2022). LP events (both shallow and deep) are interpreted to be triggered from pressure transients through an extended magmatic plumbing system beneath the group of volcanoes. Peaks in DLP seismicity are followed by shallow LP events and sometimes eruption. Similar links between volcanic activity and DLPs have been observed beneath volcanoes in Japan (Nakamichi et al. 2003; Kurihara and Obara 2021). Alternative models of DLP activity focus on the exsolution and rapid growth of volatile bubbles as magma rises through the crust (Melnik et al. 2020).

Some dormant volcanoes, where the active transport of magma is deemed unlikely, have also had DLPs observed beneath them. Beneath Mauna Kea (Hawaii), more than two decades of highly repetitive DLPs have been interpreted as being caused by fluids, sourced from secondary boiling of a deep magma body, moving through the crust (Wech et al. 2020). In Japan, Aso and Tsai (2014) suggest a model for the generation of DLP events in which high strain rates can be caused by the stagnation and cooling of magma at the Moho. Others cite the presence of DLP seismicity as indicating either the presence of magma within, or reactivation of a previously dormant volcano (Hensch et al. 2019). In each case, DLPs are seen to reflect fluid transport through the lower-to-mid crust, although no consensus has yet been made as to whether fluids are melt, volatiles or both (Hotovec-Ellis et al. 2018).

In Iceland, DLPs have been observed beneath Askja (Soosalu et al. 2010), Bárðarbunga (Hudson et al. 2017) and possibly Torfajökull (Soosalu et al. 2006). Other volcanoes have not had any DLPs detected despite the deployment of dense seismic networks (Krafla, Schuler et al. 2015). The best imaged DLPs in Iceland are located beneath the volcano Askja, where a long-term, dense deployment of seismometers (Greenfield et al. 2020b) has revealed at least four DLP clusters beneath Askja itself and two beneath adjacent interglacial shield volcanoes. The source of lower crustal earthquakes has been attributed to the transport of magma through the crust, although, other mechanisms such as thermal cooling and secondary boiling of stagnant magma bodies in the crust could equally well be used to explain the seismicity. Deep earthquakes have also been directly linked to lower crustal intrusions and magma flow in Iceland (Jakobsdóttir et al. 2008). A deep earthquake swarm beneath the mountain Upptyppingar, 20 km east of Askja, was associated with uplift modelled by the intrusion of a 1-m-wide inclined dyke (Hooper et al. 2011; White et al. 2011). Importantly, these earthquakes contained the high-frequencies common to upper crustal VT earthquakes rather than the DLP seismicity we observe beneath Fagradalsfjall and Askja. These observations indicate that the low-frequency content of DLP seismicity is not due solely to the temperature and rheology of the ductile lower crust in Iceland. Instead, the source mechanism of the earthquakes must be the cause.

The location and timing of the DLP seismicity we observe beneath Fagradalsfjall suggest that the seismicity is linked to the dyke intrusion and eruption in March 2021. A tempting interpretation is that the location of the DLPs represents the storage area of the erupted lava. However, early barometry results from petrological analysis of the lavas erupted at Fagradalsfjall indicate that the magma was stored at depths of 15 km or more, close to the base of the crust; before rapid transport to shallower levels and eruption (Halldórsson et al. 2022). Nevertheless, it is unlikely that the cluster of deep seismicity is completely disconnected from the transport of magma from the base of the crust to the surface. The nature of the deep seismicity changes after the onset of the eruption, with two intense seismic swarms occurring south-west of clusters active before the eruption. While the cause of this change may well remain unknown, the timing suggests that the deep seismicity does interact with the flow of magma from the source region at the base of the crust to the surface.

Given the lack of evidence for any long-term storage of magma at the location of the deep earthquakes, it would seem unlikely that secondary boiling or thermal contraction could be the cause of the deep seismicity, although we cannot discount other mechanisms which cause high pore-fluid pressure or transport of volatile-rich fluids. From the intense swarm we observe in December 2020, we know that the mechanism which generated the DLP seismicity was non-destructive and capable of generating earthquakes as quickly as two seconds apart. In the shallow crust, repeating LP events are sometimes called ‘drumbeat’ events and can be associated with dome growth in andesitic or dacitic volcanoes (Iverson et al. 2006), gas flux and resonance through a shallow crack network (Chouet 1996; Bell et al. 2017) or brittle failure at the edges of a magma plug in a vent (Neuberg et al. 2000). While DLPs should not be referred to as drumbeat events because they can occur during non-eruptive phases and at much greater depth than shallow repeating LP sequences, there are some similarities between the swarm we observe beneath Fagradalsfjall and drumbeat events. In particular, variations in the interevent time were observed beneath Mount St. Helens (Matoza and Chouet 2010) and the positive correlation between interevent time and magnitude (Fig. 10) is similar to that observed beneath Tenerife (D’Auria et al. 2019). However, it is noticeable that in both the Mount St. Helens and Tenerife case, variations in interevent time occur on longer timescales (hours to days) than the 30-min-long subswarms observed beneath Fagradalsfjall.

A DLP earthquake swarm beneath the interglacial shield volcano Vaðalda, in central Iceland, shows similar properties to those we observe beneath Fagradalsfjall (Greenfield and White 2015). Here, inter-event time and magnitude were correlated during a series of swarms in December 2012 at depths of 20–25 km. To produce such similar, repeating earthquakes and the observed non-double-couple source mechanisms, Greenfield and White (2015) proposed a model in which failure is caused along faults at the tips of sills in the mid-crust. Intrusion of CO2-rich fluids along these fault planes increases the pore fluid pressure and causes brittle failure in the otherwise aseismic crust. The almost constant moment release rate observed for some swarms implies a constant driving stress with the size of each earthquake being governed by the extent of the fault plane exhibiting stick–slip failure.

Beneath Fagradalsfjall, a similar mechanism could be invoked, although we are conscious that we have fewer constraints on the source mechanism than beneath Vaðalda. In particular, it has not been possible to generate any focal mechanisms due to poor focal sphere coverage. Other models such as volatile flow through the fracture network (D’Auria et al. 2019) or active transport of the melt itself could equally well generate the DLP seismicity we observe. We do know that magma is actively migrating upwards from the melt storage region at depths greater than 15 km to the upper crust (Halldórsson et al. 2022). However, it is not clear why such transport should only trigger seismicity at a particular depth in the mid-crust rather than at a range of depths through the crust.

Our observations suggest a model where at 10–12 km depth there is a barrier to upward transport of magma (Fig. 11a). Following the intrusion of a fresh batch of melt into the near-Moho magma storage region, the magma was driven upwards through the crust. The upwards transport was aseismic, in common with much of the magma transport through the Icelandic lower crust (Tarasewicz et al. 2012; Sigmundsson et al. 2015). At 10–12 km depth, a barrier, possibly an olivine-rich cumulate or other rheological change prevents or slows further magma intrusion (Taisne et al. 2011). Lower crustal reflectors observed at a 10–12 km depth but offshore the Reykjanes Peninsula could represent similar features (Weir et al. 2001). Seismicity at these depths could be induced by the intrusion of magma near the barrier at higher strain rates than that during the earlier magma transport. Alternatively, CO2-rich fluids generated by degassing of the magma could be concentrated at this depth and induce seismicity due to either bubbles within the magma (Melnik et al. 2020) or mass transport through a fracture system (Wech et al. 2020).

Fig. 11
figure 11

Schematic summary plot showing our interpretation for the mechanism that could be inducing seismicity in the lower crust. Before the onset of the eruption, DLP seismicity is triggered at mid-crustal depths at a barrier to further transport (a). During the dyke intrusion and phases I and II of the eruption, a continuous magma channel to the surface means no seismicity is triggered (b). Changes in the magma flux through the magma plumbing system causes seismicity to be triggered in a new location, 1 km SW of the original clusters, and changes to the eruption rate at the surface (c). The large magma body in the upper crust is a representation of the dyke intrusion which, in reality, would be approximately 1 m wide

Full size image

Carbon dioxide fluxes through the Icelandic crust induce DLP seismicity

The pore-fluid pressure increase associated with the exsolution of CO2 has been used to explain the occurrence of DLP seismicity in Iceland (Greenfield and White 2015; White et al. 2019) and other deep seismic swarms around the world (Shelly and Hill 2011; Hotovec-Ellis et al. 2018). Under Iceland, CO2 is likely to exsolve at depths of between 5 and 25 km (Pan et al. 1991; Lowenstern 2001) and because of the low viscosity and high temperature of basaltic melts, the volatile phase can become decoupled from the source magma (White et al. 2019) and thus ascend through the crust.

On the Reykjanes Peninsula, a series of uplift and subsidence events under the Svartsengi high-temperature geothermal field located 8 km west of Fagradalsfjall (Fig. 1a) during the first half of 2020 has been modelled as cyclic intrusions of CO2-rich fluids into a permeable aquifer below the geothermal field at 4 km depth (Flóvenz et al. 2022). The fluids were hypothesized to be transported vertically from a deep magma storage region beneath Fagradalsfjall to the brittle ductile boundary. After this, three pulses of fluid were transported west and one east to the Svartsengi and Krýsuvík high-temperature geothermal fields respectively. The model we propose does not exclude CO2 pulses through the crust as driving the DLP seismicity we observe. However, the lack of any shallow LP or migrating VT events in our earthquake catalogue would indicate that any transport of CO2-rich fluids along the brittle-ductile transition would have been aseismic. Nevertheless, release of CO2-rich fluids from lower-crustal magma storage bodies has caused seismicity at the brittle-ductile transition in central Iceland (Martens and White 2013) and throughout the lower crust beneath Mammoth Mountain (Hotovec-Ellis et al. 2018).

A clear indication of whether the CO2-rich fluids that intruded beneath Svartsengi or Krýsuvík were sourced beneath Fagradalsfjall and drive the DLP seismicity would be to see when the DLPs first occur and whether that correlates with the first intrusion event during January 2020 (Flóvenz et al. 2022). However, our current dataset only begins in June 2020, so such comparisons cannot be made. We note though that CO2-rich fluids could be sourced from the presumably more mature magma plumbing systems directly beneath Svartsengi and Krýsuvík. The stresses associated with a rifting episode along the entire plate boundary could be causing the synchronous activity observed across the multiple volcanic systems on the Peninsula (Cubuk-Sabuncu et al. 2021; Geirsson et al. 2021).

Linking DLPs to magma ascent

At some point, the barrier must have been breached and magma found its way to intrude to shallower depths in the crust (Fig. 11b). There is a reduction in DLP event rate in November 2020 which speculatively could be related to the breach of the barrier, but there is little evidence in the seismicity or petrological data to argue conclusively for this date and the intense swarm in December 2020 occurred after this date. DLPs continued until the onset of the dyke intrusion in February 2021, but no DLPs are detected during the dyke intrusion itself. The high number of shallow VT events associated with the dyke intrusion (and therefore the higher noise levels) could be the reason for no DLPs being detected, but the lack of DLPs immediately post-eruption does suggest that, initially, the release of stress due to the eruption changed the conditions at depth such that triggering of DLP events was inhibited. A reason could be that the near-Moho magma storage region was connected to the surface via a melt channel through which magma could flow freely. Diffusion chronometry provides direct evidence for timescales of melt transport from the Moho to surface in days; this is a maximum bound, and physical models of flow in open conduits can transport melt across the crust in just a few hours (Mutch et al. 2019b). Importantly, the barometry of melt storage depths based on liquid compositions provides evidence for equilibration at near-Moho depths and that some melts also stalled and re-equilibrated in the shallow crust (Halldórsson et al. 2022). The petrological evidence indicates that storage and equilibration at ~ 10 km depth was not an important process under Fagradalsfjall. Thus, any magmas intruded at a depth of 10–12 km must be a small volume or, could be present in the magma which forms the cooled dyke.

Other than two scattered DLPs at the beginning of April 2021, DLP seismicity returned in May 2021, located approximately 1 km south-west of and shallower than the pre-eruption clusters (Figs. 5, 7 and 11). During this period, the eruption rate doubled from 6 to 12 m3/s (Phase III, Figs. 5 and 7), presumably due to increased flux from the mantle (Pedersen et al. 2022). The increased magma pressure from below or the development of new magma channels through the lower crust could explain the re-emergence of DLP seismicity. Another DLP swarm occurred at the end of June 2021 when the eruption transitioned from short-term episodes (< 15 min duration) to a phase of longer duration episodes (9–26 h) separated by equally long repose intervals (Phase IV, Figs. 5 and 7) (Eibl et al. 2022). The eruption rate fell at this point, implying some change in the flux of magma into the near-Moho magma storage region. It may seem unusual for magmatic pressure to increase or additional intrusions to develop while connected to the surface. However, recent eruptions at Holuhraun (Sigmundsson et al. 2015; Eibl et al. 2017; Woods et al. 2019), Kīlauea (Owen et al. 2000; Neal et al. 2019) and Krafla (Harris et al. 2000) show that magma can intrude into other regions despite being connected to the surface by active conduits.

Conclusions

Between June 2020 and August 2021, we have detected a cluster of deep long period (DLP) earthquakes at 10–12 km depth in the mid-crust beneath Fagradalsfjall on the Reykjanes Peninsula, Iceland. The low frequency content of the detected seismicity was similar to other DLP events observed beneath volcanoes in Iceland (e.g., Askja and Bárðarbunga) and around the world (e.g., Mauna Kea and Klyuchevskoy). The epicentre of the DLP seismicity was ~ 2 km north-east of the location of the Fagradalsfjall eruption. However, when comparing the depth of the DLP cluster to the storage depths calculated from petrological modelling, it is clear that the DLP earthquakes cannot be due to active magma transport close to the magma storage region. Instead, we suggest that these DLP earthquakes were triggered close to a barrier to upward transport of magma in the lower crust either by CO2-rich fluids being transported through a fracture network or high strain-rates caused by magma intrusion. DLP seismicity stopped at the onset of the dyke intrusion but reappeared in a new location close to periods when the discharge rate and fountaining style of the eruption changed. Such observations indicate a link between DLP earthquake activity, changes to the magma plumbing system and supply rate. Real-time monitoring of DLP earthquakes, while small and difficult to detect, could therefore be an important tool in detecting changes in the magma pathway that could affect eruption style when direct observations are difficult.