Entirely covered by the Vatnajökull ice cap, Grímsvötn is among Iceland’s largest and most hazardous volcanoes. Here we demonstrate that fiber‐optic sensing technology deployed on a natural floating ice resonator can detect volcanic tremor at the level of few nanostrain/s, thereby enabling a new mode of subglacial volcano monitoring under harsh conditions. A 12.5 km long fiber‐optic cable deployed on Grímsvötn in May 2021 revealed a high level of local earthquake activity, superimposed onto nearly monochromatic oscillations of the caldera. High data quality combined with dense spatial sampling identify these oscillations as flexural gravity wave resonance of the ice sheet that floats atop a subglacial lake. Although being affected by the ambient wavefield, the time–frequency characteristics of observed caldera resonance require the presence of an additional persistent driving force with temporal variations over several days, that is most plausibly explained in terms of low‐frequency volcanic tremor. In addition to demonstrating the logistical feasibility of installing a large, high‐quality fiber‐optic sensing network in a subarctic environment, our experiment shows that ice sheet resonance may act as a natural amplifier of otherwise undetectable (volcanic) signals. This suggests that similar resonators might be used in a targeted fashion to improve monitoring of ice‐covered volcanic systems.

Grímsvötn is among Iceland’s largest and, on a decadal time scale, most active volcanoes. It is located within a dense cluster of major volcanoes beneath the Vatnajökull ice cap where the Iceland hotspot is centered (Wolfe et al., 1997). Grímsvötn ranks among the most productive volcanic systems on Earth, being responsible for the 1783 Laki eruption, the largest ever witnessed by mankind, which affected climate on a global scale (Thordarson and Self, 2003). With a caldera of ∼10 km in diameter, major eruptions occurred on average every 10 yr during the past century (Global Volcanism Program, 2012). The volcanic system is marked by the complex interaction of geothermal heating and ice melting, which leads to the formation of a subglacial lake. Occasional outburst floods cause inundations of the coastal plains and destruction of infrastructure (Gudmundsson et al., 1997). Grímsvötn’s caldera eruptions are explosive, as the interaction of basaltic magma, gases, and meltwater produces ash plumes of tens of kilometers height, which affect human lives, agriculture, livestock, and aviation.

The last major eruption started on 21 May 2011 and lasted for 8 days. The ash cloud, peaking at 20 km altitude, led to the cancellation of 900 flights and the closure of Iceland’s main airport, Keflavik (Wunderman, 2011). In December 2021, Grímsvötn’s showed increased activity, including elevated seismic activity, ice sheet subsidence by nearly 80 m, and drainage of more than 0.8  km3 of water from the subglacial lake (Global Volcanism Program, 2021).

The ongoing activity of Grímsvötn and the size of its historic eruptions call for detailed geophysical monitoring, which is, however, complicated by the remoteness of the site and its harsh environmental conditions. Within this context, emerging Distributed Acoustic Sensing (DAS) technologies offer new opportunities for densely spaced deformation measurements in a broad frequency band from millihertz to kilohertz (Hartog, 2017; Lindsey et al., 2020; Paitz et al., 2021; Lindsey and Martin, 2021). In particular, the relative ease of deploying a fibre‐optic cable on a glacier yields a measurement density that would be difficult to achieve with conventional seismic instruments (Walter et al., 2020; Klaasen et al., 2021).

In spring 2021, we trenched 12.5 km of a Solifos BRUfield fiber‐optic cable around and within the caldera of Grímsvötn using a custom‐made sled, towed by a Snowcat (Klaasen et al., 2022). In a research hut, located at Grímsfjall, the highest point of the caldera rim, we connected a single‐mode fiber to a Silixa iDAS v2 interrogator. It collected strain rate measurements with 8 m channel spacing, 10 m gauge length, and 200 Hz sampling frequency for 3 weeks, starting on 6 May 2021. A summary of the experimental setup is shown in Figure 1. Although the thickness of the ice cap in the area is roughly constant around 300 m, the lake level varies over time scales of several years. Prior to flooding events as in December 2021, it may reach 100–150 m (Bethier et al., 2006).

Despite its geologic complexity, we will show that many aspects of Grímsvötn’s low‐frequency seismic response can be understood in terms of simple effective physical systems.

Figure 2 shows a typical record section that focuses on the last 4.5 km of the cable, largely located inside the caldera. Deformation with a dominant frequency above 10 Hz, induced by a local earthquake and reaching peak strain rates above 100 nanostrain/s, is clearly visible in Figure 2a, which shows a 40 s long record section starting at UTC01:00:30 on 7 May 2021. During the 3 weeks of the experiment, nearly 2000 similar local events could be detected (Thrastarson et al., 2021). The earthquake signal is superimposed onto a nearly harmonic oscillation at lower frequency and maximum strain rates on the order of 10 nanostrain/s. As shown in Figure 2b, these oscillations persist during the morning hours of 7 May 2021, with slight variations on a one‐minute time scale but overall constant peak amplitudes around 10 nanostrain/s.

The densely sampled DAS data reveal that the low‐frequency oscillations are not a localized phenomenon but affect the floating part of the ice sheet as a whole. As indicated, for instance, by the dashed vertical line in Figure 2a, the parts of the caldera between ∼9.0 and ∼11.5 km along the cable periodically extend, while adjacent parts compress, and vice versa. This suggests that we observe flexural gravity waves, i.e., oscillatory bending of the ice sheet coupled to gravity waves within the water layer, but also to the underlying solid Earth. Though being subject to different boundary conditions, these flexural gravity waves are similar to those in Antarctic ice shelves (Bromirski et al., 2017; Lipovsky, 2018; Olinger et al., 2022). With the limited coverage of the caldera, the complete mode shape can unfortunately not be estimated. It therefore remains unclear if we observe the fundamental mode or an overtone.

Spectral analyses of the DAS channels within the caldera, shown in Figure 3, quantify the bandwidth of the oscillations at half peak amplitude to Δf0.02  Hz around the dominant frequency f0=0.22  Hz. The width of the observed amplitude spectrum can largely be explained by that of a 1D damped harmonic oscillator,
in which A denotes amplitude in strain/s2, ω circular frequency, and ω0=2πf0. The quality factor Q controls the width of the peak and has a value of around 20. This suggests that the system as a whole, despite being 3D and geometrically complex, effectively behaves as a forced resonator that dissipates around 2πQ131% of its total energy per oscillation cycle. Figure 3a also confirms that the spectral amplitudes of the ice sheet oscillations are approximately constant in the morning hours of 7 May 2021.

These spectral properties are in contrast to the characteristics of the ocean wave‐generated secondary ambient field (Peterson, 1993; Gualtieri et al., 2013; Nakata et al., 2019), shown for the coastal station KVI in Figure 3b, in which this contribution to the seismic wavefield is expected to be the largest. Compared to the ice sheet oscillations, ambient vibrations have a similar central frequency of ∼0.28 Hz but with significantly larger bandwidth of Δf0.07  Hz and spectral amplitude variations of nearly 300% within only 12 hr. As shown in Figure S1, available in the supplemental material to this article, the spectral characteristics at station KVI are very similar to those at other stations in southern Iceland.

Figure 3c extends the spectral analysis of the DAS channels and station KVI to the whole duration of the experiment. The time‐dependent spectra display visible correlations, suggesting that ambient vibrations are a significant driver of the ice sheet oscillations. However, there are several time periods when strong ice sheet oscillations could be observed, despite low ambient vibration amplitudes. These include, as seen before, the morning hours of 7 May 2021, but also others, extending from few hours to several days.

Though station KVI is only one among many in southern Iceland, it is still remarkably representative, as illustrated in Figures S1 and S2. The spectral properties of ambient vibrations recorded at broadband seismometers in southern Iceland are nearly identical. This leads to the important preliminary conclusions that the ambient field in the region is spectrally largely homogeneous, and that the anomalous oscillations in the Grímsvötn caldera require some residual driving forces.

Despite the complexity of the volcanic system, the nature of the data permits plausible simplifications to estimate the effective residual forces. At frequencies below ∼0.3 Hz and seismic phase velocities in the range of several kilometers per second for volcanic rocks (Lowrie and Fichtner, 2020), the wavelengths of the ambient and residual deformations that drive ice sheet oscillation are on the order of 10 km. Consequently, the observations are largely insensitive to spatial variations of the forces at scales below the dimension of the Grímsvötn caldera. This allows us to represent them as space‐independent quantities, as illustrated schematically in Figure 4a. Furthermore, the low deformation amplitude on the order of 10 nanostrain/s implies that the system is effectively linear. Hence, denoting by Fk(ω) the total effective force in time interval k, we have
in which uk(ω) is the observed ambient velocity spectrum, and rk(ω) the yet unknown residual forcing. The proportionality factor F, with unit s/m, states that the effective force Fk that drives the system should be linearly related to ground acceleration ω[ukobs(ω)+rk(ω)]. Owing to the previously noted spatial homogeneity of the ambient field, we may use any vertical‐component seismometer spectrum in southern Iceland for ukobs(ω). In the following, we use the spectrum of station BJK (Fig. 1b), which is the closest to the Grímsvötn caldera. Taking the effective system response g(ω) from equation (1), we may model the oscillation spectrum as

For given F and Q, the residual forcing can be found through the solution of a linear least‐squares problem (Tarantola, 2005; Fichtner, 2021). Combining this with a grid search, yields values for F and Q that minimize the norm of rk(ω), thereby providing the minimum residual forcing as a function of frequency and time interval, displayed in Figure 4b. The minimum residual forcing produces a root mean square misfit between observed spectra skobs(ω) and modeled spectra sk(ω) of 1.65, for an observational error standard deviation of 20×103 nanostrain, estimated conservatively from hourly amplitude spectra of the DAS channels. Furthermore, the amplitudes of the minimum residual forcing are significantly above the noise level of the seismometer recordings and in fact comparable to the ambient seismic noise amplitudes. This indicates that our simplistic model explains the observations to within their errors, and that no additional model complexities are required. A more detailed, time‐ and frequency‐dependent analysis of the remaining difference between skobs(ω) and sk(ω) can be found in Figure S3.

As shown in Figure 4b, the inferred minimum residual forces acted continuously during the experiment and at least within the frequency band covered by the caldera resonance peak. Temporal variations occurred over several days, with a pronounced activity peak around day 17 (24 May 2021).

Local seismicity is unlikely to be a major contributor to the observed resonance for two reasons. First, local event signals have frequencies significantly above 1 Hz and therefore cannot excite resonance at 0.22 Hz. Second, even around 0.22 Hz, the damping of the resonator attenuates transient signals to below the noise level within less than a minute. This is significantly shorter than the average interevent duration of ∼15 min (Thrastarson et al., 2021). A similar argument excludes glacial flow, amounting to ∼5 cm/day as a source, because it acts in a far lower frequency band.

In the context of earlier observations of low‐frequency continuous seismic signals at other volcanoes (Kawakatsu et al., 2000; De Martino et al., 2005; Wassermann, 2012), the residual forces can plausibly be explained in terms of volcanic tremor. In the absence of visible magmatism in May 2021, the tremor is likely to be caused by geothermal activity, which has clear surface expressions in the form of fumaroles and cauldrons near the eastern flank of the caldera rim.

It would clearly be desirable to model the ice shelf oscillation in greater detail. Although this would technically be possible using spectral‐element wavefield simulations that permit complex geometries (Komatitsch and Vilotte, 1998; Afanasiev et al., 2019), the required pieces of information are not available. This includes, the most importantly, detailed maps of ice thickness and lake depth, which can vary substantially over weeks, that is, more rapidly than any 2D geophysical survey could be repeated.

In the same context, it is important to note that more complex physical models with a larger number of free parameters than the 1D damped harmonic oscillator are not required by the available data. This is because the observations can already be explained to within the errors by a single residual forcing term. Therefore, Bayes’ theorem would ensure that models with more parameters are automatically less plausible. As a consequence, the residual forcing must be understood as an effective force that encapsulates individual forcing mechanisms that our data cannot resolve.

Regardless of the interpretation, it is evident that the floating ice sheet acts as an amplifier of signals that may otherwise be overwhelmed by ambient and instrumental noise. This suggests that similar natural resonators might be used in a targeted fashion for the monitoring of subglacial volcanoes, including, for example, Katla and the Skaftár cauldrons. Yet, while amplifying certain information, the resonance within a narrow frequency band also eliminates other frequencies that may be needed to infer more detailed source properties (Chouet, 1996; Lipovsky and Dunham, 2015).

Finally, our experiment demonstrates the logistical feasibility of installing a large fiber‐optic sensing network in a subarctic environment in which the deployment of a similar array of conventional seismic instruments would not be possible. Trenching the cable ∼30 cm into the ice provides a data quality that permits the detection of deformation as small as few nanostrain/s at frequencies from ∼0.1 to 20 Hz. This study thereby paves the way toward dense seismic monitoring of remote, glacier‐covered volcanoes under harsh environmental conditions.

The Grímsvötn Distributed Acoustic Sensing (DAS) dataset has a volume of around 1 TB, making it too large for remote access and download. However, the dataset will be made available on hard drives upon request. A video showing some fieldwork impressions can be found at https://www.youtube.com/watch?v=J7cxVZvgyWQ&t=0s (last accessed June 2022).

The authors acknowledge that there are no conflicts of interest recorded.

The authors gratefully acknowledge phenomenal support from Silixa throughout all stages of this experiment. This work was partially funded by the Real‐time Earthquake Risk Reduction for a Resilient Europe project (RISE) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement Number 821115). Constructive comments by Brad Lipovsky and an anonymous reviewer helped us to improve the article.