Temporal Seismic Velocity Changes During the 2020 Rapid Inflation at Mt. Þorbjörn‐Svartsengi, Iceland, Using Seismic Ambient Noise

Repeated periods of inflation‐deflation in the vicinity of Mt. Þorbjörn‐Svartsengi, SW‐Iceland, were detected in January–July, 2020. We used seismic ambient noise and interferometry to characterize temporal variations of seismic velocities (dv/v, %). This is the first time in Iceland that dv/v variations are monitored in near real‐time during volcanic unrest. The seismic station closest to the inflation source center (∼1 km) showed the largest velocity drop (∼1%). Different frequency range measurements, from 0.1 to 2 Hz, show dv/v variations, which we interpret in terms of varying depth sensitivity. The dv/v correlates with deformation measurements (GPS, InSAR), over the unrest period, indicating sensitivity to similar crustal processes. We interpret the velocity drop to be caused by crack opening triggered by intrusive magmatic activity. We conclude that single‐station cross‐component analyses provide the most robust solutions for early detection of magmatic activity.

additional intrusions occurred between March 6-April 17 and May 15-July 22, 2020 based on deformation observations. At the time of revising this report, intense seismic activity and deformation is observed in Fagradalsfjall (a part of KVS), which resulted in an effusive eruption in Geldingadalir valley on March 19, 2021 (see Figure 1a).
Monitoring the evolution of this volcanic unrest is crucial for effective early warning of volcanic hazards, predominantly lava flows (Tarquini et al., 2020), volcanic gas, rifting, and earthquakes. The center of the Mt. Þorbjörn-Svartsengi uplift is only 2 km from the town of Grindavík (population ∼3,300) and the Svartsengi Power plant, providing electricity, and cold and hot water; and a few tens of kilometers away from the national airport of Keflavík and main roads, posing a significant risk to the society (Figures 1a and 1b). Traditional volcano monitoring programs benefit from seismic networks, ground and satellite-based deformation CUBUK-SABUNCU ET AL. . Events related with the first, second, and third intrusion periods are shown with green, red and pink circles respectively. Gray and purple colors indicate earthquakes that occurred before January 20 and after July 23, 2020, respectively. Black triangles mark seismic stations with an analysis period spanning from February to November 2020, while blue triangles indicate permanent seismic stations. The green square shows the location of the cGPS station SKSH. Black circles represent important infrastructure locations, and yellow stars denote M ≥ 4.8 earthquakes. The purple box shows the area affected by intrusions. The seismicity compared to the dv/v curve of GRV is confined to the black box. A light blue square marks the location of the Geldingadalir volcanic eruption, (b) InSAR map showing deformation in the satellite's line-ofsight (LOS) from the Sentinel-1 PS-InSAR processing for T16. The red colors indicate movement toward the satellite (inflation). observations, and geochemical measurements. Advances in seismic noise-based methods have provided new data processing approaches (e.g., Brenguier et al., 2008) and real-time monitoring tools (Lecocq et al., 2014). Considerable research has been carried out in the field of ambient noise seismic interferometry (ANSI), including volcano monitoring (e.g., Brenguier et al., 2016). Previous studies of Jónsdóttir (2018) and Donaldson et al. (2019) have provided essential constraints on the variations of seismic velocities in Iceland associated with the 2011 unrest at Katla and the 2014 Bárðarbunga-Holuhraun intrusion, respectively.
We focus on the recent unrest (January-July 2020) in Mt. Þorbjörn-Svartsengi, which offers a tangible case study for examining seismic velocity changes over time in response to a series of repeated magmatic intrusions in the RP. Our primary objective is to use ANSI in near-real-time to characterize temporal velocity changes to track the physical changes in the medium related to the unrest.

Data Analysis
We analyzed three-component continuous waveforms of five broadband seismic stations, a part of the permanent seismic network (Bodvarsson et al., 1996) operated by the Icelandic Meteorological Office (IMO). We also used additional data from three seismic stations (ISS, LAG, and LSF) provided by the Czech Academy of Sciences and two stations co-operated with Iceland Geosurvey (LFE, RAH) ( Figure 1a). The data set spans from April 2018 to November 2020 ( Figure S1). We performed a visual quality assessment on data to exclude large gaps and clock errors and selected a period (August-November 2019) to construct a reference without significant volcanic/tectonic activity. All seismological data are processed using the MSNoise software package. Our approach is similar to the methodology described in Lecocq et al. (2014) and references therein. We analyzed noise cross-correlations at four frequency bands (f 1 :0.1-1.0 Hz, f 2 :0.5-1 Hz, f 3 : 0.4-0.8 Hz, and f 4 :1-2 Hz) to compute seismic velocity changes ( Figure S1).
The Persistent Scatterers (PS) InSAR technique (Hooper, 2008;Hooper et al., 2004) was used to generate line-of-sight (LOS) deformation maps (Figures 1b and S4) covering the RP using satellite data acquired from Sentinel-1 ascending track 16 and descending track 155. Time series of LOS deformation were extracted close to the center of the deformation signal (at cGPS station SKSH) resulting from the magmatic intrusions ( Figure 2). We used the vertical component GPS time series of the SKSH station (Geirsson et al., 2010), processed using the GIPSY/OASIS II software (Zumberge et al., 1997). The GPS/InSAR observation period is chosen to be long enough to compare with the dv/v time series. The reference frame used for SKSH is based on the ITRF-2008. The InSAR, daily GPS and dv/v time series are smoothed with a rolling mean of 3 days.

Cross-Correlations and dv/v Measurements
We computed the relative seismic velocity changes taking place in the RP as a function of time, using ANSI. Figure 2 summarizes the temporal evolution of dv/v for the seismic station GRV, which showed the largest seismic velocity change during the unrest. The three-component-based correlation functions were computed for single stations (NZ, EZ, and EN). We also calculated the cross-correlation functions (CCFs) between station pairs for vertical component records (ZZ, Figure S2). We compare our results with the cumulative number of earthquakes and deformation data. For a detailed description of the methodology and other results, the reader is referred to the supporting information.
Short-term subsurface velocity changes during the inflation are successfully detected using a single station (SC) cross-component analysis. The nearest seismic station (GRV) is located at the southern part of the uplifted area and detected velocity changes in f 4 :1-2 Hz frequency band. On January 21, 2020, we observed a rapid dv/v decrease in the north-vertical cross-component (NZ) of the GRV seismic station (Figure 2a). After reaching its minimum value (−1%) on February 17, relative velocity changes for the NZ cross-component progressively recovered to the pre-unrest values.
A second distinct velocity drop in the NZ cross-component was detected between early March to April, with a minimum of −0.8%, suggesting reinflation during this period. A small decrease in dv/v (∼0.2%) was observed in EZ and EN components before the third intrusion in mid-May. In the NZ cross-component, a 0.35% dv/v drop was detected coinciding with the third intrusion which terminated in June. Approximately four months after the first intrusion, seismic velocity changes returned to pre-unrest background levels ( Figure 2a), possibly indicating a pause in magmatic activity around Mt. Þorbjörn, and rapid recovery of the stress in the medium. Similar characteristics were observed between all components, but the magnitude of changes in the EZ and EN components was smaller (−0.3%) than the NZ. Differences in the amplitudes of component pairs can be attributed to the heterogeneities and seismic anisotropy (e.g., Machacca et al., 2019) that may exist in the upper crust. We suggest that the structural changes below GRV should be larger for the intrusion before March, as the amplitude of dv/v drops to much lower values (NZ, −1%) in comparison with the other intrusion periods.
In Figure 2a, the black over purple-to-green color scale shows the correlation coefficient. We observed a high similarity (correlation coefficient >0.9-1), indicating a good match between the reference and daily seismograms, before the first intrusion (Figures 2a and S5). However, lower values (∼0.7) were persistent among all components for the unrest period. Development of new scatterers and fluids (e.g., Larose et al., 2010;  . The gray boxes correspond to periods related to Mt. Þorbjörn inflation. Yellow filled boxes mark the timeline of Fagradalsfjall events. The color scale displays the correlation coefficient between the reference and daily cross-correlation functions (CCFs). Error is plotted on the dv/v curve with short gray lines. Purple arrows indicate dv/v changes that are observed before the deformation signal. (b) The catalog of earthquakes in the analysis period (see black box in Figure 1a). The red curve shows the cumulative number of events. (c) Vertical component of GPS (black) and the InSAR profiles (blue, purple) from the SKSH station. The orange crosses mark the timing when the uplift surpasses the previous uplift signal.
with time suggests the upper crust in the vicinity of Mt. Þorbjörn has been permanently deformed and damaged due to local intrusive activity (e.g., Feng et al., 2020;Viens et al., 2018).

Overview of the Seismicity
We summarize the seismic activity associated with the recent inflation in Svartsengi (Figure 1a). We use the IMO earthquake catalog to review seismicity, and standard location errors are within the order of 1 km in the horizontal direction. We only consider events with M ≥ 1. RP is a seismically active region known for intense swarm activity. Approximately ∼3,500 earthquakes with magnitudes M ≥ 1 (see S8) were detected between January and November 2020 north of Grindavík and the central part of the RP (west of −22.14°W).
We distinguished two intense seismic swarms in the Mt. Þorbjörn-Svartsengi and Fagradalsfjall (see Figure 1a). The cumulative number of earthquakes showed increased swarm activity as of January 23, and magnitude distribution further revealed that events concentrated into three groups over time (Figure 2b). Approximately 500 earthquakes (M ≥ 1) were recorded between late January-March (green circles in Figure 1a). The timing of the swarm correlates with a substantial seismic velocity drop (average dv/v, 0.5%) at the GRV station during the first intrusive period. Following 12 March, the earthquakes continued to occur in the vicinity of Mt. Þorbjörn (red circles in Figure 1a), coinciding with the second inflation period. Of the ∼700 events from mid-March to the end of April, one earthquake was recorded with a magnitude M 4.8 (March 12). The seismicity abated during May, but more than 400 earthquakes occurred across the inflating area (purple box in Figure 1a) in the third intrusion period. The earthquakes occurred at depths down to ∼6-8 km during the unrest period. Swarm activity in the region rapidly migrated to ∼10 km NE of Grindavík on July 19 (pink circles, Figure 1a). In this period, the largest recorded seismic events were two M 4.9 earthquakes in Fagradalsfjall on July 19-20. On October 20, a M 5.6 earthquake occurred near Krýsuvík, and the coseismic change due to this event has been captured by dv/v at ISS (−0.5%), LSF and KRI (−0.4%) ( Figures S1 and S7). The dramatic increase of seismicity indicates ongoing active deformation spreading in a relatively broad zone (∼25 km) between Reykjanes, Svartsengi, Fagradalsfjall-Krýsuvík. It is likely the result of the interaction of stresses caused by a magma injection and the regional crustal stresses.

Changes in dv/v Compared With GPS and InSAR
We compare the dv/v time series obtained at GRV with the deformation changes demonstrated by Geirsson et al. (2020). Their geodetic modeling results suggest that the deformation from three inflation sources is roughly collocated, centered ∼1-1.5 km west of Mt. Þorbjörn at depths between 3.2 and 4.3 km. The total LOS deformation derived from PS-InSAR analysis caused by these multiple intrusions is displayed in Figure 1b. Vertical crustal motion recorded at the SKSH cGPS station displays clear inflation signals, which correlate with the periods of seismic velocity decrease (Figures 2a-2c). SKSH station displayed ∼60 mm of uplift between January 21 and February 6, a period also accompanied by a high rate of seismicity. We observed ∼50-60 mm increment between March 6 and April 17, and clear deflation after all inflation periods. Over the observed time range, the PS-InSAR LOS patterns correlate with the GPS displacements ( Figure 2c). We identified a sudden increase (∼40 mm) in the LOS measurements corresponding to the first intrusive period. Similar ground uplift was observed during the second and third inflation periods from LOS deformation (∼30 mm).
Inflation and seismicity ceased in May, but again, a third inflation episode was picked up by deformation data (May 15), whereas seismicity accelerated in June (Figures 2b and 2c). A comparison of the deformation and seismicity suggests that earthquake activity due to the onset of a new uplift, intrusive period, is only triggered once the deformation has surpassed its previous maximum uplift (see Figure 2c, orange cross). Seismic velocities continued decreasing throughout this period (purple arrows in Figure 2a). The overall changes in seismicity, dv/v, GPS, and InSAR could be explained by the response of the crust to repeated intrusions. Dv/v and deformation observations might suggest that intrusions that occurred before April 17 contributed more to the opening of new cracks, and smaller anomalies after this date could reflect an already fractured medium (less resistant to magma migration) and/or a lower rate of magma flow. Indeed, the deformation rate was the fastest in the first inflation episode (∼3-4 mm/day), then slower during the second inflation, and even slower, but with a longer duration, during the third inflation period. Thus, the strain rate may also affect how deformation and seismic velocity variations relate (e.g., Rivet et al., 2013).

Seasonal Effects
Identifying the seasonal signals arising from diverse environmental factors (e.g., precipitation, snow loading, and temperature) and eliminating their effects from dv/v contributes to more acute monitoring of velocity changes during active volcanic periods. In Iceland, strong seasonal contrasts occur with snow accumulation and melting. Donaldson et al. (2019) reported a striking seasonal cycle of dv/v in central Iceland, interpreted as the elastic loading in the shallow crust and groundwater level (GWL) changes. According to their study, dv/v values increase through October-April and decrease during the summer and fall. The average magnitude of these annual changes was ±0.05% and ±0.1% for 0.1-0.4 and 0.4-1.0 Hz frequency bands, respectively. The dv/v changes were larger (±0.25%) for 1-2 Hz. Analysis of GPS data (Drouin et al., 2016) also revealed a similar seasonality (winter-spring and summer-fall) interpreted in terms of snow and glacial loading effects. Deformation, in connection with these large seasonality changes, is more representative of the central highlands of Iceland (27 mm) than the coastal areas (4 mm) in Iceland (Drouin et al., 2016). Our analyses revealed a seasonal cycle at the f 1 :0.1-1 Hz frequency band ( Figure S1a) with maximum changes less than ±0.2%, and even smaller at the f 4 :1-2 Hz frequency band (<±0.1%). We find that positive dv/v tends to peak in the RP during the summer ( Figure S1a) since the changes in GWL and snow load contribution to stresses are considerably less in this coastal region compared to ice-covered areas. During the study period, the magnitude of the seasonal changes at high frequency is much lower than the dv/v induced by the magmatic intrusion (−1%). We thus interpret dv/v above ±0.15% at the 1-2 Hz frequency band as the result of physical changes induced by the recent magmatic processes. Understanding the peninsula's longterm seasonal effects would be further improved by introducing more data to our analysis.

Frequency Dependence of dv/v
The lapse-time dependence of dv/v provides an estimate on the depth sensitivity of our measurements Obermann, Planès, Larose, Sens-Schönfelder, & Campillo, 2013). Figure S6 displays that the variation of dv/v (f 4 :1-2 Hz) along the coda for the GRV seismic station is nearly constant. This finding also implies that body waves contribute more to the coda.
Scattering and attenuation of seismic waves affect the energy propagation in the crust (Lawrence and Prieto, 2011). High-frequency signals are mostly sensitive to shallow changes in the medium due to strong scattering and attenuation. Low-frequency signals are less attenuated and scattered, allowing them to depict velocity changes from deeper parts of the subsurface (Aki & Chouet, 1975;Rivet et al., 2011). In the high-frequency range (f 4 :1-2 Hz), dv/v is sensitive to changes in the shallow crust near the GRV seismic station. Following each dv/v decrease, the seismic velocity rapidly recovers (Figure 2a). In contrast, permanent scatterers may have developed in the medium as indicated by the permanent drops of the cross-correlation coefficients. We also find that the dv/v drops at GRV (Figures S1a and S1b) at lower frequency ranges (f 2 :0.5-1 Hz, f 3 :0.4-0.8 Hz) are permanent. They may have been induced by magma intrusions at relatively deeper parts in the crust, consistent with the deformation modeling (Geirsson et al., 2020), which indicates intrusions at depths of 3.2-4.3 km.

Seismic Velocity Changes
Analyses in the 1-2 Hz frequency band revealed that the rapid dv/v drop is very consistent among five receivers (GRV, ISS, LFE, LSF, and LAG). We describe our findings for two station groups (see Figures 1a  and S1). The stations in Figures S1c and S1d, encompass the area controlled by intrusions (<8 km), and good coherence between them suggests that velocity changes are originating from similar sources in the crust.
The intrusion effect was large during the second uplift episode at station ISS in the east, where the dv/v dropped from 0.5% to 0%. We find that dv/v continued to decrease at stations ISS, LSF, and LAG until May, and when the second intrusion terminated, we further observed a similar drop at two other stations (LFE, GRV) during May. However, the deformation did not clearly show inflation, as it was mostly stable during the same period. This small difference in temporal characteristics between closely spaced seismic stations (4-10 km) might imply that the physical changes initially took place close to the NE part of the inflating region, and the impact of the related lateral stresses only picked up approximately a week later in western stations (e.g., LFE, GRV). We suggest that the crustal structure beneath the RAH station was not strongly modified since almost no variation in velocity was observed, and this points out less change in dv/v toward RAH in the medium.
Dv/v dramatically decreased to −0.6% (LSF, LFE- Figure S1d, yellow boxes) since swarm activity has increased in Fagradalsfjall on July 22, 2020. Simultaneously, we also observed a velocity decrease (0.2%) in the eastern station KRI ( Figure S1b, yellow boxes), possibly indicating a stress change in the Krýsuvík volcanic system. Following the M 5.6 earthquake (Figures 1a and S7), the dv/v values at KRI further dropped about 0.4% approximately a month after the first observed velocity change at this station ( Figure S1b). Hence, these overall temporal changes in dv/v might be associated with propagating volcano-tectonic events between Svartsengi, Fagradalsfjall to Krýsuvík, and possibly the dv/v velocity drop in Krýsuvík is associated with a local magmatic intrusion which started in July 2020.
For frequencies below 1 Hz, the average dv/v at GRV dropped to −0.3% only for 0.5-1 Hz and 0.4-0.8 Hz frequency bands, corresponding to the timeline of intrusions (gray boxes in Figures S1a and S1b). We also note a 0.2% decrease of dv/v for stations KRI and ISS in response to the Fagradalsfjall-Krýsuvík area´s reactivation in late July.
Different characteristics of dv/v between local stations and those far from the uplift (e.g., VOG, NYL) imply that the stress changes were limited to distances of < 8-10 km. The vertical component cross-correlation analyses ( Figure S2) for different receiver pairs showed that the average dv/v was mainly oscillating within a narrow range (dv/v %, <±0.2) without a clear velocity change. The station spacing in the RP ranges from 4 to 34 km, and the medium is characterized by strong attenuation down to 4 km (Menke et al., 1995) as a result of the heterogeneous crust (e.g., fissures, faults, geothermal activity). Thus, high attenuation in the uppermost crust likely cancels out the higher frequencies within the noise wavefield between far apart receivers. We conclude that dv/v based on individual station analyses are more reliable/relevant than cross-correlations between station pairs for this study.

Volumetric Strain Changes Under GRV Seismic Station
Several studies have emphasized the sensitivity of seismic waves to the opening/closing of cracks in the medium (e.g., Nur, 1971;Yamamura et al., 2003). The closing of fractures during magmatic inflation can explain the increase in seismic velocity (O'Connell & Budiansky, 1974). Mechanism of crack opening also assists inflation of the magmatic material (Goulty & Schofield, 2008;Johnson & Pollard, 1973) and decreases seismic velocities. Recently, Donaldson et al. (2017) carried out strain modeling for the Kílauea volcano and demonstrated how elucidating volumetric strain variations contributed to a better understanding of changes in dv/v.
To gain insight into the spatial distribution of crustal strain and crack opening around the GRV station and its relationship with dv/v, we further compared our data with the volumetric strain caused by the intrusions. We used the sill model for the first intrusive episode at Svartsengi from Geirsson et al. (2020) and used the Coulomb software (Toda et al., 2011) to compute the resulting changes in volumetric strain (sum of the directional strains in each direction, x, y, and z) in the subsurface assuming a uniform elastic halfspace. The increased dilatational strain should correspond to the opening of cracks and therefore lowering of dv/v. The strain model (Figures 3a and 3b) shows contraction directly above and below the sill and two dilatational lobes at the sides. Furthermore, a region of dilatational strain is concentrated near the surface, above ∼1.5 km. The GRV seismic station is located toward the edge of this near-surface extensional area, and we observe that large dv/v variations (−1%) are related to maximum strain changes of ∼10 microstrain at the surface and around ∼30 microstrains at 4 km depth. The strong relationship between positive strain changes (dilatation) and decreasing dv/v is very evident in the 1-2 Hz frequency band since dv/v mostly samples shallow crustal parts in this range. Further below, at 2-4 km depth, there is a sharp gradient in the strain field such that deformation model uncertainties may easily flip the computed strain sign locally at this depth below GRV. We find, however, a striking decrease in dv/v observed at the frequency bands of 0.5-1 and 0.4-0.8 Hz (Figures S1a and S1b), indicating that indeed the net volumetric stress change at deeper parts of the crust is dilatational as indicated by our model.
Recently, Takano et al. (2019) investigated the strain sensitivity in the 1-2 Hz frequency band using different measuring time windows in the coda of the CCFs (lapse-time). They indicated that dv/v respond to dilatation and contraction at early lapse-times (i.e., 2-7 s). However, the influence of strain changes on seismic velocities is less with increasing lapse-times after 15 s, and wave velocities do not change considerably with strain for late measurement times (i.e., 25 s). These observations support that our dv/v results, which have an analysis time between 5 and 25 s, are sensitive to dilatational strain variations. The volumetric strain model reveals that the extensional stresses and opening of fractures play a role in decreasing dv/v in the vicinity of Mt. Þorbjörn.

Conclusion
We present temporal dv/v variations based on cross-correlations of seismic ambient noise for monitoring the January-July, 2020 volcanic unrest period in the RP. Our study provides a multidisciplinary analysis of how velocity changes emerge due to repeated magmatic intrusions in a complex plate boundary zone. We identify co-intrusive seismic velocity reduction in the 1-2 Hz frequency band by analyzing ∼3 years of continuous waveforms. Results confirm a noticeable change in the velocity structure beneath Mt. Þorbjörn coinciding with the unrest period. Crustal magma intrusions and concurrent opening of new cracks to accommodate magma propagation may have resulted in the observed seismic velocity decrease. We find that the dilatational strain below the GRV station plays a critical role in reducing seismic wave velocities. Similar temporal history of changes pointed out by seismic and deformation data supports that velocity decrease is closely linked to new magmatic intrusions. Capturing changes in dv/v contributed to operational monitoring for authorities as a complementary approach to provide timely warnings to the public. The capability to monitor these changes in real-time would facilitate timely detection of future intrusive events/changes in volcanic activity in Iceland. For future work, we suggest modeling to improve our understanding of how magma migration, cooling and relaxation, and seismic anisotropy affect temporal changes of dv/v.