Fiber‐Optic Airplane Seismology on the Northeast Greenland Ice Stream

Firn is the transition material between snow and ice formed by a sequence of processes, including the settling and rounding of crystals, recrystallization, and viscoplastic deformation (Gow, 1968). Although firn often accounts for only a small fraction of ice sheet mass and volume, detailed knowledge of its properties is essential for ice sheet mass balance calculations (e.g., Helsen et al., 2008), ice core climatology (e.g., Parrenin et al., 2012), estimates of surface melt (e.g., van den Broeke, 2005), and the correction of seismic reflection data for near‐surface effects (e.g., Schlegel et al., 2019).

For several decades, firn properties have been estimated noninvasively using seismic refraction data (Thiel and Ostenso, 1961; Kohnen and Bentley, 1973; Schlegel et al., 2019; Hollmann et al., 2021) and horizontal‐to‐vertical spectral ratios (Lévêque et al., 2010), sometimes combined with empirical conversions from seismic wavespeeds to density (Kohnen, 1972; Diez et al., 2014). Alternatively, the distribution of S wavespeed may be estimated from seismic noise interferometry (Zhou et al., 2023).

New opportunities for seismic studies of firn come from multiple directions. Distributed acoustic sensing (DAS) combines unprecedented spatial coverage with a large bandwidth (Lindsey et al., 2020; Lindsey and Martin, 2021; Paitz et al., 2021), and has proven to be logistically feasible and beneficial in glacial environments (Walter et al., 2020; Hudson et al., 2021; Klaasen et al., 2021, 2022; Fichtner et al., 2022). Also recently, unconventional seismic sources related to human activity have been recognized more widely as valuable carriers of information on subsurface structure (Dou et al., 2017; Sager et al., 2022).

In this work, we harness DAS data from an airplane landing near the EastGRIP ice core drilling site located in the upstream part of the Northeast Greenland Ice Stream (NEGIS; Fig. 1a). With a length of nearly 600 km, NEGIS accounts for ∼12% of Greenland’s total ice discharge (Rignot and Mouginot, 2012), thereby attaching special relevance to detailed firn models. The strength and duration of the airplane source produced seismic data of remarkable clarity, containing 15 visually discernable wave propagation modes, many of which can be identified and used to constrain S wavespeed with meter‐scale resolution and uncertainties of few meters per second depending on depth.

Realizing that an airplane landing may be the strongest available seismic source, we deployed 3000 m of Solifos BRUfield fiber‐optic cable (light‐weight tactical, nonmetallic loose‐tube) for DAS recording. To maximize sensitivity to waves excited near the expected landing spot, the cable runs perpendicular to the skiway, except for the first few hundred meters as shown in Figure 1. One of four single‐mode fibers was connected to a Silixa iDAS with 10 m gauge length, and interrogated with 1 kHz sampling frequency and 2 m channel spacing. To ensure good coupling, a snowcat produced a trench for the cable, which was then covered by snow. We estimate the average burial depth of the cable at >50 cm. For our analyses, we use the last 2500 m of the cable, which are nearly linear, orthogonal to the skiway and at sufficient distance from the EastGRIP camp to avoid the presence of significant anthropogenic noise, mostly originating from the generator.

The C‐130 Hercules touched down on a groomed skiway at UTC 17:02:06 on 26 July 2022 at a position that makes an angle between ∼3° and 7° to the linear cable segment. Its weight and sinking speed are estimated around 60 t and 1 m/s, respectively. A record section of the raw strain rate recordings is displayed in Figure 2a. Although individual components of the wavefield are difficult to discern in the time–space domain, a transformation to the frequency–phase velocity domain via a simple 2D Fourier transform shown in Figure 2b reveals ∼15 visually distinguishable wave propagation modes. These include several Rayleigh wave modes with cutoff velocity around 1900 m/s, a leaky mode propagating faster than ∼1900 m/s, as well as pseudoacoustic modes that correspond to P waves trapped in the firn layer, which acts as a waveguide. Thanks to the exceptional clarity of the data, no additional preprocessing of the raw data was needed.

The modes were identified through direct synthesis of the frequency–phase speed spectrum for a surface source (Kennett, 2023) that we employed for the construction of an initial model, described in more detail subsequently. This approach allows mapping of modal extension into the leaking mode regime and so includes the P mode behavior.

The unusual bandwidth combined with the large number of identifiable modes carries detailed information about the seismic properties of the upper few hundred meters of the ice sheet and the firn layer in particular. Radar soundings constrains the topographic variations of these layers to be on the order of only few meters over several kilometers distance (e.g., Franke et al., 2022), which enables the use of modeling and inversion techniques for stratified media without significant approximation errors. To account for the nonzero incidence angle of ∼5° on average, we multiply all apparent phase velocities with cos(5°) ≈ 0.996.

To avoid biases introduced by subjective discretization and regularization, we implement a Backus–Gilbert inversion of the dispersion data (Backus and Gilbert, 1968, 1970; Fichtner, 2021). This comes with the additional benefit of automatically providing objective model errors and resolution (averaging) lengths. A prerequisite for Backus–Gilbert inversion is the availability of an initial model that explains the main features of the data and thereby permits a reasonable linearization. In the following paragraphs, we first describe the initial model construction before presenting the results of the Backus–Gilbert inversion.

The initial model, displayed in Figure 3a, is the result of a trial‐and‐error procedure primarily aiming to explain the fundamental‐mode Rayleigh wave, and the P‐wave modes $P0$ and $P1$ to within the observational standard deviations, estimated as the half‐width of the phase velocity–frequency peaks in Figure 2b. For this, we took advantage of the long‐established prior knowledge that the depth dependence of seismic velocities in firn approximately follows an exponential (Brockamp and Pistor, 1967; Kohnen and Bentley, 1973). Using constant instead of exponential velocities just within the top 1 m significantly improves the fit to the fundamental Rayleigh mode. This feature plausibly reflects the more rapid compaction of snow around the EastGRIP camp in response to human activity. The initial density model is scaled to the initial $vp$ model following the relation proposed by Kohnen (1972).

Figure 3b illustrates that the initial model, by design, reproduces the fundamental Rayleigh mode and the P‐wave mode dispersion measurements. This match enables the identification of the other modes, including the realization that $R1$ is only weakly excited. Discrepancies between observed and modeled phase velocities remain for the higher Rayleigh modes, especially for the fourth and fifth overtones. These are the sources of information about more structural details in the $vs$ model that we consider in the section on Backus–Gilbert inversion.

The initial model serves in the calculation of sensitivity kernels $Gi$⁠, which play a double role in Backus–Gilbert theory. They link the difference $δci=ciobs−ci$ between observed and calculated phase velocities for the ith measurement to the $vs$ distribution via $δci=∫Gi(z′)δvs(z′)dz′$⁠. In addition, they are basis functions of averaging kernels, $A(z′,z)=∑iai(z)Gi(z′)$ that produce a depth‐averaged version of the S wavespeed, $v‾s(z)=∫A(z′,z)vs(z′)dz′$⁠. Though $vs(z′)$ is not directly observable, the average can be computed directly and uniquely from the data as $v‾s(z)=∑iai(z)δci$⁠. Basis coefficients $ai(z)$ that produce averaging kernels with optimal concentration around a target depth z can be found by maximizing the deltaness $λ(z)=12∫(z−z′)2A2(z′)dz′$⁠, which simultaneously provides a conservative estimate of the averaging or resolution length. Standard deviations of the averages, $σ(z)$⁠, can be computed from the standard deviations of the observations, $σobs$ via $σ2(z)=∑i[ai(z)σiobs]2$⁠.

For the inversion we employ all available pairs of observed and calculated phase velocities with 1 Hz frequency spacing shown in Figure 3b. Omitting the first overtone $R1$⁠, because of its weak excitation and large observational errors results in 94 data points. In the construction of S‐wavespeed models, the averaging length $λ$ and the standard deviation $σ$ are antagonistic. Increasing $λ$ decreases $σ$ and vice versa, thereby giving rise to infinitely many models with a different balance between accuracy and resolution.

Two of these models are displayed in Figure 4. In accord with the mostly positive phase velocity residuals $δci$⁠, the less conservative one in Figure 4a shows several high‐velocity deviations from the initial model, especially around 3 and 20 m depth. These are less pronounced in the more conservative version in Figure 4b, which has a larger averaging length but lower standard deviations. Within the upper 10 m, both the models have a vertical resolution or averaging length of few meters, increasing to few tens of meters around 100 m depth. Standard deviations of the inferred $v‾s$ reach the minimum well below 10 m/s near the surface and approach 100 m/s at several hundred meters depth. Code to produce a suite of alternative models is available in the model website referenced in the Data and Resources section.

The averaging length scales of ∼1 to few tens of meters within the upper 100 m can be considered conservative, because $λ$ is a mathematically conservative measure that gives large weight to distant side lobes in the averaging kernels. In fact, ice core analyses show that firn and ice properties below several meters depth vary smoothly over seismic wavelength scales (Vallelonga et al., 2014) and can be approximated by linear functions within the width of the averaging kernels. Though it is difficult to turn this into a quantitative argument, the smoothness suggests that the difference in depth in which $v‾s(z)$ and $vs(z)$ reach a specific value is well below the averaging length.

Although dispersion curves and seismic travel times are primarily sensitive to wavespeeds, most glaciologically relevant information derives from density profiles. Empirical scaling relations have been estimated primarily from the comparison of seismic refraction profiles and direct density measurements on drill cores. Despite inherent limitations detailed in the Discussion section, the resulting density profiles carry useful approximate information, for instance, about the depth of the firn–ice transition at a density of $830kg/m3$⁠.

We employ the empirical scaling derived from SH‐wave diffraction by Diez et al. (2014). To make their relation applicable to our SV‐wave model, we multiply SH wavespeeds by a factor of 0.95 that accounts for the typical 5% radial anisotropy in firn layers (e.g., Schlegel et al., 2019). In the modified scaling, visualized in Figure 5a, a density of $830kg/m3$ corresponds to an S wavespeed of 1685 m/s. Using the model and uncertainties from Figure 4a, the averaged S wavespeed $v‾s$ exceeds $1685m/s±σ$ at depths between 65 and 71 m (Fig. 5b). This interval differs only at submeter scale when we choose other models with a different balance between accuracy and resolution, suggesting that the result is robust in this sense. Direct density measurements on firn cores at EastGRIP place the firn‐density transition at ∼65 m (Vallelonga et al., 2014), in good agreement with our inference.

The scaling from wavespeed to density also translates the pronounced high‐velocity anomaly near 20 m depth into a high‐density anomaly. Although such an anomaly is not visible in the currently available firn core density measurements, the depth range around 20 m is characterized by anomalously thin annual layers (Vallelonga et al., 2014). Because layer thickness affects the strength of induced radial anisotropy, this suggests that the seismically inferred high‐density anomaly may be an artifact introduced by a scaling relation that is unaware of seismic anisotropy.

Much of this work rests on the exceptional data quality, which has two major contributors in addition to low noise and the sheer strength of the source:

First, the strong velocity increase by a factor of ∼10 over less than 100 m depth leads to wave trapping that is particularly effective. For comparison, P and S wavespeeds in the Earth’s crust differ from those in the mantle only by a factor of ∼2 (e.g., Kennett et al., 1995). As a consequence, the character of the spectrum in Figure 2 resembles the displacement spectrum of the Sun (Kosovichev et al., 1997), which has a similar wavespeed contrast as a firn layer (Christensen‐Dalsgaard et al., 1996). Furthermore, the lateral invariance of the medium produces clearly defined modes that hardly vary along the cable.

Second, the availability of densely spaced DAS data is critical. With 2 m channel spacing, we achieve unaliased sampling of the wavefield, which has wavelength as low as ∼15 m for the fundamental Rayleigh mode at 30 Hz. Because the cable is ∼10–100 times longer than a wavelength, we obtain the resolution in wavenumber space that allows us to observe a large number of individual modes and measure their phase velocities with errors as low as a few meters per second. A direct consequence are the short averaging lengths and low uncertainties of the inferred S wavespeed profile.

The success of the Backus–Gilbert inversion rests on the availability of a useful initial model, and in particular on the realization that the second strongly excited Rayleigh mode is $R2$ and not $R1$⁠. The latter derives from the prior knowledge that firn velocity profiles are approximately exponential. Ignoring this prior knowledge and trying to fit observed $R2$ to calculated $R1$ phase velocities with a global simulated annealing algorithm produced entirely implausible velocity models with a fast sequence of low‐velocity layers and extreme P‐wave anisotropy. That modes of lower order are less strongly excited than modes of a higher order is unusual but not impossible. It has previously been observed, for instance, in low‐velocity sedimentary basins (Cruz‐Atienza et al., 2016), and it can easily be explained mathematically by an approximate orthogonality of a modal eigenfunction and the forcing term (Gilbert, 1971).

We tested the effect of linearization by rerunning the inversion with some of the inversion results as new initial model. However, such a one‐step iterative improvement produced results that differed from the ones shown in Figure 4 by less than 1%, thereby suggesting that the problem is insignificantly nonlinear.

The empirical scaling from wavespeed to density inherits fundamental limitations from vastly different resolution lengths. Most scaling relations (e.g., Kohnen, 1972; Diez et al., 2014) link in situ or drill core measurements of density at centimeter scale to seismically inferred velocities with resolution lengths of many meters. As a consequence, each scaling depends on the seismic data used to derive it, and it may therefore not easily be transferred to other seismic experiments with different averaging kernels. This inconsistency is exacerbated by the presence of various types of anisotropy, which is neither represented in the scaling relations nor sufficiently resolvable with most seismic data. In summary, it does not seem unlikely that errors in the scaling relations dominate errors related to the insufficiency of the seismic data. Hence, small‐scale details in the inferred density profiles, such as the high‐density anomalies in Figure 5b, should be interpreted with caution.

We presented DAS data produced by an airplane landing near the EastGRIP ice core drilling site on the Northeast Greenland Ice Stream. The data are of exceptional clarity, especially in the frequency–phase velocity domain, in which around 15 wave propagation modes are easily visible. An inversion of seven identifiable modes produces an S‐wavespeed profiles with resolution lengths as low as a few meters and standard deviations in the range of only 10 m/s. Empirically scaling to density provides an estimate of the firn–ice transition, which agrees with direct firn core measurements. Our experiment as a whole illustrates the potential of DAS combined with unconventional sources to provide detailed inferences of firn structure.

Yet, the work presented here is only a step towards a more comprehensive inversion. Possible improvements include (1) the incorporation of the leaky‐mode observations, using, for instance, adjoint‐based sensitivity kernels from numerical wave propagation, (2) the extension of the Backus–Gilbert method to multiple parameter classes, such as P wavespeed and density, which may make significant contributions at shallow depth, and (3) a fully nonlinear inversion for all relevant parameter classes using Monte Carlo sampling.

All data as well as Python Jupyter notebooks to reproduce Rayleigh‐wave dispersion and inversion are available on the Earth model website of the ETH Seismology and Wave Physics Group (https://www.swp.ethz.ch → Models). Code for the synthesis of the frequency‐phase speed spectrum is available at Zenodo (doi: 10.5281/zenodo.7384357). A video with field work impressions can be found at https://www.youtube.com/watch?vCF7nZ7mxzSo. All websites were last accessed in May 2023.

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

Andreas Fichtner and Coen Hofstede gratefully acknowledge support by the whole EastGRIP team, including the provision of all necessary infrastructure, help with trenching the cable, as well as logistic and technical advice. Invaluable technical support before and during the experiment was provided by Silixa (Athena Chalari and support team) and Solifos (Andrea Fasciati). Olaf Eisen and Dimitri Zigone were supported by the CHIPSM grant of the University of Strasbourg Institute for Advanced Studies. EastGRIP is directed and organized by the Centre for Ice and Climate at the Niels Bohr Institute, University of Copenhagen. It is supported by funding agencies and institutions in Denmark (A. P. Moller Foundation, University of Copenhagen), the United States (U.S. National Science Foundation, Office of Polar Programs), Germany (Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research), Japan (National Institute of Polar Research and Arctic Challenge for Sustainability), Norway (University of Bergen and Trond Mohn Foundation), Switzerland (Swiss National Science Foundation), France (French Polar Institute Paul‐Emile Victor, Institute for Geosciences and Environmental Research), Canada (University of Manitoba), and China (Chinese Academy of Sciences and Beijing Normal University).