Temporal Seismic Velocity Changes Associated With the Mw 6.1, May 2008 Ölfus Doublet, South Iceland: A Joint Interpretation From dv/v and GPS

In South Iceland, populated and agricultural areas are at risk of earthquakes due to their location within the South Iceland Seismic Zone (SISZ). In 2008, two moderate‐sized earthquakes (M5.8 and M5.9) occurred in Ölfus, the western end of this highly active transform zone. We analyze temporal seismic velocity variations (dv/v) related to the Ölfus earthquake doublet, using cross‐correlations of ambient noise in the frequency range of 0.1–3.0 Hz. The two mainshocks decrease the average velocity by 0.8% at the nearest stations. The co‐seismic changes are most noticeable from 0.7 to 1.7 Hz and affect a 40 km wide region. We present a first‐time comparison of dv/v to crustal deformation, seismicity, co‐seismic volumetric stress changes and reported PGA distribution for the Ölfus doublet. Ground accelerations caused by mainshocks at intermediate distances suggest that strong shaking‐related damage may contribute to the co‐seismic dv/v decrease. A rapid velocity increase (0.3%) in a month after the co‐seismic drop indicates crustal rock healing. We find 3‐months of post‐seismic decorrelation, followed by a nearly permanent velocity decrease (0.2%) confined to a shallow layer (1 km) until the end of the observation period. Afterslip and pore fluid effects in the near‐source region are likely to influence post‐seismic dv/v. We demonstrate that seismic interferometry can contribute to future fault‐zone monitoring operations in the SISZ by detecting small changes in velocity.


Introduction
Iceland is intersected by the mid-Atlantic divergent plate boundary, separating the Eurasian Plate to the East and the North American plate to the West.The boundary across Iceland can be categorized by zones of active tectonic extensional and transform motion.Volcanic zones coincide with extensional regimes while regions of transform tectonic motion define the seismic zones: the South Iceland Seismic Zone (SISZ) in the south and the Tjörnes Fracture Zone, mostly offshore north Iceland.
The relative plate motion across the ridge is estimated to be 2 cm per year (cm/yr) in direction ∼100° (DeMets et al., 1990(DeMets et al., , 1994;;Sella et al., 2002).The SISZ extends some 75-100 km from east to west and is ∼15 km from north to south (Bjarnason, Cowie, et al., 1993).The SISZ is a left-lateral transform zone that connects to the Eastern Volcanic Zone and the Hengill triple junction, where the Reykjanes Oblique Rift and Western Volcanic Zone meet.Rather than a single left-lateral E-W transform fault, the plate motion in the SISZ is primarily accommodated by bookshelf tectonics on numerous right-lateral N-S strike-slip faults (Einarsson, 2008).
On 29 May 2008 (15:45 UTC), two earthquakes on adjacent N-S aligned faults struck the Ölfus region in the SISZ, separated by 3 s and 5 km, first on the Ingólfsfjall fault (to the east) and later on the Kross fault that runs through Hveragerdi (Hreinsdóttir et al., 2009).The estimated moment magnitudes of the two main events, based on geodetic data, were M w 5.8 and M w 5.9, respectively, which is equivalent to a total moment of M w 6.1 (Decriem et al., 2010).The 2008 Ölfus doublet is one of Iceland's largest instrumentally recorded event, and the largest one yet to be recorded together with continuous seismic data spanning over months preceding and following the event.Earthquake aftershock activity was high until 9 June 2008, and outlined both N-S oriented and E-W trending faults (Brandsdóttir et al., 2010).Deformation modeling (Decriem et al., 2010) suggested that the area of highest slip on the Ingólfsfjall fault was between 2 and 4 km depth, while on the Kross fault it was between 3 and 6 km depth.The doublet caused high intensity ground motion (Halldórsson & Sigbjörnsson, 2009), also damage in nearby towns and was the single most expensive natural catastrophe in Iceland since the Heimaey (Vestmannaeyjar archipelago) eruption in 1973 (H.Árnadóttir et al., 2014).
The relationship between variations in seismic wave velocities and moderate earthquakes in Iceland is still largely unknown.We use ambient noise-based seismic interferometry, which is a powerful tool for providing timedependent properties of subtle variations of the crust.Various studies have widely discussed temporal seismic velocity changes (dv/v) in fault zones using noise correlations (i.e., Brenguier et al., 2008;Chaves & Schwartz, 2016;Rivet et al., 2011;Wegler et al., 2009).Yet, Iceland has exceptional conditions, that is, frequent volcano-tectonic activity, broadly distributed strong ambient noise sources, and multidisciplinary monitoring network, and can provide a unique opportunity to broaden our current knowledge in the field.Previous research (Donaldson et al., 2019) has shown the seismic velocity changes induced by a dike intrusion and graben formation in Holuhraun, north of Vatnajökull glacier, as well as the significant environmental influence on dv/v.The seismic velocity changes have been also reported in connection with repeated magma intrusions on the Reykjanes Peninsula and have since been used as a complementary tool for volcano monitoring in Iceland (Cubuk-Sabuncu et al., 2021).
We investigate earthquake-induced seismic velocity changes around the epicenters of the 2008 Ölfus earthquakes and we jointly interpret the main characteristics of the seismic signals together with the surface deformation and other available observations.To the best of our knowledge, this is the first time that temporal changes of seismic wave velocities are evaluated for a significant earthquake in Iceland, using ambient seismic noise.

Seismic Data
Continuous data were collected from three permanent seismic stations operated by the Icelandic Meteorological Office (IMO) around the source region of the 2008 Ölfus earthquakes.At that time, most seismic stations of the network were operating in triggered mode and were not configured to collect continuous data (see Figure 1).The seismic site SOL is located about 8 and 12 km east of the epicenters, respectively (Figure 1).SAN and KAS, the two other stations used in this study are located 25 and 35 km west of the eastern epicenter on the Ingólfsfjall fault, respectively.The stations were equipped with Lennartz 3D, 5-s sensors.We checked the data set for gaps, clock, and other instrumental errors to ensure waveform quality.We used the IMO earthquake catalog (https://skjalftalisa.vedur.is) to examine the overall seismicity.

Geodetic Data
The crustal deformation associated with the 2008 Ölfus events has been analyzed in detail in several studies (e.g., T. Árnadóttir et al., 2018;Decriem et al., 2010;Hreinsdóttir et al., 2009).A sparse network of continuous GPS (hereafter "cGPS") stations were operating in southwest Iceland in May 2008, with two stations in the immediate epicentral area (HVER and SELF, Figure 1).The whole GPS campaign network in the SISZ was re-measured following the 2008 mainshocks, and instruments were deployed in semi-continuous mode at several of the benchmarks (Decriem et al., 2010).Here, we use the time series of the horizontal components from three GPS stations, that is, two cGPS stations closest to the epicentral area (HVER, and SELF) and a campaign GPS station (HAAL) which is co-located with one of the seismic stations (SOL, see Figure 1).The data processing is described by T. Árnadóttir et al. (2018).We present detrended time series of GPS displacements from HVER, SELF, and HAAL, assuming an average station velocity for 2001-2008 in the ITRF08 reference frame (Altamimi et al., 2011).

Seismic Ambient Noise Correlations
Using a single-station approach, we computed ambient noise cross-correlations between all three components of the seismic records.The single-station cross-correlation analysis is a robust method that has been widely applied to quantify earthquake-induced temporal seismic velocity variations (e.g., Caudron et al., 2021;Cubuk-Sabuncu et al., 2021;De Plaen et al., 2016, 2019;Wegler & Sens-Schönfelder, 2007).The MSNoise software (Lecocq et al., 2014) was used to calculate noise cross-correlation functions (NCFs) and relative seismic velocity variations (dv/v) of its coda, which represents multi-scattered data in a volume beneath the station.The trend and mean of the data were removed prior to applying pre-processing band-pass filtering between 0.01 and 8.0 Hz.Seismograms were recorded at a rate of 100 samples per second and were downsampled to 20 Hz for our analyses.Cross-correlation calculations were carried out in four different frequency bands (first 0.1-3.0Hz, then 0.4-0.7,0.7-1.7,and 1.7-3.0Hz).We normalized the waveforms with an amplitude clipping threshold equal to three times the root mean square level (Tukey, 1962) to identify and remove the effects of intermittent strong signals on waveforms and applied spectral whitening.To obtain more stable NCFs, we constructed daily NCFs by stacking 10-min slices and then we computed 10-day stacked traces from 1-day-long data.We demonstrated the performance of four different stacking durations (1-, 2-, 5-, and 10-day) in Figure S1a in Supporting Information S1.Our stacking tests suggested that using a longer stack improved the signal-to-noise ratio (SNR) and the overall Continuous seismic data were only available at the sites denoted by black triangles during the analysis period.Stations highlighted in green were in trigger mode during operation and utilized to locate earthquakes; they were not used in our dv/v analyses.Yellow stars represent M ≥ 5.8 earthquakes, and circles represent background seismicity in 2008 (M ≥ 1) both before (blue) and after (red) the Ölfus doublet.The weather stations are also highlighted in orange color.We used earthquake focal mechanism solutions from the ISC Bulletin.The aftershocks inside the black box are compared with seismic velocity changes also used to calculate the total number of earthquakes (see Figure 3a1).A red box indicates the study region on the upper left inset map.correlation in the three frequency ranges we employed.Therefore, to enhance temporal stability, we used 10-day stacks and carefully interpreted the dv/v fluctuations within a 10-day time window range to avoid potential processing effects on shorter time scales.We calculated NCFs for the NZ, EZ, and EN component combinations by cross-correlating different channels of a single three-component seismic station (N, E, and Z are the north, east and vertical components, respectively, of the seismic recording).We analyzed data from January 2008 to January 2009, using the period January-May 2008 as a background reference prior to the onset of the earthquakes.

Measurement of dv/v Using Seismic Interferometry
We applied two methods, both of which rely on previously calculated cross-correlations and are known to be robust in determining seismic velocity variations.We first employed a frequency-time-based procedure to better assess the frequency-dependence of the dv/v measurements.We then implemented a time-domain technique that is widely used in noise-based seismic interferometry.Utilizing both techniques for the analysis of data allows for a more accurate assessment of seismic velocity changes.
The wavelet method (Mao et al., 2019) is a recently developed technique to determine seismic wave travel time differences with high time and frequency resolution.We provide a brief overview of the workflow; the reader is referred to Mao et al. (2019) and Mordret et al. (2020) for a complete description of the technique.Our analyses were based on calculating the continuous wavelet transforms (CWT) of reference and current waveforms, which provided phase and amplitude information in the time-frequency domain.For each station, we calculated the reference CCF by stacking 5 months of data before the doublet, and the current was the daily data.We used the cross-wavelet transform (XWT) (Grinsted et al., 2004) to calculate the time-frequency domain phase difference between these two signals and derived the time shift (dt) at various frequencies.We computed the time delays for different cross components (NZ, EN, EZ) over a specified period and a frequency range (0.1-3.0 Hz), we then used linear regression on time shifts generated by the wavelet approach for determining dv/v.The method allowed the quantification of multifrequency dv/v (0.1-3.0 Hz) at once, and we presented our results as the average dv/v of three cross-components for a single station (Figure 2).Since the cross-wavelet transform analysis enabled us to determine frequencies at which the significant velocity changes take place, we then extracted band averaged dv/v time series in chosen frequency ranges based on our observations (0.4-0.7 Hz, and 0.7-1.7 Hz).
The stretching method (e.g., Sens-Schönfelder & Wegler, 2006) was further used to obtain temporal velocity variations after calculating the cross-correlations at the same three frequency bands.The stretched (or compressed) daily correlations were compared to the reference stack, and dv/v were allowed to change up to ±2 percents to investigate the larger variations.We used scattered seismic waves in the coda, with a time window of 5-25 s, as coda waves are less sensitive to noise source variations (Colombi et al., 2014;Hadziioannou et al., 2009).Figure S1b in Supporting Information S1 demonstrates the variations in relative seismic velocity as a function of various lag times.Before selecting the coda lag time window for our overall analysis (5-25 s), we checked the highest correlation coefficient values and consistency in time-series at different lags.After stretching, seismic velocity variation time series for each frequency range were averaged over different cross-components (EZ, NZ, EN) to reveal the common features among data.A rolling mean of 3 days was applied to both GPS and dv/v to smooth-out small fluctuations in the data series.

Spatial Sensitivity of dv/v Measurements
The seismic coda consists of scattered body and surface waves due to heterogeneities in the crust (Aki, 1969;Aki & Chouet, 1975).The contribution of body and surface waves to depth sensitivity varies with time in the coda (Obermann, Planès, Larose, Sens-Schönfelder, & Campillo, 2013;Obermann et al., 2016).Most of the surface wave energy is concentrated in the early part of the coda and is sensitive to changes in shallow layers (Obermann, Planès, Larose, Sens-Schönfelder, & Campillo, 2013).Surface-to-body wave sensitivity transition occurs at late lag times in the coda, around six times the mean-free time (Obermann, Planès, Larose, & Campillo, 2013;Obermann, Planès, Larose, Sens-Schönfelder, & Campillo, 2013).The mean-free time (t * ) is the average amount of time a wave spends in the medium before it scatters.The mean-free time can be calculated using a mean wave velocity (v) of 2 km/s and t * = l * /v, where l * is the mean-free path.We can anticipate a mean-free path of 24-30 km, which is comparable to previous studies (e.g., Sánchez-Pastor et al., 2019) and therefore 6t * ∼ 72-90 s.This illustrates that our coda analysis windows for dv/v calculations (5-25 s) are in the early part of the coda.
Hence, to estimate the depth of observed seismic velocity changes, we calculated sensitivity kernels of surface wave phase-velocities (Figure S2a in Supporting Information S1), assuming a dominant surface-wave contribution in the coda windows we analyzed.The crustal S-wave velocity model of Tryggvason et al. (2002), a representative average 1D velocity structure in SW Iceland, was used to compute Love and Rayleigh wave sensitivity kernels.Figure S2a in Supporting Information S1 displays the Love and Rayleigh wave sensitivity at individual periods ranging from 0.3 to 10 s (0.1-3.0 Hz).Based on this surface-wave assumption, the frequencydependence of the dv/v results indicate that the higher frequencies are sensitive to shallower parts of the crust.Longer periods (or lower frequencies) exhibit larger wavelengths and enable sampling through the medium at greater depths.In the high-frequency band of 1.7-3.0Hz, dv/v is sensitive to depths less than a kilometer.Depth sensitivity of frequencies 0.7-1.7 Hz can reveal changes at 1-1.5 km depths, whereas frequencies 0.4-0.7 Hz appear to be more sensitive to depths of around 2.5-3 km.
To identify dv/v-sensitive regions around seismic stations, we computed simple 2D scattering kernels of coda waves at various lags (10-25 s) and mean-free paths (e.g., 5-200 km) (Text S2, Figures S2b and S2c in Supporting Information S1) and employed uniform scattering assumption for efficiency following the methodology of Paasschens (1997) and Planes (2013).Ongoing research indicates that more complex techniques (e.g., van Dinther et al., 2021) may offer a detailed representation of laterally varying scattering properties than simplified uniform models.However, such models require extensive (often unavailable) medium information and are beyond the scope of this work.Visual inspection on the generated 2D spatial kernels for a homogenous medium (Figures S2b and S2c in Supporting Information S1) showed that both mainshock epicenters are located within the region of sensitivity for SOL.This confirms the station's sensitivity to small, earthquake induced medium changes, with higher sensitivity (K:0.1)near the station location, decreasing (K:0.01) with distance.Considering this sensitivity and the observed apparent dv/v of 0.6% at SOL, the true velocity changes at the epicenter should be in the order of tens of percent yet accurate quantification remains challenging without complex scattering assumptions in the 2D model.

Results
We summarize the temporal variations in relative seismic wave velocities for the SOL seismic station (Section 4.1) using both the wavelet (Section 4.1.1)and stretching (Section 4.1.2)methods.The dv/v results for the other two stations (SAN, KAS) are presented in Section 4.2.We report average seismic velocity variations over ±0.2% at the time of the earthquake, and we consider changes within the next 10 days as co-seismic dv/v to avoid possible post-processing effects on the signal onset (e.g., Figure S1a in Supporting Information S1).The time period covering the recovery of the co-seismic drop and decorrelation in dv/v until September 2008 is referred to as post-seismic dv/v.In the supplementary material, we further present the seasonal background dv/v variations (<0.25%) in the study area for longer-time series (Figures S3a-S3e in Supporting Information S1), the dv/v variations for individual cross-components (Figures S4a-S4c in Supporting Information S1), and the exemplary CCFs (Figures S5a-S5c in Supporting Information S1).The cross-wavelet analysis for KAS and SAN are also provided in Figures S6a and S6b in Supporting Information S1.

Wavelet-Based Seismic Velocity Variations at SOL
We applied cross-wavelet analysis in the frequency band of 0.1-3.0Hz and investigated how relative seismic wave velocities changed at each frequency.As shown in Figure 2, we detected dramatic changes in seismic velocities caused by the 2008 Ölfus doublet using wavelet transforms.A rapid co-seismic velocity decrease up to 0.7% in average dv/v values was noticeable in a broad frequency range from 0.4 to 2.0 Hz, which lasted for about few weeks until mid-June.This suggests that the uppermost 3 km of the crust in the epicentral area were affected by the short-term velocity variations that followed the doublet.The formation of new crustal fractures and strong earthquake activity may account for this seismic velocity decrease seen over wide range of frequencies.
The magnitude of the seismic velocity drop varied with frequency.For instance, in June, seismic velocity decrease was large (>0.5%) at mid-frequencies between 0.7 and 1.7 Hz (1-1.5 km depth), but it was around 0.2% above 2 Hz.The weaker velocity decrease at the highest (2.5-3.0Hz) frequency range than at 1 Hz suggested that changes to the near-surface were rather small.The velocity drop was about 0.35% at crustal depths around 3 km, as indicated by 0.4-0.7 Hz.
Another intriguing finding was that the seismic velocity decrease in the frequency range of 0.7-1.7 Hz was not fully recovered over the time of our observation period (Figure 2).We found that the high frequencies above 1.7 Hz showed signs of velocity recovery post-seismically until mid-September.The lowered velocities increased by about 0.4% in the 0.7-1.7 Hz frequency range in the 3 months following the co-seismic period, but then again decreased by 0.3%.Sensitivity kernels (Figure S2 in Supporting Information S1) implied that the depth at which dv/v drop occurred for frequencies between 0.7 and 1.7 Hz was most likely limited to a thin layer around 1 km in the upper part of the crust.As there were no other large earthquakes or aftershocks throughout this period, this velocity change may indicate the presence of some long-term effects of the earthquakes in the medium until the end of the observation period.Such shallow quasi-permanent changes may result from temporal changes in groundwater levels or water diffusion properties in the presence of new fractures in the structure after intense earthquakes.Localized groundwater level changes observed after the 2008 Ölfus doublet support a possible link between quasi-permanent dv/v variations and hydrological changes.For instance, water level rise in a borehole near the SOL station (Þorbjörnsson et al., 2009), consistent with the compression predicted by the stress model and focal mechanism, potentially contributing to a reduction in wave speeds.However, limitations in existing groundwater level data (e.g., non-continuous measurements, variable borehole depths and lack of geological corrections) preclude any further correlation between dv/v and groundwater levels.
To compare the resulting dv/v with the long-term seasonal variations, we used an additional seismic data set and extended our wavelet-based dv/v calculations for the years 2019-2021.We also compared dv/v to precipitation and temperature (Text S3 in Supporting Information S1).The Figure S3c in Supporting Information S1 indicates that periods of seismic velocity increase at SOL occur from May to October (2019-2021), which is followed by a decrease in seismic velocity between October to May (2019-2021) (based on cross-wavelet analysis).This demonstrates that the strong decrease (>0.2%) in seismic velocities from June to the end of 2008 (Figure 2) is not seasonal but primarily emerges from the medium after the earthquakes, and a possible seasonal contribution may only occur on the order of ±0.1%-0.2%.We further notice the influence of the sources in the frequency range of 0.1-0.3Hz (e.g., oceanic waves).Thus, dv/v for frequencies below 0.3 Hz is excluded from our overall interpretations and its further investigation is beyond the scope of this paper.

Overview of Seismic Velocity Change Time-Series at SOL
We detected a co-seismic velocity decrease (0.65%) at the seismic station closest to the main shock epicenters (SOL) following the M w 6.1 doublet on 29 May. Figure 3 summarizes (a) seismicity and earthquake rate, (b) seismic velocity changes, and (c) deformation.Wavelet-based dv/v time-series displays seismic velocity changes as red thick lines, whereas scatter plots reveal stretching-based dv/v (Figure 3).We present low and intermediate frequency band dv/v in Figures 3b1 and 3b2, also high frequency dv/v in Figure S7 in Supporting Information S1.
The correlation coefficient (Corr.Coeff.) of the stretching changes between the selected reference period (January to mid-May 2008) and the stacked seismic NCFs.Bright colors with correlation coefficients close to 0.9-1.0indicate highly correlated waveforms.In contrast, dark colors with correlation coefficients below 0.7 indicate decorrelation, which means less similarity between daily and reference NCFs.The strength of correlation declines with decreasing values (Table S1 in Supporting Information S1).Throughout the paper, we address decorrelation in relation to structural changes within the medium, whereas changes in dv/v amplitudes are attributed to local velocity variations (temporal phase-shifts).
The seismic velocity variations at SOL are likely to reflect changes in the medium near the source area due to its proximity to the rupture zone (8-12 km, Figure S2b in Supporting Information S1).In the low frequency band (0.4-0.7 Hz), the co-seismic velocity decrease at SOL was relatively small (∼0.35%), and the earthquake-induced dv/v drop recovered rapidly, returning to pre-quake values less than a month (Figure 3b1).The most significant velocity changes were observed in the 0.7-1.7 Hz frequency band calculations (Figure 3b2).The average dv/v was reduced by 0.65-0.7%after the two mainshocks (Figure 3b2, pink shaded area).The amplitude of the coseismic dv/v drop was restored by 0.4% between mid-June to July (blue dashed line), suggesting velocity recovery in the medium.However, an offset of 0.2% remained between the dv/v time series before and after the earthquake (0.7-1.7 Hz) and dv/v fluctuated below zero for the rest of the observation period (Figure 3b2).After September, another striking decrease in velocity of 0.3% was observed mainly in the 0.7-1.7 Hz frequency band, as also seen in the Figure 2.These changes in the 0.7-1.7 Hz frequency range should occur in a volume around 1 km depth beneath the station as indicated by sensitivity kernels (Figure S2 in Supporting Information S1).We suggest that this long-term post-seismic velocity change is limited to a crustal layer at 1 km depth and can be attributed to redistribution of crustal fluids over time after the earthquakes.
Earthquakes caused temporary loss of correlation between the reference and the current waveforms.In the 0.4-0.7 Hz frequency band, the decline in the correlation coefficient was relatively lesser (0.4), than in other frequencies, and it restored back to pre-earthquake levels much faster (before July) (Figure 3b1).As lower frequencies are more sensitive to deeper crustal parts, this may imply that any structural changes occurring at depths around 2.5 km at SOL were not permanent (e.g., Figure S2a in Supporting Information S1).Similarly, at the time of the earthquake, the correlation rapidly dropped from ∼1.0 to 0.7 within the 0.7-1.7 Hz range but returned to its pre-earthquake levels (>0.7) in September.This recovery indicated no significant permanent structural damage around the station.We defined this distinct 3-month period as the post-seismic period, based on our observation of waveform decorrelation in the dv/v time-series (Figures 3b1 and 3b2).
In the high-frequency band (1.7-3.0Hz) stretching calculations, dv/v dropped by 0.4% at SOL (Figure S7 in Supporting Information S1).The correlation coefficient also rapidly declined from 0.8 levels to 0.1 co-seismically and it restored to 0.7 within 3 months.The dv/v variations in the 1.7-3.0Hz band were similar to those in other bands.However, detailed interpretations of near-surface changes are limited due to low correlation coefficients below 0.7.
Figure 3 further demonstrates that the temporal variations of two series obtained by the stretching and the wavelet methods are highly similar.Since dv/v time-series represent the average response of frequencies within a given band, which is defined by upper and lower frequency limits, using the wavelet method to select relevant frequency bands to extract dv/v curves offers a significant advantage for reliable dv/v interpretations.We compared the seismic velocity variations to meteorological data to evaluate the seasonality in the data set.Figure 3d1 displays the daily rainfall and temperature variations observed at station KALFH, located 21 km east of SOL, for the year 2008.The region experienced its highest annual rainfall between late September and October in 2008, as the average daily rainfall surpassed 1.2 mm in mid-September.It is known that dv/v can decrease after heavy rainfall (or snow meltwater) due to pore-pressure increase and saturation (e.g., Donaldson et al., 2019).At SOL, precipitation-induced decrease can be relevant in mid-September based on data, and dv/v amplitudes slightly decreased at that time, around 0.1% at the lowest frequency range and 0.2% in the 0.7-1.7 Hz frequency band.The temperature data exhibited a consistent pattern, with a rapid increase beginning in April, peaking in July and August, and then progressively declining in fall and winter.Notably, the long-term data (Figures S3a and  S3b in Supporting Information S1) also showed this seasonal trend in temperature.We observed that the annual increase of dv/v at SOL coincided with the increase in temperature over the summer months, up to 0.25% (Figure S3b in Supporting Information S1), likely to be related with thermally induced stresses near the surface in the dry season.Particularly, in the lowest frequency range (0.4-0.7 Hz), positive dv/v amplitudes reaching up to 0.25% were observed between April and August, while in the intermediate frequencies (0.7-1.7 Hz) the maximum dv/v amplitude were less than 0.25%.The time series of dv/v in the highest frequency range (1.7-3.0Hz) did not exhibit any notable patterns in relation to precipitation and temperature (Figures S3a and S3b in Supporting Information S1).Therefore, overall data set (2008 and long-term dv/v) suggests that annual changes were on the order of or less than ±0.25% in our dv/v calculations.

Temporal Seismic Velocity Changes With Distance
We analyzed dv/v from two other available seismic stations, SAN and KAS, to explore spatial variations in velocity following the Ölfus doublet (Figure 4 and Figures S6a and S6b in Supporting Information S1).At lower frequencies (0.4-0.7 Hz), station SAN experienced relatively larger co-seismic velocity decrease (0.3%) compared to KAS (0.2%) (Figures 4a and 4b).This decrease rapidly recovered until July, reaching near pre-quake levels by September (Figures 4a and 4b).Following the end of the post-seismic period in mid-September, we observed another 0.2% velocity drop at SAN (0.4-0.7 Hz).Both stations experienced slightly larger co-seismic velocity drop in the intermediate frequencies (0.7-1.7 Hz) than in the low frequencies (Figures 4c and 4d).For instance, the stretching time-series for SAN revealed a 0.8% decrease in velocities.KAS revealed a lower co-seismic dv/v decrease of 0.4% compared to SOL and SAN (0.65-0.8%) likely due to its distance to epicenters.The dv/v amplitudes at KAS fully recovered after the earthquake until mid-September that same year, however those at SAN only partially returned to pre-quake levels, showing an increase about 0.6% (Figure 4d, yellow shaded area).Our high-frequency wavelet calculations showed a velocity decrease of 0.8% at SAN in the 0.7-2.0Hz frequency range, and a decrease of less than 0.4% for KAS (Figures S6a and S6b in Supporting Information S1).The stretching calculations within the 1.7-3.0Hz band revealed a 0.4% dv/v drop at KAS (which recovered over time), however interpreting dv/v changes for near-surface at both stations is not straightforward due to low correlation coefficients (below 0.7).
The decorrelation at low frequencies (0.4-0.7 Hz), was slightly higher at SAN compared to KAS, both stations experienced a decrease in the correlation coefficient from around 1.0 to 0.86 and 0.88 levels, respectively (Figures 4a and 4b).The decreased correlation coefficients in the low frequencies recovered in the early postseismic period until July, much faster than the recovery observed in the intermediate frequencies.The earthquake-induced decorrelation was more pronounced in the 0.7-1.7 Hz compared to the low frequency band calculations.For instance, we found a rapid decrease in the correlation coefficients from 0.9 to 1.0 levels to 0.75 (KAS) and 0.7 (SAN) co-seismically.The co-seismic correlation decrease (also velocity) at KAS fully restored by September (Figure 4c and Figure S6b in Supporting Information S1).Notably, the decorrelation for SAN was not recovered to pre-quake levels (0.95-1.0) completely, indicating substantial changes within the subsurface.
Assuming that the seasonal daily rainfall was in its seasonal highest between late September to October (2.5 mm HELLS), the dv/v decrease of 0.2% at SAN in the low frequency range of 0.4-0.7 Hz, is likely to coincide with the precipitation increase in late September (Figures 4e and 4f).The local precipitation recorded at REITR (0.6 mm) was less than half of the amount recorded at HELLS, with a peak also occurring in December.Therefore, in winter, dv/v amplitudes exhibited minor fluctuations around 0.1%, suggesting that slight changes may be associated with precipitation, yet no strong correlation between precipitation and dv/v at KAS was observed in the 2008 data (Figures 4e and 4f).The potential dv/v reduction caused by rainfall appeared to be less pronounced within the 0.7-1.7 Hz frequencies.The periods of heaviest rainfall (0.6 and 2.5 mm) over a few days did not coincide with the dv/v reduction observed for two and a half months between late September and December 2008 (Figures 4c-4f).Based to the dv/v calculations obtained at KAS and SAN, the relationship between the seasonal seismic velocity change and temperature was not obvious.However, small seasonal positive dv/v amplitudes less than 0.2% were seen at KAS in the lowest frequency range (0.4-0.7 Hz) particularly in the 2016 time-series data.
The analyses of all three seismic stations used in this study revealed that the earthquake induced relative seismic velocity change occurred over a 40 km spatially wide zone.Both SAN and SOL stations showed semi-permanent crustal changes at a shallow layer (around 1 km), whereas the rapid recovery of velocities and correlation at KAS, particularly, suggests that the impact of earthquakes on the crust decreases with distance.

Sensitivity of dv/v to Aftershocks
The May 2008 Ölfus mainshocks were followed by 3,200 earthquakes with local magnitudes greater than M ≥ 1 (Figure 1).Earthquake aftershocks were clustered in two parallel N-S oriented faults (depths 8 and 9 km) (Brandsdóttir et al., 2010).Seismicity was also triggered along an E-W trending fault which marks the plate boundary and even further west on the Reykjanes Peninsula (see Figure 1) (Brandsdóttir et al., 2010).This extensive aftershock triggering also indicated that the spatial triggering was well above what was expected from an M w 6.1 doublet.Within an hour of the main earthquakes, 40 foreshocks (M ≤ 3.5) were registered.Approximately 40 large aftershocks with magnitudes up to M3-4.75 were reported on the western fault (Kross, Figure 1) shortly after the two main earthquakes (Brandsdóttir et al., 2010).
We compare the dv/v time-series obtained from SOL (0.7-1.7 Hz) with aftershocks of the 2008 Ölfus doublet.In Figure 3a1, we show the daily earthquake rate (M ≥ 0.1), magnitude distribution (M ≥ 1), and the cumulative number of earthquakes over a year (M ≥ 1). Figure 3a2 displays the histogram of aftershocks (M ≥ 0.1), in daily bins for a shorter time window; 8 days before the mainshock until July 1 and compares it with dv/v time-series at frequencies with the largest changes (0.7-1.7 Hz).
The relative seismic velocities began to decrease after the doublet on 29 May.More than 600 events were observed per day within a few days following the mainshock, and aftershocks decayed with time (Figure 3a2).Starting on 9 June, the dv/v values steadily increased and the number of aftershocks per day further declined to 100 in mid-June, possibly indicating the beginning of crustal recovery.We found that seismic velocities increased between mid-June and August at lower frequencies (0.4-0.55 Hz), possibly coinciding with declining aftershock activity at the same timescale (Figures 2 and 3).There were 2,700 events (M ≥ 1) recorded between January and September 2008 (IMO earthquake catalog), however seismic activity around Ölfus began to diminish dramatically in the fall.Between September 2008 and January 2009, only 450 earthquakes (M ≥ 1) were observed (Figure 1, black box).With reduced earthquake rate in September, the correlation coefficient increased above 0.7 level (Figures 3b1 and 3b2).
We further notice triggered seismic activity along the plate boundary in the Reykjanes Peninsula in our data set (Figure 1, west of 21.8°W), starting at the same time as the 2008 Ölfus activity.The seismicity suggests that aftershocks do not only occur on the main faults but the mainshocks may have generated transform motion also in the west along the plate boundary of the Reykjanes Peninsula.This remote deformation (>20 km) may not be noticed by GPS data since it may be below the observation thresholds.It does, however, support the dv/v findings observed by the westerly seismic sites (SAN and KAS).
The concentration of seismic energy is around 3 Hz or higher for most of our data, according to aftershock spectrograms of earthquakes (M > 3) recorded at SOL.Large aftershocks (e.g., M > 4.5) can affect correlation calculations, but we mitigate this by using normalization techniques (Bensen et al., 2007).As described in Section 3.1, the daily cross-correlations largely consist of background seismic noise and are not expected to be contaminated by transient noise due to earthquakes.As a result, we suggest that the complex physical changes in the medium related to mainshocks are likely the main reason for velocity reduction (more than 0.2%) rather than aftershocks interfering with the noise field.

Deformation Associated With the Ölfus Doublet and Seismic Velocity Changes
Several studies reported the crustal deformation caused by the 2008 Ölfus earthquakes.Hreinsdóttir et al. (2009) provided the first GPS-based fault models of the earthquake sequence.Decriem et al. (2010) performed a joint analysis of co-seismic surface displacements observed by GPS and InSAR, to obtain fault models (location, geometry and slip) for the two main events.They also estimated the static Coulomb stress changes from these models.T. Árnadóttir et al. (2018) indicated spatial changes in geodetic strain rate across a large area near Ölfus as well as longer-term evolution of GPS velocities in the SISZ.Our objective is to compare the main features of dv/v variations with reported co-seismic and post-seismic deformation at selected near-field GPS sites (Figures 1  and 3c1-3c5).
Significant deformation occurred within a 10-15 km epicentral zone (Decriem et al., 2010), and due to proximity to ruptures, the cGPS stations HVER and SELF recorded the largest co-seismic offsets (about 200 mm in Decriem et al., 2010;Hreinsdóttir et al., 2009).We show data from GPS sites 20 km northwest (HVER, Figures 3c1-3c3) and 4.5 km west (SELF, Figure 3c5) of the nearest seismometer SOL.A campaign GPS station, HAAL was colocated with the seismic station SOL, but as this area was close to nodal planes in the deformation pattern, the coand post-seismic displacements at HAAL were small (Figure 3c4).The co-seismic velocity decrease induced by the M w 6.1 earthquake doublet correlated with co-seismic horizontal motion recorded by HVER and SELF (pink shaded areas in Figures 3b and 3c).Following the main earthquakes, the horizontal GPS time series at HVER exhibited a distinct post-seismic surface displacement that lasted until August 2008 (total displacement of about 30 mm), and the signal strength abated by the end of the year.The post-seismic decorrelation observed at SOL and remote sites (<40 km) over three months was consistent with the period of post-seismic deformation (Figures 3  and 4).Beginning in mid-June, the rapid post-seismic decay at HVER slowed down; similarly, dv/v began to recover around the same time, implying that those changes might be driven by similar crustal processes.Decriem (2011) modeled the post-seismic deformation of the 2008 Ölfus earthquakes and deduced that the postseismic deformation observed in the first two months at sites close to the mainshock fault (e.g., HVER, SELF) was caused by shallow afterslip in the upper crust (<2 km).Their modeling also suggested viscoelastic relaxation to explain surface deformation recorded post-seismically at far-field locations, on a longer time scale.However, models did not indicate any significant deformation signal that could be attributed to poroelastic relaxation.Continued slip on the main fault(s) (i.e., afterslip) changes the stress in the crust, although the afterslip is much less, and hence promotes smaller stress changes than the mainshocks.The near-field seismic and deformation data indicated significant crustal changes caused by the 2008 Ölfus doublet, with similar patterns emerging at comparable time scales.We suggest that dv/v increase at all stations is likely to reflect the post-seismic crustal healing process of rocks at various depths (<3 km).The post-seismic afterslip in the fault zone may also have an impact, particularly on the nearest seismic station (SOL).

Changes in Volumetric Stress and dv/v
The seismic velocity can be affected by the change in stress or strain due to crack opening/closing (Nishimura et al., 2000;Poupinet et al., 1984).To explore possible factors causing the observed velocity decrease at seismic stations, we compared dv/v with modeled volumetric stress changes based on fault models of the two mainshocks (Hreinsdóttir et al., 2009).The spatial distribution of volume stress caused by the 2008 mainshocks is shown in the map and N-S cross-section view (Figure 5).The change of scale from blue to yellow indicates a transition from negative volume stress (contraction) to positive stress (dilatation), respectively.Modeling showed that the Ölfus doublet caused a contraction in the crust beneath SOL by 0.2 MPa up to depths of 8-10 km.The dilatation occurred in the northeast (also southwest) quadrants of the N-S oriented right-lateral strike-slip fault model, and the dilated region was centered 9-10 km north of SOL (Figures 5a and 5b).
In our calculations, the spatial distribution of the dilatational area (crack opening) does not overlap with the location of SOL.The model also predicts very small stress changes around SAN and KAS, since they are located further away (25-35 km) from the fault zones.There were no seismic sensors providing continuous data in the NE-SW dilatational zones at the time of the earthquake to allow for further comparison.Based on our model (Figure 5), we observe that there is no direct correlation between the measured co-seismic velocity decrease and the estimated co-seismic static stress changes.Therefore, it is likely that dv/v variations are a result of another dominant mechanism in the crust.

Possible Driving Mechanisms and Seismic Velocity Changes: Insights From Major Earthquakes
In this section, we compare seismic velocity changes caused by SISZ earthquakes to important findings from earthquakes of comparable magnitude in different tectonic areas.Temporal changes in crustal seismic wave velocities are driven by a variety of physical mechanisms, the most common of which are stress and strain changes, hydrological factors, damage within the crust due to crack opening and strong shaking.We review each of these components to better understand the temporal seismic velocity change caused by the Ölfus doublet in Iceland.
The co-seismic dv/v variations associated with moderate earthquakes have been widely documented over the past decade.These changes are frequently detected as a decrease in velocities, related to the opening of microcracks and damage near the fault zone.For instance, Wegler and Sens-Schönfelder ( 2007) and Wegler et al. (2009), detected a rapid 0.5% seismic velocity drop during the M w 6.6 Mid-Niigata, Japan earthquake at shallow and deep media (5 km) and attributed this decrease to crustal stress variations.Their findings supported the presence of a nonlinear site response in the shallow subsurface layer and fault zone damage.They stressed that water level changes couldn't explain the observed velocity reductions, and the exact mechanism remained elusive.
In California, the Parkfield section of the San Andreas Fault stands as one of the widely studied regions regarding temporal seismic velocity changes caused by earthquakes.The dv/v was reported to decrease by 0.04% and 0.08% immediately following the 2003 M6.5 San Simeon, and the 2004 M6.0 Parkfield earthquakes, respectively (Brenguier et al., 2008).Chunquan et al. (2016) further investigated the depth of changes and found S wave velocity reductions (0.02%) at 2.3 km depths caused by the San Simeon earthquake and more significant decreases (0.2%) at 1.2 km depths due to the Parkfield earthquake.These studies indicated to co-seismic material damage in the shallow layers caused by strong ground shaking (0.15 g).
Similarly, in 2019 the M w 7.1 Ridgecrest earthquake also resulted in dv/v drops ranging from near zero to 2.0% in a large network (Boschelli et al., 2021).As indicated by high-frequency autocorrelations, dv/v changed proportionally with its closeness to the fault rupture and the peak dynamic strain.Notably, stations located more than 50 km from the fault did not experience any reductions in dv/v.These observations are consistent with our results.For example, the most substantial dv/v changes we found are 0.65-0.8%reductions in co-seismic dv/v at nearby stations (12 and 25 km).We emphasize that despite the environmental challenges and sparse station coverage in the region, the distinct reduction in dv/v amplitudes (0.4%) as a function of distance (40 km) can be successfully observed.These findings are similar to those of Boschelli et al. (2021), who used a better coverage for their study.Our dv/v calculations provide depth sensitivity (<3 km) comparable to aforementioned studies.However, unlike some previous research, we avoid utilizing frequencies lower than 0.3 Hz to exclude possible source-related effects.
One distinctive aspect of our study is our use of cross-wavelet analysis (CWT) to evaluate daily dv/v changes at each frequency.This approach contributed to fully capturing broad changes within the 0.1-3.0Hz range, as well as to precisely selecting frequency intervals for time-series analysis, such as the 0.4-0.7 Hz or 0.7-1.7 Hz frequency bands.As a result, our analysis offers more insights into the frequency (and therefore depth) distribution of dv/v changes in relation to the Ölfus earthquake.

Relation to Stress Changes
Previous research has shown a correlation between dv/v variations and crustal dilatation and compression (e.g., Donaldson et al., 2019).New cracks frequently form in the crust as a result of dilatation, and opening of cracks leads to seismic velocity decrease, whereas contraction causes crack closure, which increases wave velocities (e.g., Nishimura et al., 2000Nishimura et al., , 2005;;Nur, 1972;Poupinet et al., 1984).Hence, it is anticipated that stations in stressincreasing areas would indicate dv/v increases, while those in regions exhibiting dilatation would experience coseismic dv/v decrease.In our volumetric stress change model, the observed co-seismic decrease in dv/v at the epicentral area (at station SOL) does not appear to correlate with the expected stress changes.Instead, we estimate contraction rather than dilatation in this area.Similar results have been reported, including findings related to the M w 6.6 Mid-Niigata earthquake, the M w 6.9 2008 Iwate-Miyagi Nairiku earthquake in Japan (Hobiger et al., 2012(Hobiger et al., , 2016;;Takagi & Okada, 2012;Wegler et al., 2009), and the 2014 M w 6.0 South Napa earthquake in California (Taira et al., 2015).These studies, based on their stress models and dv/v time series, identified co-seismic dv/v reductions in both compression and dilation areas, with no velocity increase in the compression zone.This suggests that the primary factor influencing the observed co-seismic dv/v is not static stress variations.According to Takagi and Okada (2012), it's likely that dv/v is influenced by damage in shallow layers caused by strong ground motion, which could hide the effects of static stress changes.In summary, these findings agree with the dv/ v detected at our stations, demonstrating that the co-seismic dv/v for the 2008 Ölfus earthquakes cannot be explained by volumetric stress changes.

Relation to Hydrological Effects
Changes in hydrological processes in the crust, whether caused by environmental factors or large earthquakes, can affect dv/v (e.g., Clements & Denolle, 2018;Lecocq et al., 2017).Large earthquakes can lead to rapid changes in groundwater levels (Sens-Schönfelder & Wegler, 2006), and pore pressure variations can also influence seismic velocity (Schaff & Beroza, 2004).In contraction areas around fault zones, cracks may close, reducing porosity and increasing seismic velocities.As porosity decreases in these contraction zones, the ground water levels can rise, resulting in increased pore-fluid pressure and decreased seismic velocities (Christensen & Wang, 1985).
For instance, Chaves and Schwartz (2016) reported a notable velocity reduction (0.6%) occurring at depths exceeding 5 km following the 2012 M w 7.6 Nicoya Peninsula earthquake in Costa Rica.These deep dv/v changes suggest that the variations in dv/v were likely influenced not by near-surface damage but by increased pore fluid pressure.In Iceland, the water level changes (around 10 m) were previously observed due to the 2000 M w 6.5 SISZ earthquake (Jónsson et al., 2003), and the recovery timing of this effect, around 10 days, is comparable to the turnaround of our relative dv/v changes.Due to a lack of water level data for the 2008 events, it is difficult to attribute the co-seismic dv/v drop at SOL (in the compressional area) to pore-pressure changes.Volumetric stress calculations further indicate negligible change around SAN and KAS (Figure 5), suggesting that near-surface hydrologic changes are unlikely to be the primary cause of such velocity reductions across a large volume of the crust.However, the spatially extensive aftershock sequence indicates that the crust on the Reykjanes Peninsula was indeed affected by the doublet.

Relation to Physical Damage
During co-seismic fracturing, structural deformation occurs as microcracks open, lowering elastic moduli in rocks and seismic velocity adjacent to the fault zone (Brenguier et al., 2008;Li et al., 1998;Rubinstein et al., 2007).Damage surrounding the fault heals over time, and velocities show reversal post-seismically (Brenguier et al., 2008;Wegler & Sens-Schönfelder, 2007).Microfracturing also declines with distance, and the spatial distribution of damage is well known to be limited to a narrow zone away from the fault (Ben-Zion et al., 2003;Mitchell & Faulkner, 2009).Recent findings indicated that the fault damage-zone can be several kilometers wide (1-3 km) (Cochran et al., 2009;Scholz et al., 1993;Shelef & Oskin, 2010).Physical damage and weakening in the fault zone caused by fault motion may contribute to observed co-seismic dv/v decreases in the Ölfus area.However, since all three stations used in the analysis are situated more than 8-12 km from the fault, co-seismic dv/v drops resulting from near-surface damage alone cannot account for the observed seismic velocity reductions.

Relation of Co-Seismic dv/v to Strong Ground Motion
The decrease in co-seismic velocity following major earthquakes exhibits a good correlation with the strong ground motion (Ikeda & Tsuji, 2018;Peng & Ben-Zion, 2006;Rubinstein & Beroza, 2005;Schaff & Beroza, 2004;Viens & Van Houtte, 2019;Viens et al., 2018).Numerious studies have established connections between dv/v changes and peak ground velocity (PGV), acceleration (PGA), and dynamic strain (PDS) parameters.The intense shaking of the ground can lead to physical damage in the near-surface layers, and substantial evidence has been observed at various sites.Soils exhibit non-linear behavior when subjected to accelerations surpassing about 0.1-0.2g threshold (Beresnev & Wen, 1996).In Japan, high peak ground acceleration values ranging from 0.2 to 0.8 g were recorded due to the 2004 Mid-Niigata earthquake (Wegler et al., 2009).Furthermore, during the 2008 M6.9 Iwate-Miyagi Nairiku earthquake, Takagi andOkada (2012), Hobiger et al. (2012) reported peak ground accelerations between 0.3 and 4.0 g.The authors noted a correlation between peak ground accelerations and high co-seismic velocity decreases, particularly in the range of 0.4-0.6% at frequencies between 0.5 and 1.0 Hz.In Chile, the M w 7.7 Tocopilla earthquake, with a peak ground acceleration of 0.25 g, led to a co-seismic velocity decrease of about 0.7% and demonstrated a linear relationship between the negative dv/v amplitude and peak ground acceleration (Richter et al., 2014).It is noteworthy that these prior studies also evidenced that the observed dv/v decrease in compressional stress regions where an increase is expected, is probably due to damage induced by strong shaking (e.g., Hobiger et al., 2012;Takagi & Okada, 2012).
Similarly, our data also indicated considerably higher sensitivity to previously reported PGA values associated with the 2008 Ölfus doublet.The peak ground acceleration was 0.4-0.66g within the near-fault region (Figure S8 in Supporting Information S1), especially within 3 km of the epicenter (Bessason, 2019;Halldórsson et al., 2010;Sigbjörnsson et al., 2009).PGA values ranging from 0.3 to 0.6 g were recorded in the towns of Selfoss and Hveragerdi (Figure S8 in Supporting Information S1), coinciding with the substantial dv/v drops at SOL and SAN (0.65%-0.8%, 0.7-1.7 Hz).With distances greater than 25 km, the strong ground motion was significantly decreased, that is, the PGA value of less than 0.1 g recorded in the capital area about 40 km west of the earthquake epicenters (Halldórsson & Sigbjörnsson, 2009;Halldórsson et al., 2010).The co-seismic dv/v amplitudes (0.4% at KAS) decreased as the peak ground acceleration declined farther from the epicenter, indicating a correlation between PGA, dv/v and the distance.
The co-seismic dv/v decrease observed in Iceland and other regions, as discussed in the studies mentioned above, is likely to be explained by strong ground motion that induce crustal damage.This results in a significant reduction in dv/v (about 0.65-0.8%),especially within the shallow subsurface layers (e.g., 1-2 km in our case).The transient strong shaking may also explain the rapid recovery of decorrelation and velocity following the Ölfus doublet, particularly in the lowest frequency range (0.4-0.7 Hz), as these frequencies correspond to deep subsurface layers (less than 3 km) that experience relatively less damage.

Post-Seismic dv/v Recovery
Following an earthquake, co-seismically decreased velocities may gradually recover or partially revert to their pre-earthquake values, indicating healing after crustal damage.For the Ölfus doublet, both velocities and decorrelation rapidly improved (within 1 month) in the lowest frequencies (0.4-0.7 Hz) especially at the most distant station KAS, indicating no permanent structural changes at layers of possibly around 2-3 km.At SOL, the closest station to the rupture area, decreased velocities also quickly recovered in the low frequency band.However, the dv/v drop was higher and didn't fully recover to pre-quake levels for the intermediate frequencies (0.7-1.7 Hz) during the calculation period.A similar observation was also made at SAN, where the dv/v drop consistently increased over time but partially recovered resembling the trend observed in SOL (0.7-1.7 Hz).The lowered correlation coefficients mostly recovered post seismically but not reached to 1.0 (e.g., 0.7-1.7 Hz, SAN).These findings suggest that the largest velocity changes took place in a particular shallow layer around 1-km depth (0.7-1.7 Hz), and the healing process possibly was much longer than the observation period in 2008.
Post-seismic healing resulting from the closure of cracks in the medium has been documented for example by Li et al. (1998), Peng and Ben-Zion (2006), Brenguier et al. (2008).However, both the amount of this recovery and the timescale of the healing process have significantly varied.For instance, after the M w 7.1 Ridgecrest earthquake, seismic velocity recovered observed over a period of 3 months, particularly in the areas with the highest slip (Boschelli et al., 2021).Some studies have suggested recovery times extending beyond several years, as no complete recovery detected during their analysis periods (e.g., Brenguier et al., 2008;Hobiger et al., 2016;Takagi & Okada, 2012).After smaller magnitude earthquakes (M5), complete post-seismic recovery within several months has been observed in Japan (Maeda et al., 2010).Regarding the previously reported post-seismic timescales, velocities for the Ölfus doublet appears to agree with other tectonic earthquakes.The fast recovery of decorrelation observed in the deeper upper crustal layers (less than 3 km) is likely related with the strong shaking, which primarily induces damage in much shallower regions.The semi-permanent decorrelation at SAN could suggest the presence of shallow and localized non-linear effects caused by shaking.The long-term dv/v decrease at SOL, occurring around 1-1.5 km depths until the end of the analysis period in 2008, could be indicative of temporal fluctuations in groundwater levels or water diffusion properties within a specific layer following the Ölfus earthquakes.

Conclusion
The May 2008 M w 6.1 Ölfus doublet caused co-seismic decrease in crustal seismic wave velocities, extending up to 40 km distance from the epicenter, indicating extensive impact far beyond the Ölfus region.Employing the cross-wavelet transform and stretching based dv/v analysis, we better characterized the co-seismic and postseismic velocity changes.Notably, we observed substantial co-seismic velocity decrease up to 0.65-0.8%at SOL and SAN, which are quasi-permanent within the 0.7-1.7 Hz frequency range.This prolonged crustal response may be associated with near-field changes in water levels or fluid flow confined to a shallow layer (around 1 km) surrounding the rupture area.
The dv/v time-series reveal a 3-month post-seismic recovery of decorrelation with a faster correlation increase at lower frequencies (0.4-0.7 Hz).This recovery in dv/v amplitudes and correlation may indicate healing of rock damage after the mainshocks.The high permeability of the upper crust also allows for a fast rebound of coseismically disturbed crustal fluids, potentially contributing to dv/v reversal in the early post-earthquake period.
The relative seismic velocity changes agree with co-seismic and post-seismic crustal displacement observed by GPS, yet no clear correlation is found between the co-seismic dv/v drop and the volumetric stress model.Crustal damage from strong ground motion is likely to explain the temporal co-seismic dv/v decreases, as indicated by the PGA distribution.
The Ölfus doublet appears to induce a transform motion along the plate boundary (westward), as evidenced by triggered far-field aftershocks on the Reykjanes Peninsula.This may also suggest the possible presence of "silent" slip, releasing tectonic stress without causing destructive earthquakes along the plate boundary in the peninsula.In conclusion, we demonstrated that fault zone monitoring in Iceland utilizing single-station analysis, with fewer monitoring stations and less data collection, is reliable and consistent with previous research.

Figure 1 .
Figure 1.Seismic stations and continuous GPS sites (pink squares) in the vicinity of the 2008 Ölfus seismic sequence.Continuous seismic data were only available at the sites denoted by black triangles during the analysis period.Stations highlighted in green were in trigger mode during operation and utilized to locate earthquakes; they were not used in our dv/v analyses.Yellow stars represent M ≥ 5.8 earthquakes, and circles represent background seismicity in 2008 (M ≥ 1) both before (blue) and after (red) the Ölfus doublet.The weather stations are also highlighted in orange color.We used earthquake focal mechanism solutions from the ISC Bulletin.The aftershocks inside the black box are compared with seismic velocity changes also used to calculate the total number of earthquakes (see Figure3a1).A red box indicates the study region on the upper left inset map.

Figure 2 .
Figure 2. Average dv/v (NZ, EZ, EN) for the SOL seismic station calculated using the cross-wavelet transforms in the 0.1-3.0Hz frequency range.Seismic velocity decrease and increase are indicated by the red and blue color scales, respectively.The time of the 2008 Ölfus doublet is marked by a black dashed line.

Figure 3 .
Figure 3. (a1) The magnitude-time distribution of M ≥ 1 earthquakes (blue circles).The red curve depicts the cumulative number of earthquakes (M ≥ 1).The yellow bars represent the daily earthquake rate (M ≥ 0.1), and the vertical scale is the same as in (a2), which compares the aftershocks in daily bins with the dv/v (0.7-1.7 Hz) between 21 May and 1 July.Panels (b1, b2) show the average seismic wave velocities for station SOL.The thick red lines represent dv/v based on the wavelet method.The correlation coefficient values from the stretching calculation are also presented at the bottom of each panel as solid black line, with the scale given on the right.The correlation coefficient color scale is further provided beside panel (b1).Daily stretching error is plotted with thin gray vertical bars on each figure.We label the analyzed frequency bands at the top of each plot.The red dashed lines mark the time of the M w 6.1 earthquake.Pink shaded areas indicate the co-seismic decrease, while blue dashed lines indicate the approximate times of velocity recovery for different frequency ranges.The post-seismic recovery timeline is also shown in yellow shaded areas.(c1, c2) Horizontal component deformation profiles from HVER.(c3-c5) GPS time series from HVER, HAAL, and SELF within the same time window as (a2).For details about the GPS processing, see T. Árnadóttir et al. (2018).(d1) The rainfall and temperature at KALFH as a function of time during 2008.

Figure 4 .
Figure 4. Temporal seismic velocity variations over a year at stations KAS (left) and SAN (right).(a, b) Show the dv/v values for both stations in the 0.4-0.7 Hz (d ≈ 2.5-3 km) frequency band.The dv/v for the 0.7-1.7 Hz frequency band (d ≈ 1-1.5 km) is given in (c, d), respectively.The red dashed lines denote the 29 May earthquakes, the pink shaded area represents the co-seismic dv/v reduction, and the yellow shaded areas indicate the post-seismic period.The dv/v color scale is the same as in Figure 3.The correlation coefficient is shown with a solid black line beneath each dv/v plot.Panels (e, f) depict 2008 rainfall and temperature data obtained from the nearest weather stations to KAS and SAN.

Figure 5 .
Figure 5. Volumetric stress changes (in MPa) due to the M w 6.1 Ölfus doublet, estimated using fault models from Hreinsdóttir et al. (2009).Dilatation is positive (warm colors) and contraction is negative (cold colors).(a) Horizontal distribution of stresses at 5 km depth near the N-S oriented Kross (W) and Ingólfsfjall (E) faults; (b) stress changes in a vertical N-S oriented depth section below SOL.Station locations are indicated by inverted black triangles (seismic) and green squares (GPS).White lines display the coastline and the location of mapped faults.