# Permanent, seasonal, and episodic seismic sources around Vatnajökull, Iceland, from the analysis of correlograms

Sylvain Nowé, Thomas Lecocq, Corentin Caudron, Kristín Jónsdóttir, Frank Pattyn

<figure><img src="image_1.png" /></figure>

## To cite this version:

Sylvain Nowé, Thomas Lecocq, Corentin Caudron, Kristín Jónsdóttir, Frank Pattyn. Permanent, seasonal, and episodic seismic sources around Vatnajökull, Iceland, from the analysis of correlograms. Volcanica, 2021, pp.135-147. 10.30909/vol.04.02.135147. insu-03594215

HAL Id: insu-03594215

https://insu.hal.science/insu-03594215

Submitted on 2 Mar 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Permanent, seasonal and episodic seismic sources around Vatnajökull, Iceland, from the analysis of correlograms

Sylvain Nowé*<sup>α,β</sup>, Thomas Lecocq<sup>α</sup>, Corentin Caudron<sup>γ</sup>, Kristín Jónsdóttir<sup>δ</sup>, Frank Pattyn<sup>β</sup>

## ABSTRACT

In this study, we locate and characterise the main seismic noise sources in the region of the Vatnajökull icecap (Iceland). Vatnajökull is the largest Icelandic icecap, covering several active volcanoes. The seismic context is very complex, with glacial and volcanic events occurring simultaneously and the classification between the two can become cumbersome. Using seismic interferometry on continuous seismic data (2011–2019), we calculate the propagation velocities and locate the main seismic sources by using hyperbolic geometry and a grid-search method. We identify and characterise permanent oceanic sources, seasonal glacial-related sources, and episodic volcanic sources. These results give a better understanding of the background seismic noise sources in this region and could allow the identification of seismic sources associated with potentially threatening events in real-time.

## RÉSUMÉ

Dans cette étude, notre but est de caractériser et localiser les sources de bruit sismique principales dans la région de la calotte glaciaire du Vatnajökull (Islande). Cette calotte est la plus vaste d'Islande et recouvre également plusieurs volcans actifs. Le contexte sismique local est très complexe, avec des événements glaciaires et volcaniques en simul-tané, rendant la classification difficile. Nous avons utilisé la méthode d'interférométrie sismique sur des données sismiques continues de 2011 à 2019. Nous avons étudié les corrélogrammes, calculé les vitesses de propagation et localisé les sources sismiques en utilisant la géométrie des hyperboles ainsi qu'une méthode de probabilités de lo-calisation. Une sismicité océanique constante, une saisonnière liée aux glaciers et une épisodique volcanique ont été identifiées et caractérisées. Ces résultats permettent d'acquérir une meilleure compréhension du contexte sismique de cette région et pourraient également permettre d'identifier plus rapidement une source sismique associée à un évènement potentiellement menaçant.

Keywords: Seismic noise; Correlograms; Seismic interferometry; Vatnajökull icecap; Seasonal seismicity

## 1 INTRODUCTION

Iceland is an Atlantic island resulting from the inter-actions between a mantle plume and plate spreading along the Mid-Atlantic Ridge, processes which have generated intense volcanism. The volcanic and seismic activity is primarily generated along the plate bound-ary across Iceland [Sigmundsson 2006]. Iceland's con-struction started approximately 24 million years ago [Thordarson and Larsen 2007]. The Icelandic hotspot is moving and currently located beneath the north-west part of the Vatnajökull icecap [Martin et al. 2011]. Approximately 30 volcanic systems are considered ac-tive on the island, each of them consisting of one cen-tral volcano associated with fissure swarms and, often, geothermal areas [Thordarson and Larsen 2007]. The

Vatnajökull icecap covers five active volcanoes, includ-ing Öræfajökull and Bárðarbunga (Figure 1). The Öræ-fajökull volcano (located beneath the southern part of the Vatnajökull icecap) had the most violent historical eruption on record in 1362C.E., devastating nearby ar-eas and killing approximately 300 people [Sharma et al. 2008]. On average, 20 volcanic eruptions occur every century in Iceland; the majority of these eruptions are explosive and basaltic [Gudmundsson et al. 2008]. This explosive pattern is often linked to the interaction be-tween magma and water produced by ice melting. This phreatomagmatism is characteristic of volcanoes cov-ered in ice, as is often the case in Iceland [Gudmunds-son et al. 2008].

The overlying glaciers are thought to have a complex influence on volcanic activity. They could decrease the volcanic activity by the pressure and the weight applied on volcanoes [Sigmundsson 2006]. However, when ice

<footnote><sup>α</sup>Seismology-Gravimetry, Royal Observatory of Belgium, Avenue Circulaire 3, 1180 Brussels, Belgium.</footnote>

<footnote><sup>β</sup>Laboratoire de Glaciologie, Université Libre de Bruxelles, B-1050 Brussels, Belgium.</footnote>

<footnote><sup>γ</sup>Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, IFSTTAR, ISTerre, 38000 Grenoble, France.</footnote>

<footnote><sup>δ</sup>Warning – Forecasting, Icelandic Meteorological Office, IS-150 Reykjavík, Iceland.</footnote>

<footnote>*Corresponding author: sylvain.nowe@seismology.be</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1:** Map of the seismic stations used in this study, with the blue layout being the Vatnajökull icecap and with the major active volcanoes within this icecap. For the volcanoes included in "Others", KVE stands for Kverkfjöll, ESJ stands for Esjufjöll, THO stands for Thor-darhyrna, and EYJ stands for Eyjafjallajökull.</figcaption></figure>

is melting, this pressure decreases, potentially resulting in increased volcanic activity [Sigmundsson 2006]. In the last century, the Vatnajökull icecap has lost approximately 10 % of its mass, with its total volume loss between 1890 and 2003 estimated to be 435 km³. It seems that new volumes of magma are produced in response to this ice thinning [Pagli and Sigmundsson 2008]. This injection of new magma has also been suggested to cause crustal changes and therefore increase seismicity [Pagli and Sigmundsson 2008]. The recent retreat of the Vatnajökull icecap, due to climate change, could cause an increase of 1.4 km³ of produced magma per century [Pagli and Sigmundsson 2008]. Jull and McKenzie [1996] and Tukey [1962] studied the effect of deglaciation on volcanic activity and mantle melting. These studies indicate that 20 to 30 times more magma was produced between 10000 and 4500 years B.P. (corresponding to the unloading of ice) compared to the period since 2900 years B.P. Indeed, Jull and McKenzie [1996] predicted that the removal of 2 km of ice thickness would allow the whole melting column to rise about 0.6 km above its initial depth due to the release of pressure on the underlying mantle. The sensitivity of Vatnajökull icecap to climate change is also of great interest. In temperate conditions (like in Iceland), icecaps are sensitive to small changes in climate. For the Vatnajökull icecap, slight changes in air temperature, for example, could cause drastic changes in its geometry and mass balance [Flowers et al. 2005]. According to simulations, the icecap would lose 31–64 % of its volume and 13–37 % of its area by 2300 with the RCP (Representative Concentration Pathway) 4.5 scenario and 51–94 % of its volume and 24–80 % of its area by 2300 with the RCP 8.5 scenario [Schmidt et al. 2020].

The complex context of Iceland with ice and volca-

noes being collocated and interacting with each other— in addition to the current evolution of climate affecting glacier dynamics and volcanism—are responsible for a very unique seismic context [Flowers et al. 2005; Svattoni et al. 2016; Sigmundsson 2006]. Seismic monitoring methods are very useful for studying Iceland and can complement deformation measurements because both glacial and volcanic activity can generate sub-surface seismic activity. For glacial seismicity, some studies on glaciers and icequakes (glacial earthquakes) have been undertaken, using passive seismic monitoring to classify and characterise these events. In Antarctica, Greenland [Nettles and Ekström 2010], and Switzerland [Canassy et al. 2016], seasonal seismicity has been observed and is associated with glacial processes being enhanced during the summer months. Water content and the transport of water into ice bodies influence glacial processes such as crevasse opening, stick-slip motion, and calving, all of which are directly associated with glacial seismicity. The increase of water content due to ice melting in summer months is therefore expected to increase these processes and therefore the associated induced seismicity [Canassy et al. 2016; Nettles and Ekström 2010]. Nettles and Ekström [2010] focussed on the calving front to characterise the seasonal pattern. Indeed, calving is more effective during summer months and the calving front therefore retreats in summer. Furthermore, calving is believed to be the principal seismic source among glacier processes: Canassy et al. [2016] concluded that glacier motions and seismicity are enhanced between April and November, corresponding with the period of ice melting and already higher velocity of motions of glaciers. In Norway, Köhler et al. [2015] focussed on seasonal glacier seismicity and calving and concluded that a delay in time exists between the increase in seismicity and the actual beginning of the warm period. In that case, calving is linked to sea water temperature (varying slower than air temperature), therefore enhanced calving is expected to lag hotter air temperature by a few days to weeks [Köhler et al. 2015]. Aster and Winberry [2017] characterised all the types of glacial seismicity and the involved processes. These processes are classified into three categories: ice-ice interactions, ice-water interactions, and ice-rock interactions [Aster and Winberry 2017]. Among these processes, surface icequakes, crevasses opening, calving and collisions between icebergs, the stick-slip motion of the glacier, icebergs grounding against the bedrock, and the drainage and transport of subglacial lakes are the most important [Aster and Winberry 2017].

In seismic monitoring, seismic interferometry is of great interest where seismic sources and processes are initially unknown, because it is a passive method that uses every signal that instruments record [Hadziioan-nou et al. 2011; Heirbrant 2006]. Such monitoring techniques are therefore relevant in complex areas like the Vatnajökull region. In recent years, seismic interferom-

<figure><img src="image_4.png" /></figure>
etry has become more important for monitoring Iceland [Donaldson et al. 2019; Konstantinou et al. 2020; Li and Gudmundsson 2020; Svattoni et al. 2016; Woods et al. 2018] and other seismically active areas around the globe [Ballmer et al. 2013; Burtin et al. 2010; Droznin et al. 2015]. Seismic interferometry has been applied in several cases. For example, it was used at Piton de la Fournaise volcano and showed velocity decreases preceding eruptions [Brenguier et al. 2011]. In Iceland, volcanic events have been studied with seismic interferometry, such as the Bárðarbunga–Holuhraun dyke intrusion and eruption in 2014–2015 [Donaldson et al. 2019; Woods et al. 2018]. Woods et al. [2018] focussed on the volcanic tremor during this eruption, on the propagation path of a dyke between the volcano and the actual site of the eruption (Holuhraun), and the associated seismic events [Woods et al. 2018]. Donaldson et al. [2019] studied the Northern Volcanic Zone in Iceland using seismic interferometry to analyse the relative changes in seismic wave velocities ($δv/v$). They used the so-called “cross-component” correlation between pairs of components but using the same station and not pairs of stations. They characterised the changes in $δv/v$ due to the dyke propagation during the 2014–2015 Bárðarbunga–Holuhraun dyke intrusion and eruption as well as a seasonal pattern observed in the yearly $δv/v$ records [Donaldson et al. 2019]. Positive and negative changes are observed after the Bárðarbunga–Holuhraun dyke intrusion and eruption, the positive changes are associated with volumetric compression (areas far from the dyke path) and the negatives changes are associated with volumetric dilatation (areas closer to the dyke path), in line with modeled Coulomb stress changes [Donaldson et al. 2019]. The $δv/v$ seasonal pattern is linked to seasonal physical processes like snow thickness, atmospheric pressure, and ground water level (GWL). Together, these factors represent the total load; however, the GWL also has another influence on velocity changes associated with the pore pressure [Donaldson et al. 2019]. The typical cycle showed $δv/v$ decreases from April to June that continued to decrease but at a slower rate between July and October and finally increased from October to April when it reached the maximum [Donaldson et al. 2019]. More recently, a study located tremor sources using probability density maps and back projection. It has been successfully tested on the Katla volcano (Figure 1) during its unrest in 2011 [Li and Gudmundsson 2020]. Other studies have used seismic interferometry to locate volcanic tremor sources, for example in Hawai'i [Ballmer et al. 2013] and in Kamchatka [Droznin et al. 2015].

Konstantinou et al. [2020] used ambient noise interferometry to study the Gjälp eruption in 1996, which occurred at the so-called fissure located between the Bárðarbunga and Grimsvötn volcanoes on 30 September. Ambient noise interferometry provided evidence for the scenario of a small subglacial eruption prior

to the main event. Furthermore, they showed an increase of $δv/v$ in the caldera of Bárðarbunga from 8 to 28 September, evidence of the pressurisation of the magma chamber. Then, on 28 September, harmonic tremor was observed, suggesting the small subglacial eruption mentioned before [Konstantinou et al. 2020].

Lastly in Iceland, Svattoni et al. [2016] and Jónsdóttir et al. [2009] worked on the region of the Katla volcano. Svattoni et al. [2016] sought to characterise the seasonality of seismicity at the Katla volcano (see Figure 1) and the glacier covering it using seismic interferometry. The seasonal seismicity was not only due to glacier processes but to the interaction between the glacier and volcanic activity. Hydrothermal activity associated with the volcano generates ice melting and therefore induces seismicity. Since ice melting is enhanced in the summer months, more water is available for the hydrothermal system, therefore increasing the observed seismicity [Svattoni et al. 2016]. Jónsdóttir et al. [2009] characterised long-period seismic events in the region of Katla volcano. These events are in fact related to ice movement in a steep outlet glacier and not, as previously thought, to intrusive volcanic activity. Furthermore, these events are consistent in character and magnitude with seasonal changes of the glacier [Jónsdóttir et al. 2009].

Cross-correlation functions have also been used to locate and characterise river-induced seismic signals [Burtin et al. 2010; Möllhoff et al. 2017]. Möllhoff et al. [2017] focussed on the location of rivers related seismic sources in the south-west part of the Vatnajökull icecap. Burtin et al. [2010] used cross-correlation of seismic noise to locate stream segments of the trans-Himalayan Trisuli River, responsible for a large high-frequency seismic noise signal. These segments present a seasonal pattern and are also associated with the most restricted and steepest portions of the river [Burtin et al. 2010].

In this article, we aim to identify, characterise, and locate the main seismic sources and processes in the region of the Vatnajökull icecap by using permanent seismic networks, cross-correlation of ambient seismic noise (continuous data), and two location methods.

## 2 DATA AND METHODS

The data used in this study come from two seismic networks, the VI network, operated by the IMO (Icelandic Meteorological Office) and the Z7 network, operated by the University of Cambridge [White 2010]. The studied period ranges from 2011 to 2019, with all the stations having their own recording periods. For the Z7 network, the stations include three types of sensors from the Güralp company (CMG-6T, CMG-3T, and CMG-3ESP), and all of them are three-component sensors. The frequency bandwidth is from 0.01 Hz to 100 Hz for the CMG-6T sensors and from 0.003 Hz to 50 Hz for the

<figure><img src="image_5.png" /></figure>
CMG-3T and CMG-3ESP sensors. For the VI network, each station is equipped with a three-component seismometer and the sampling rate is 100 Hz.

## 2.1 Seismic interferometry

The method used in this study is seismic interferometry, also called cross-correlation of seismic noise. The idea is that a seismic source generates seismic waves, which are recorded by two seismic stations. The station closer to the source will record the signal before the second one with a delay in time, usually called “lag time”. This lag time can be used to locate the seismic source and to calculate the velocity of seismic waves [Heirbrant 2006]. The main advantage of this method over more traditional ones (P- or S-waves based methods) is its independence to earthquakes, meaning that it does not require any individual earthquake (natural or artificial), although they can be analysed and located if they occur [Hadziioannou et al. 2011].

The daily Cross-Correlations Functions (CCFs) were calculated using MSNoise [Lecocq et al. 2014] at a sampling rate of 20 Hz as a stack (average) of 48 CCFs computed on 1800 s windows. A standard three RMS (Root Mean Square) winsorizing was applied (in the time domain), followed by a spectral whitening (amplitudes set to 1 for every frequencies within the band pass of the filter and sharply decreasing to zero around them (squared half cosine)). Winsorizing is a process taking the extreme values out of the sample of data by replacing them by the nearest unaffected values of the same sample [Tukey 1962], also often referred to as “clipping”.

The CCFs were calculated for three frequency bands: 0.5–1 Hz, 1–5 Hz, and 4–8 Hz. The 0.5–1 Hz band is dominated by the ocean swell and the 4–8 Hz band is dominated by volcano-tectonic events and earthquakes. The 1–4 Hz band therefore offers a continuity of observation characterised by ocean swell as well as volcanism and earthquakes [Nakata et al. 2019].

## 2.2 Location Methods

### 2.2.1 Hyperbolic geometry

For the location of seismic sources using seismic interferometry, the simplest method is to use hyperbolic geometry. Indeed, the source location of any differential time of arrival, or (lag time), recorded by two stations is located anywhere along a hyperbola. In the case of simultaneous arrival (lag time = 0 seconds), the location of the source lies anywhere on a straight line perpendicular to the interstation line. For this method, we assume that the seismic waves detected are predominantly surface waves, therefore restricting our calculations to two dimensions. Our procedure for solving the location involves the following steps. First, geograph-

ical grids are created around the seismic stations and the distance between each point of the grid and each seismic station are calculated after Equation 1:

$$D_A = \sqrt{(x_i - x_A)^2 + (y_i - y_A)^2} \quad (1)$$

where $D_A$ is the distance between station A located at $(x_A, y_A)$ and each $i$th point $(x_i, y_i)$ of the grid. Then, for one pair of stations, the differential distance between every point of the grid and the station pair is calculated (Equation 2). These differential distances can be converted into differential times by dividing the distances by the propagation velocity (Equation 3):

$$\Delta D_{AB} = D_B - D_A \quad (2)$$

$$\Delta T_{AB} = \frac{\Delta D_{AB}}{v} \quad (3)$$

where $\Delta D_{AB}$ is the differential distance for stations A and B, $\Delta T_{AB}$ is the differential orlag time, and v is the propagation velocity in m/s.

There are different ways to solve the non linear relation between time lag and location along an hyperbola, here we simply use the $\Delta T_{AB}$ maps and find the points where $\Delta T_{AB} = \Delta T_{obs}$, the observed lag time. This is done using the contour package of Matplotlib [Hunter 2007].

By taking the lag times observed in the cross-correlations functions and combining several hyperbolas, we can derive the seismic source location from their intersections or convergences. Intersections are assessed in theory with artificial stations and sources but in practice, we assess convergences of hyperbolas (discussed later in the manuscript).

### 2.2.2 Grid-search method

Secondly, we used a grid-search location method—hereafter named “GS method”—originally proposed by Ballmer et al. [2013]. This approach based on ambient seismic noise interferometry allowed them to locate volcanic tremor in Hawai'i. In this method, Ballmer et al. [2013] firstly assume predominant surface wave propagation due to the shallowness of the tremor sources. In our study, oceanic and glacial sources are surface to subsurface sources and therefore this assumption is valid. It is a grid-search method assuming constant propagation velocities. A 0.05° spaced geographical grid is constructed within and around the seismic network with each point being a potential seismic source. The theoretical lag times are then calculated for each of these points and for every pair of stations (Equation 3). Then, in the cross-correlation functions, a window at ±1 s is sliced around the expected lag times. The amplitudes are converted to absolute values, the sliced time-series are summed and the results from all station pairs are finally stacked [Ballmer

<figure><img src="image_6.png" /></figure>
et al. 2013]. These stacked amplitudes give the location probability of seismic sources. The method also outputs the most probable location (highest summed stacked amplitude) defined as the zone with minimum misfit, i.e. the global minimum.

Our base products are daily CCFs, and thus we produce daily probability maps using the GS method. These maps can be stacked for specific filter and propagation velocity over the whole period. This allows us to produce global maps in which each point and probability of location is the mean value for this point over the whole period. Therefore, we can see which region is globally active over a given period of time and decrease the influence of the few days where an unusual source could be dominant for various reasons.

# 3 RESULTS

## 3.1 Propagation velocities

For each of the three frequency bands, we computed the minimum, mean, and maximum propagation velocities from the CCFs sorted by interstation distance (Figure 2). Extracting the lag time corresponding to the maximal amplitude, for each pair of stations and their interstation distance, allows the computation of a linear regression to obtain an estimate of the propagation velocity for each frequency band (Figure 3). Since not all values are close to the linear regression, we can use the 5th and 95th percentiles trends and calculate the associated linear regressions, corresponding to maximum and minimum velocity trends (Table 1). Figure 2 and Figure 3 are examples for frequencies between 0.5 and 1 Hz, but the workflow is the same for other frequency bands.

<table><caption>**Table 1:** Propagation velocities as minimum, average and maximum velocities for each frequency band. Minimum is for the 95th percentile and maximum is for the 5th percentile of the data.</caption><thead><tr><th>Velocity (km s<sup>-1</sup>)</th><th>0.5–1 Hz</th><th>1–4 Hz</th><th>4–8 Hz</th></tr></thead><tbody><tr><th>Minimum</th><td>1.39</td><td>1.47</td><td>2.12</td></tr><tr><th>Average</th><td>1.89</td><td>2.79</td><td>2.79</td></tr><tr><th>Maximum</th><td>2.86</td><td>3.78</td><td>3.92</td></tr></tbody></table>

## 3.2 Identification of different sources

Using correlograms between all station pairs from 2011 to 2019, we identify three major types of seismic signatures: permanent, seasonal, and episodic, described here.

<figure><img src="image_7.png" /><figcaption>**Figure 2:** Distance plot for frequencies between 0.5 and 1 Hz, showing lag times on horizontal axis and interstation distances on vertical axis. Solid lines are the theoretical velocities and the coloured points represent the maximum correlation for each pair of stations.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 3:** Absolute values of lag times against interstation distances. The slope is equivalent to $1/v$. The 5th and 95th percentiles stand for maximum and minimum velocity trends, respectively.</figcaption></figure>

### 3.2.1 Permanent source

Firstly, we observe a very constant signature with the same pattern over the years on the correlograms between 0.5 and 1 Hz. Although also usable, the correlograms from 1 to 4 Hz did not add any information and were less clear. Figure 4A is a typical cross-correlation function showing this permanent seismic signature (blue outline) as well as an episodic one (red outline) described later in the manuscript.

Figure 4A shows that most of the seismic signals arrive at the station IEY approximately 20 seconds before they arrive at station JOK, indicating a seismic source

<figure><img src="image_9.png" /></figure>
<figure><img src="image_10.png" /></figure>

**Figure 4:** [A] Plot of the Cross-Correlation Function versus time (CCFtime) for the stations IEY and JOK and for frequencies between 0.5 and 1 Hz. The blue outline shows the permanent signature and the red outline shows the episodic signature. [B] Plot of the Cross-Correlation Function versus time (CCFtime) for the stations JOK and KAL for frequencies between 4 and 8 Hz. The orange outline shows the seasonal signature.

located closer to IEY than JOK.

### 3.2.2 Seasonal source

Secondly, we observe a seismic signature with a very clear seasonal pattern. The seismicity is visible during the summer months (usually May to September) and for frequency bands between 1 and 8 Hz (orange outline in Figure 4B). It was identified by observing either the 1–4 Hz or the 4–8 Hz cross-correlation functions. Figure 4B shows a typical example of this seasonality.

Here, the seismic waves arrive at station KAL 17 seconds before they arrive at station JOK. It is consistent with a seasonal source located closer to KAL.

### 3.2.3 Episodic sources

Finally, we observed a seismic signature lasting from mid-August 2014 to end of February 2015 (red outline in Figure 4A). It is visible at every frequency between 0.5 and 8 Hz and is the most energetic source since it is dominant during this period (its amplitude is the largest in the correlograms during this period), hiding every other seismic signature. This seismic signature is characterised by discontinuous patterns in the cross-correlation functions. This signature presented in Figure 4A for pair IEY-JOK, is associated with the Bárðarbunga–Holuhraun dyke intrusion and eruption

<figure><img src="image_11.png" /></figure>

**Figure 5:** Map of hyperbolas for a lag time of 12.2 s and an interstation distance of 31.3 km, showing the differences between two propagation velocities, 1.39 and 1.89 km s⁻¹.

between 16 August 2014 and 27 February 2015 [Woods et al. 2018]. Figure 4A also exhibits an issue with one of the sensor clocks, resulting in an apparent “slide” of the CCF at the end of the eruption period.

## 3.3 Location of seismic sources

### 3.3.1 Hyperbolic method

In Figure 5, hyperbolas for two different propagation velocities are shown for the station pair BHAF-KSK. Three velocities were used, namely the minimum, average, and maximum velocities calculated for the frequency band 0.5–1 Hz (Table 1). The maximum velocity (2.86 km s⁻¹) is not plotted because for such a high velocity, no hyperbolas exists for this specific pair of stations based on the interstation distance and differential time. We discuss this aspect later in the discussion section.

We then show the convergences of hyperbolas for the permanent and the seasonal seismic sources. When calculating hyperbolas for a given set of station pairs, it is important to use the same propagation velocity. The choice between the minimum, average, and maximum velocity trends is based on which velocity exists for these pairs, as mentioned previously.

Figure 6 shows that for the permanent seismic source and for a velocity of 1.39 km s⁻¹ (minimum velocity trend for frequencies between 0.5 and 1 Hz), hyperbolas converge on a source located perpendicular to the shoreline and approximately 20 km in the ocean. Figure 7 shows that for the seasonal seismic source and for a velocity of 2.79 km s⁻¹ (average velocity trend for frequencies between 1 and 8 Hz), hyperbolas converge on a source located a few km south-west of the glacier tongue in the area of glacial rivers. We assessed convergences of hyperbolas and not intersections—as

<figure><img src="image_12.png" /></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 6:** Map showing three hyperbolas for stations BHAF, BJK and KSK for a propagation velocity of 1.39 km s⁻¹. Lag times are: 36.65 s for BHAF-BJK, 12.2 s for BHAF-KSK, and -22.65 s for BJK-KSK.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 7:** Map showing four hyperbolas for stations FAG, FAL, HUS, IEY, KAL and SKAF for a propagation velocity of 2.79 km s⁻¹. Lag times are: -20.05 s for FAG-IEY, -19.3 s for FAL-IEY, -8 s for HUS-IEY, and 10.55 s for KAL-SKAF.</figcaption></figure>

mentioned before—because we used trends of velocity, meaning that light differences of propagation velocities between pairs of stations prevent assessing an exact intersection.

### 3.3.2 GS method

The first three panels (A, B and C) in Figure 8 are examples of the daily solutions with the colour scale showing the probability of locations (summed stacked amplitudes of the cross-correlation functions) for the seismic sources. Figure 8A displays the permanent oceanic sources locations. We can see that the azimuth of the source is very well constrained but not the distance. Figure 8B represents the seasonal glacial-related sources that are well constrained. Figure 8C focuses on the Bárðarbunga–Holuhraun eruptive event in 2014–2015. We can see that the northern part of the Vatnajökull icecap is seismically active. However, the oceanic signature (south of the Vatnajökull icecap) is visible as well as it persists throughout the year. Since this method uses all pairs of stations, we used the average trend of propagation velocity for both frequency bands; 1.89 km s⁻¹ between 0.5 and 1 Hz and 2.79 km s⁻¹ between 1 and 8 Hz.

We stacked the daily maps to obtain maps in which each point and probability of location is the mean value for this point over the whole period (depending on the source). Figure 8D–F shows maps of the stacked probabilities of location assessed with the GS method, combined with previous hyperbolas to compare results.

In Figure 8D, the azimuthal direction of points with highest probability of location is consistent with the hyperbolas. However, the distance between the source and the shoreline is not well constrained. We used different propagation velocities for hyperbolas and the GS method because the latter is based on all the station

pairs, it is therefore important to use the average propagation velocity (see Table 1). In Figure 8E, the points with the highest probability of location exactly coincide with the convergence area of hyperbolas. This seismic source is probably linked to glacial rivers and fluxes in the summer months. This location is validated by a study on back-azimuth projections and slowness curves combined with hydrologic data [Möllhoff et al. 2017] in which the authors located seasonal seismic sources in the Hverfisfljot river at the location of a waterfall (Figure 9). In Figure 8F, we zoomed in to the area of the Bárðarbunga–Holuhraun eruption. The dyke path and the eruption fissures [Woods et al. 2018] are included in the map to compare with the probabilities of location assessed with the GS method. We can see that the path of the dyke and the two eruption fissures are in an area of high probability of location according to the GS method. However, the highest probabilities (yellow on Figure 8F), are located a few kilometers north of the dyke path. This is due to the disposition of seismic stations around this eruptive event and the presumed velocity. We discuss this aspect later in the discussion section.

# 4 DISCUSSION

## 4.1 Location methods

When calculating hyperbolas some velocity trends do not match the observed lag times. This occurs for the highest velocity, which we relate to the picking of wiggles originating from outside of the end-fire lobes for specific station pairs, due to a non-homogeneous distribution of sources around the given station pair. Therefore, when combining hyperbolas to assess a location, it is important to use one of the three velocities that ex-

<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /></figure>

**Figure 8:** [A] to [C] Daily geographical grids assessed with the GS method showing the seismic stations and the probability of source location (summed stacked amplitudes) represented by the colour scale. [A] 25th May 2015 for frequencies between 0.5 and 1.0 Hz and a propagation velocity of 1.89 km s⁻¹. [B] 16th May 2014 for frequencies between 4.0 and 8.0 Hz and for a propagation velocity of 2.79 km s⁻¹. [C] 14th October 2014 for frequencies between 0.5 and 1.0 Hz and a propagation velocity of 1.89 km s⁻¹. [D] to [F] Maps of the stacked daily probabilities of location assessed with the GS method. Seismic stations are not represented for clarity. [D] Map of the Oceanic source for frequencies between 0.5 and 1.0 Hz and from 1st May to 30th July 2015. The propagation velocity is set at 1.89 km s⁻¹. Hyperbolas are the same as in Figure 6. [E] Map of the glacial source for frequencies between 4.0 and 8.0 Hz and from 10th May to 15th June 2014. The propagation velocity is set at 2.79 km s⁻¹. Hyperbolas are the same as in Figure 7. [F] Map for the Bárðarbunga–Holuhraun eruption for frequencies between 0.5 and 1 Hz and from 15th August 2014 to 1st February 2015. The propagation velocity is set at 1.89 km s⁻¹. The path of the dyke and the eruption fissures are shown as well and come from Woods et al. [2018].

<figure><img src="image_17.png" /></figure>
ist for every used pair. Furthermore, when this velocity is chosen, the theoretical intersection between hyperbolas is not necessarily assessed and convergences are assessed instead, as seen in Figure 6 and Figure 7. This can be due to slight variations in propagation velocities from one pair to another or to uncertainties of the lag time measurement. Indeed, the hyperbola method relies on a single lag time value for the maximum of the correlation, and a ±1 s uncertainty around that value would change the shape of the hyperbola. This limitation is overcome when using the GS method, as it uses the full CCF waveform information and not a single Δt value. The method uses the information located ±1 s (constant) around the predicted differential time, effectively allowing for small velocity variations to occur within the sampled medium.

Both methods are, by design, not able to provide accurate locations outside of the network, as the location of hyperbolas branches overlap and do not cross with large angles as is the case for solutions within the network. The azimuth, in turn, is correctly determined for these sources.

It is important to note that by stacking the amplitudes to obtain the maps of stacked probability of location (Figure 8D–F), we suppress some more episodic and irregular sources that could have been studied. However, our goal in this study was to locate and characterise the main seismic sources in the region.

## 4.2 Seismic sources

### 4.2.1 Permanent, oceanic source

The location of the oceanic source outside of the network is not precise in terms of location but their azimuth relative to the network is well constrained. We

<figure><img src="image_18.png" /></figure>

**Figure 9:** Map showing the stacked probabilities of location at 2.79 km s⁻¹ for the seasonal source as well as the path of the Hverfisfljot river and the associated waterfall.

used one seismic station from the Faroe Islands (SOFL) to compute the cross-correlation functions between this new station and the Icelandic seismic network (see Figure 8A). If a strong oceanic source lies between Iceland and the Faroe Islands or Scotland, the CCFs should contain coherent arrivals. However, the cross-correlation functions did not show any coherent signals between the Faroe Islands and Icelandic stations. We conclude that the oceanic source is located much closer to Iceland, therefore too far from the Faroe Islands and Scotland to be recorded there.

The source seems to be associated with very unique coastal settings. This source (Figure 8D), located east of the Öræfajökull volcano, seems to be associated with an old glacier bed currently forming an underwater canyon of nearly 170 m deep below sea level [Angel Ruiz Angulo and Halldór Björnsson, pers. comm., and NCEI (National Centers for Environmental Information) bathymetric data from NOAA (National Oceanic and Atmospheric Administration)]. This focusses waves in this depression making the coupling between these waves and the ocean floor more efficient. Combined with the parabolic form of the shoreline focussing waves by reflection, it can explain the observed seismic noise.

### 4.2.2 Seasonal, glacial-related source

The location of the seasonal source—inside the network—is precise, and located approximately 10 km SSW of the Síðujökull glacier's tongue in the western margin of the Vatnajökull icecap (Figure 8E). This area is characterised by the Hverfisfljot glacial river starting at the margin of the glacier [Möllhoff et al. 2017] and flowing southwards towards the ocean (Figure 9).

In glacial seismology, many processes can generate seasonal seismicity. Processes giving rise to glacial seismicity can usually be divided into three categories: ice-ice interactions, ice-water interactions, and ice-rock interactions [Aster and Winberry 2017]. The collisions between icebergs, calving, the basal stick-slip motion of a glacier (known to generate important seismicity through the friction between the glacier and the bedrock), and the drainage and transport of subglacial lakes are examples of these potential processes [Aster and Winberry 2017]. However, in our study, these processes cannot be directly linked to the source we observed since our location is outside of the glacier tongue.

The observed seasonal increase of glacier-related seismicity during the summer is validated by several studies. It has been observed in Switzerland where glacier motions are enhanced in summer months and a corresponding burst in seismicity is observed [Canassy et al. 2016]. In Antarctica and Greenland, seasonal seismicity is also present and linked to the added water amount enhancing glacier motions and glacial processes such as crevasse opening and stick-slip motion,

<figure><img src="image_19.png" /></figure>
therefore increasing seismicity [Nettles and Ekström 2010]. Finally, in Norway, the burst in glacial seismicity due to calving occurs later than the rising in air temperature because of the sea water [Köhler et al. 2015]. To summarise, processes usually involved with seasonal seismicity are the surface icequakes, the opening of crevasses, basal stick-slip motion, and transport and drainage of subglacial lakes, each of which are believed to be enhanced in the summer months when ice melts.

Since we located the seismic source outside of the glacier, the enhanced flow of the glacial rivers during the warm period of the year seems to be the only possible explanation. This generates seismic waves against the bedrock of the rivers, especially in the area of a waterfall. Usually, seismic sources associated with rivers would not be spatially located but distributed along the path of the river depending on the river flow (higher seismicity where the flow is higher) [Einarsson 2018]. However, along the route of glacial rivers some enhanced seismicity can be generated, for example at waterfalls or at the convergence between several glacial rivers. Our location is validated by Möllhoff et al. [2017], who located seasonal seismic sources using seismic analysis, back-azimuth projections, and slowness curves combined with hydrologic data. Möllhoff et al. [2017] located seismic sources in the Hverfisfljot river, in particular at the associated waterfall. These locations are exactly the same as the area we located in our study (Figure 9). The advantage of our method over that of Möllhoff et al. [2017] is that we did not need to install a new array of sensors to reach the same results. By only using permanent seismic networks and continuous seismic noise data, we were able to identify this seismicity and locate the area of the Hverfisfljot river and the associated waterfall. Locating seismic noise originating from a waterfall with such accuracy by using permanent networks has never been done before. Nevertheless, at this point, it is not possible to completely eliminate the indirect influence of proper glacial processes on our located source since the usage of another velocity assumption could change the source location, closer to or further from the glacier's tongue. We plan to further explore this in the future by installing an array of seismometers around the glacier tongues and down the path of glacial rivers.

### 4.2.3 Episodic volcanic events and unrest

The Bárðarbunga–Holuhraun dyke intrusion and eruption of 2014–2015 at Holuhraun is the only volcanic event observed, characterised, and analysed in this work. It was the largest eruption in Iceland for 230 years and was associated with more than 30,000 volcano-tectonic earthquakes [Woods et al. 2018]. Tremor was also intense during the event and was studied and located in a previous study using cross-correlation as well as the GS method [Woods et

al. 2018]. It is important to mention that periods of weak tremor were difficult to locate and isolate from other seismic events. These seismic events were linked to dyke propagation between the volcano and the site of the actual eruption (Holuhraun) as well as to the collapse of the caldera [Woods et al. 2018]. Three aspects need to be discussed regarding this event and the associated seismic signatures in our study.

The first aspect is the broad lag times observed in every CCF time plot including this eruption. The lag times spread from the negatives to the positives as can be seen in Figure 4A for example, from -45 seconds to +45 seconds. There are many potential eruption processes that could be associated with different seismic velocities. Each of these processes can be defined as different sources, even if located in the same area. The correlation itself can generate positive and negative lag times when P- and S-waves (having different velocities) correlate together. Lastly, the distribution of the seismic sources is important. In the case of the Bárðarbunga–Holuhraun dyke intrusion and eruption, a dyke propagated from the volcano to the site of the eruption, Holuhraun, with a total spread of more than 40 km. This motion was discontinuous but nonetheless results in a large distribution of seismic sources all along the path of the dyke [Ágústsdóttir et al. 2016]. Different source locations mean different distances between sources and seismic stations, eventually resulting in different travel times and lag times.

The second aspect is the difference between the time period of the eruption seen in the cross-correlation functions in comparison with the real time period found in the literature, from 16 August to 27 February [Woods et al. 2018]. The former is usually shorter than the latter and is observed between 30 September 2014 and 20 January 2015. The first factor to explain this variation is the normalisation of the energy of the seismic traces performed automatically by MSNoise when computing the cross-correlation functions. It means that only the most energetic seismic signatures are represented in the plots and are dominant in the CCFs even if they last only for a few hours. This is a limitation due to the usage of the permanent seismic network and parameters suited for general observations, directly linked to our workflow and objectives. Therefore, the shorter period of eruption in the cross-correlation functions suggests that at the very beginning and at the very end of this event, it is less energetic than the other seismic signatures. For example, in Figure 4A, the eruption is less energetic than the oceanic seismic source at the beginning and the end of the eruptive period, and therefore is not visible.

The last aspect to discuss is the locations assessed with the GS method for this volcanic event in comparison with the study of Woods et al. [2018]. In their study, Woods et al. [2018] used waveform cross-correlation to analyse and characterise seismicity (mainly long-period earthquakes and tremor) during the Bárðar-

<figure><img src="image_20.png" /></figure>
bunga-Holuhraun dyke intrusion and eruption of 2014–2015. For this purpose they used 72 seismometers from the University of Cambridge seismic network, 14 stations form the IMO network as well as one British Geological Survey station and one University College Dublin station. For the waveform cross-correlation and the location of seismic tremor, two frequency bands were used: 0.01 to 8 Hz and 0.5 to 2 Hz. The location was performed based on the approach of Ballmer et al. [2013] and the propagation velocity was set to 1.2 km s⁻¹ [Woods et al. 2018]. In comparison with our study, the general approach is the same but some parameters differ. Our frequency band is narrower (0.5 to 1 Hz) and our propagation velocity is higher (1.89 km s⁻¹). Our velocity calculation was based on all the seismic signatures found between 0.5 and 1 Hz, where the oceanic source is dominant, meaning that this velocity could be less suited to this volcanic event. We believe that this difference in propagation velocity can generate changes of a few km in the resulting locations because a velocity of 1.2 km s⁻¹ would induce locations closer to the network compared to a velocity of 1.89 km s⁻¹. Another important point is the geometry of the seismic network. Indeed, Woods et al. [2018] used seismic instruments located all around the dyke path and eruption fissures, resulting in a best location when using a GS-like method. In our case, all the instruments are located south of the volcanic events. As seen for the oceanic source, this aspect is responsible for good azimuthal determination but poor accuracy and precision. However, by using fewer specific parameters (seismic network geometry, propagation velocity, and frequency ranges), we assessed locations located only a few km north relative to the results of Woods et al. [2018].

Unrest at Öræfajökull volcano occurred in 2017–2018 but could not be seen in the CCFs. This unrest was characterised by geothermal activity and earthquakes. Airborne surveillance has also shown circular cauldrons due to the melting of ice on the volcano. All this could suggest magma up-flow in the conduits and therefore prelude an eruption in the next months or years [Trausti 2017]. The reason this unrest was not visible in this work could be due to the method itself. It does not mean that seismic interferometry is not adequate but that our parameters are not suitable for that case study. The normalisation of the energy of the seismic traces has potentially hidden the seismic signatures of the unrest, meaning that they are not energetic enough in comparison with other seismic sources, such as glacial or oceanic ones. The unrest would be visible if we change parameters and if we optimise the method.

# 5 CONCLUSIONS AND PERSPECTIVES

This work was part of an Icelandic project called IS-Noise, aiming at evaluating the potential of ambient

noise-based methods to monitor different targets in Iceland (volcanoes, fault zones, and geothermal areas). The main objective was to reach a good understanding of the background of seismic events and seismic sources in the region of the Vatnajökull icecap using seismic interferometry with a specific set of parameters. Three frequency bands were investigated from 0.5 to 8 Hz and cross-correlation functions for every pair of stations provided by two networks between 2009 and 2019 were computed. We were able to characterise and locate three primary natural seismic sources in the region of the Vatnajökull icecap.

The oceanic seismic source, generating permanent and clear seismic signatures, is associated with a very specific coastal context. The glacial source generating seasonal seismicity has also been identified and located. It seems to be linked to glacier rivers flow, the Hverfisfljot river, and the associated waterfall, enhanced by the ice melting during summer months. Finally, the Bárðarbunga-Holuhraun dyke intrusion and eruption in 2014–2015 was identified. We were able to highlight the complexity of this event using seismic interferometry.

For the location of seismic sources, we calculated propagation velocities and used the geometry of hyperbolas as well as the GS method. The results from both methods are similar and we discussed the slight differences occurring. For the oceanic source, azimuthal direction is similar but the distance is not well constrained by the GS method. However, we concluded that the source is located close to the shoreline. For the glacial source, locations from both location methods are the same and are validated by a previous study, confirming this area as being a principal seismic source in the region. The other seismically active areas around the Hverfisfljot river could be associated with glacial processes as well as all kinds of seismic noise but we did not analyse this aspect and focussed on the main seismic source. This will be the focus of future research.

This study shows that seismic interferometry is efficient to characterise and locate seismic processes and sources in geologically and physically active regions such as Iceland. Characterising and locating different seismic sources in the same region with the same set of instruments, method, and workflow has never—to the authors’ knowledge—been done before. Furthermore, this study shows that passive methods based on ambient seismic noise can be more suitable than other seismic methods (from the instrumental cost to the continuous recording), depending on the case study.

Further work is required to investigate the physical processes generating the observed seismic signatures in the region of Vatnajökull and to discriminate between glacial and volcanic sources of seismic noise. Indeed, we know that these two categories of sources can act together. It will be worth working on the regional scale of the Vatnajökull icecap to understand how glaciers and volcanoes are linked using the same method of seismic

<figure><img src="image_21.png" /></figure>
interferometry. Lastly, going sub-daily in the workflow to avoid the dominant source associated with the daily average is necessary to optimise the resolution.

## ACKNOWLEDGEMENTS

We thank the developers of Obspy [Krischer et al. 2015], Python [Van Rossum and Drake 2009], MSNoise [Lecocq et al. 2014] and QGIS [QGIS Development Team 2021] for the software used in this study. We also thank the FNRS (Fonds de la Recherche Scientifique) for funding the PhD project associated with this study (contract number FC 36411) as well as the IS-NOISE project (https://is-noise.earth/), which is supported by the Icelandic Research Fund, Rannis (https://www.rannis.is/) (contract number 185209-051), for providing the data and funding the field work in Iceland that allowed us to discuss with experts at the IMO. We sincerely thank these same experts from the IMO, Angel Ruiz Angulo, Halldor Bjornsson, Bergur Einarsson, and Yesim Cubuk Sabuncu.

## AUTHOR CONTRIBUTIONS

Thomas Lecocq and Corentin Caudron provided the MSNoise method and workflow to create the cross-correlation functions. Grid-search location method was provided by Thomas Lecocq. Corentin Caudron and Kristín Jónsdóttir provided the seismic data from the two seismic networks. Kristín Jónsdóttir allowed us to communicate with glaciologists and oceanographers during a meeting at IMO (Iceland) to discuss our results. Thomas Lecocq, Corentin Caudron, and Frank Pattyn supervised the research.

## DATA AVAILABILITY

The seismic data were obtained from the IMO via Kristín Jónsdóttir and from the University of Cambridge via Corentin Caudron. The MSNoise method and workflow used in this study is free and in open access, as well as all the documentation. It can be found at http://msnoise.org/doc/. The stations BHAF, KALF, LAUF, SKAF, SVIN, and SYLG from the Z7 network from 2011 to 2015 are freely available on IRIS. Data from January 2016 will be available in 2022/2023 with the network code 8K. Data from other stations should be publicly available on ORFEUS in the future.

## COPYRIGHT NOTICE

© The Author(s) 2021. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

## REFERENCES

Ágústsdóttir, T., J. Woods, T. Greenfield, R. G. Green, R. S. White, T. Winder, B. Brandsdóttir, S. Steinthórs-son, and H. Soosalu (2016). "Strike-slip faulting during the 2014 Bárðarbunga-Holuhraun dike intrusion, central Iceland". *Geophysical Research Letters* 43 (4), pp. 1495-1503. doi: 10.1002/2015GL067423.

Aster, R. C. and J. P. Winberry (2017). "Glacial seismology". *Reports on Progress in Physics* 80 (12), p. 39. doi: 10.1088/1361-6633/aa8473.

Ballmer, S., C. J. Wolfe, P. G. Okubo, M. M. Haney, and C. H. Thurber (2013). "Ambient seismic noise interferometry in Hawai'i reveals long-range observability of volcanic tremor". *Geophysical Journal International* 194 (1), pp. 512-523. doi: 10.1093/gji/ggt112.

Brenguier, F., D. Clarke, Y. Aoki, N. M. Shapiro, M. Campillo, and V. Ferrazzini (2011). "Monitoring volcanoes using seismic noise correlations". *Comptes rendus - Géoscience* 343 (8), pp. 633-638. doi: 10.1016/j.crte.2010.12.010.

Burtin, A., J. Vergne, L. Rivera, and P. Dubernet (2010). "Location of river-induced seismic signal from noise correlation functions". *Geophysical Journal International* 182 (3), pp. 1161-1173. doi: 10.1111/j.1365-246X.2010.04701.x.

Canassy, P. D., C. Röösli, and F. Walter (2016). "Seasonal variations of glacier seismicity at the tongue of Rhonegletscher (Switzerland) with a focus on basal icequakes". *Journal of Glaciology* 62 (231), pp. 18-30. doi: 10.1017/jog.2016.3.

Donaldson, C., T. Winder, C. Caudron, and R. S. White (2019). "Crustal seismic velocity responds to a magnetic intrusion and seasonal loading in Iceland's Northern Volcanic Zone". *Science Advances* 5 (11), eaax6642. doi: 10.1126/sciadv.aax6642.

Droznin, D. V., N. M. Shapiro, S. Y. Droznina, S. L. Senyukov, V. N. Chebrov, and E. I. Gordeev (2015). "Detecting and locating volcanic tremors on the Klyuchevskoy group of volcanoes (Kamchatka) based on correlations of continuous seismic records". *Geophysical Journal International* 203 (2), pp. 1001-1010. doi: 10.1093/gji/ggv342.

Einarsson, B. (2018). *Subglacial hydrology of the Icelandic ice caps: Outburst floods and ice dynamics*. Bergur Einarsson. ISBN: 978-9935-9306-9-9.

Flowers, G. E., S. J. Marshall, H. Björnsson, and G. K. C. Clarke (2005). "Sensitivity of Vatnajökull ice cap hydrology and dynamics to climate warming over the next 2 centuries". *Journal of Geophysical Research: Earth Surface* 110 (F2). doi: 10.1029/2004JF000200.

Gudmundsson, M., G. Larsen, A. Hoskuldsson, and Á. Gylfason (2008). "Volcanic hazards in Iceland". *Jökull*. 58, pp. 251-268.

Hadzioannou, C., E. Larose, A. Baig, P. Roux, and M. Campillo (2011). "Improving temporal resolution in ambient noise monitoring of seismic wave speed".

<figure><img src="image_22.png" /></figure>
Journal of Geophysical Research: Solid Earth 116 (B7). DOI: 10.1029/2011JB008200.

Heirbrant, K. (2006). "Identification du bruit sismique ambiant sur le réseau de Parkfield en Californie". MA thesis. Laboratoire de Géophysique Interne et de Tectonophysique, Grenoble.

Hunter, J. D. (2007). "Matplotlib: A 2D Graphics Environment". Computing in Science & Engineering 9 (03), pp. 90-95. DOI: 10.1109/MCSE.2007.55.

Jónsdóttir, K., R. Roberts, V. Pohjola, B. Lund, Z. H. Shomali, A. Tryggvason, and R. Böðvarsson (2009). "Glacial long period seismic events at Katla volcano, Iceland". Geophysical Research Letters 36 (11). DOI: https://doi.org/10.1029/2009GL038234.

Jull, M. and D. McKenzie (1996). "The effect of deglaciation on mantle melting beneath Iceland". Journal of Geophysical Research: Solid Earth 101 (B10), pp. 21815-21828. DOI: 10.1029/96JB01308.

Köhler, A., C. Nuth, J. Schweitzer, C. Weidle, and S. J. Gibbons (2015). "Regional passive seismic monitoring reveals dynamic glacier activity on Spitsbergen, Svalbard". Polar Research 34 (1), p. 26178. DOI: 10.3402/polar.v34.26178.

Konstantinou, K. I., I. W. Utami, D. Giannopoulos, and E. Sokos (2020). "A Reappraisal of Seismicity Recorded During the 1996 Gjálp Eruption, Iceland, in Light of the 2014-2015 Bárðarbunga-Holuhraun Lateral Dike Intrusion". Pure and Applied Geophysics 177 (6), pp. 2579-2595. DOI: 10.1007/s00024-019-02387-x.

Krischer, L., T. Megies, R. Barsch, M. Beyreuther, T. Lecocq, C. Caudron, and J. Wassermann (2015). "ObsPy: a bridge for seismology into the scientific Python ecosystem". Computational Science & Discovery 8 (1). DOI: 10.1088/1749-4699/8/1/014003.

Lecocq, T., C. Caudron, and F. Brenguier (2014). "MSNoise, a Python Package for Monitoring Seismic Velocity Changes Using Ambient Seismic Noise". Seismological Research Letters 85 (3), pp. 715-726. DOI: 10.1785/0220130073.

Li, K. L. and O. Gudmundsson (2020). "A Probabilistic Tremor Location Method". Geophysical Research Letters 47 (4), e2019GL085538. DOI: 10.1029/2019GL085538.

Martin, E., J. L. Paquette, V. Bosse, G. Ruffet, M. Tiepolo, and O. Sigmarsson (2011). "Geodynamics of rift-plume interaction in Iceland as constrained by new 40Ar/39Ar and in situ U-Pb zircon ages". Earth and Planetary Science Letters 311 (1), pp. 28-38. DOI: 10.1016/j.epsl.2011.08.036.

Möllhoff, M., E. Eibl, C. Bean, and K. Vogfjörd (2017). "Remote monitoring of water-flow induced seismic noise". Ohlstadt, Germany.

Nakata, N., L. Gualtieri, and A. Fichtner (2019). Seismic Ambient Noise. Cambridge University Press. ISBN: 978-1-108-41708-2.

Nettles, M. and G. Ekström (2010). "Glacial Earthquakes in Greenland and Antarctica". Annual Review of Earth and Planetary Sciences 38 (1), pp. 467-491. DOI: 10.1146/annurev-earth-040809-152414.

Pagli, C. and F. Sigmundsson (2008). "Will present day glacier retreat increase volcanic activity? Stress induced by recent glacier retreat and its effect on magnetism at the Vatnajökull ice cap, Iceland". Geophysical Research Letters 35 (9). DOI: 10.1029/2008GL033510.

QGIS Development Team (2021). QGIS Geographic Information System. QGIS Association.

Schmidt, L. S., G. Ådalgeirsdóttir, F. Pálsson, P. L. Langen, S. Gudmundsson, and H. Björnsson (2020). "Dynamic simulations of Vatnajökull ice cap from 1980 to 2300". Journal of Glaciology 66 (255), pp. 97-112. DOI: 10.1017/jog.2019.90.

Sgattoni, G., Z. Jeddi, Ó. Gudmundsson, P. Einarsson, A. Tryggvason, B. Lund, and F. Lucchi (2016). "Long-period seismic events with strikingly regular temporal patterns on Katla volcano's south flank (Iceland)". Journal of Volcanology and Geothermal Research 324, pp. 28-40. DOI: 10.1016/j.jvolgeores.2016.05.017.

Sharma, K., S. Self, S. Blake, T. Thordarson, and G. Larsen (2008). "The AD 1362 Öræfajökull eruption, S.E. Iceland: Physical volcanology and volatile release". Journal of Volcanology and Geothermal Research 178 (4), pp. 719-739. DOI: 10.1016/j.jvolgeores.2008.08.003.

Sigmundsson, F. (2006). Iceland Geodynamics: Crustal Deformation and Divergent Plate Tectonics. Springer Praxis Books. Berlin, Heidelberg: Springer Berlin Heidelberg. ISBN: 978-3-540-24165-2.

Thordarson, T. and G. Larsen (2007). "Volcanism in Iceland in historical time: Volcano types, eruption styles and eruptive history". Journal of Geodynamics. Hotspot Iceland 43 (1), pp. 118-152. DOI: 10.1016/j.jog.2006.09.005.

Trausti, A. (2017). Unrest in Mt. Öræfajökull. URL: https://lavacentre.is/620-2/ (visited on 06/16/2019).

Tukey, J. W. (1962). "The Future of Data Analysis". Annals of Mathematical Statistics 33 (1), pp. 1-67. DOI: 10.1214/aoms/1177704711.

Van Rossum, G. and F. L. Drake (2009). Python 3 Reference Manual. Scotts Valley, CA: CreateSpace. ISBN: 1441412697.

White, R. (2010). Northern Volcanic Zone. DOI: 10.7914/SN/Z7_2010. [Data set].

Woods, J., C. Donaldson, R. S. White, C. Caudron, B. Brandsdóttir, T. S. Hudson, and T. Ágústsdóttir (2018). "Long-period seismicity reveals magma pathways above a laterally propagating dyke during the 2014-15 Bárðarbunga rifting event, Iceland". Earth and Planetary Science Letters 490, pp. 216-229. DOI: 10.1016/j.epsl.2018.03.020.

<figure><img src="image_23.png" /></figure>