Distributed Acoustic Sensing of Seismic Properties in a Borehole Drilled on a Fast‐Flowing Greenlandic Outlet Glacier

Distributed Acoustic Sensing (DAS) is a new technology in which seismic energy is detected, at high spatial and temporal resolution, using the propagation of laser pulses in a fiber‐optic cable. We show analyses from the first glaciological borehole DAS deployment to measure the englacial and subglacial seismic properties of Store Glacier, a fast‐flowing outlet of the Greenland Ice Sheet. We record compressional and shear waves in 1,043 m‐deep vertical seismic profiles, sampled at 10 m vertical resolution, and detect a transition from isotropic to anisotropic ice at 84% of ice thickness, consistent with the Holocene‐Wisconsin transition. We identify subglacial reflections originating from the base of a 20 m‐thick layer of consolidated sediment and, from attenuation measurements, interpret temperate ice in the lowermost 100 m of the glacier. Our findings highlight the promising potential of DAS technology to constrain the seismic properties of glaciers and ice sheets.

Glaciological seismic reflection surveys typically install receivers at, or close to, the glacier surface. Although this is logistically practical, the observed seismic response must be corrcted for propagation between the source, the target horizon, and receviers. This is typically done using hyperbolic traveltime approximations and, in the absence of ground-truth constraint, impacts the accuracy of physical property estimates. The lack of ground-truth data can be mitigated using attributes of the seismic data set itself, including calibrating reflection coefficients using the amplitude ratio of primary and multiple arrivals (Booth et al., 2012;King et al., 2008) and measuring source strength using direct-wave amplitude decay (Holland & Anandakrishnan, 2009;Peters et al., 2008Peters et al., , 2012, but seismic quantities often remain underconstrained. The problem is exacerbated for passive seismic data sets, where source positions are unknown and location algorithms must assume background velocity/attenuation structures (e.g., Podolskiy & Walter, 2016;Röösli et al., 2014).
Borehole seismic deployments are advantageous because they directly sample in situ seismic properties, thus enabling surface seismic responses to be calibrated. Interpretations of "Vertical Seismic Profile" (VSP) data invoke fewer traveltime assumptions than for surface surveys; hence, properties can be more accurately defined . These properties can extend to velocity and velocity anisotropy (Diez et al., 2015), attenuation (Beckwith et al., 2017), and reflection coefficients (Lira et al., 2012). Furthermore, VSP surveys can recover englacial seismic properties even where there is no internal reflectivity .
Recently, Distributed Acoustic Sensing (DAS) has been applied in a variety of seismic settings, notably where long-term variations in seismicity are monitored (e.g., Ajo-Franklin et al., 2019;Miller et al., 2016;Sladen et al., 2019). DAS is a novel technology for recording seismic responses, involving the transmission of laser pulses through a fiber-optic cable from an interrogator unit. As it travels, a fraction of the laser pulse is backscattered via Rayleigh scattering, and the phase lag between backscattered pulses changes as the cable is strained (e.g., by the passage of a seismic wavelet). The interrogator monitors the changing phase lag and thereby reconstructs the real-time seismic response along the cable length (Hartog, 2017). The fiber-optic cable therefore acts as a near-continuous geophone string; although the backscattered pulse is sampled at high spatial resolution, the seismic wave is reconstructed across a larger interval termed the "gauge length," which represents the effective spatial resolution of the cable. Gauge lengths are typically set to be greater than~½ of the shortest expected seismic wavelength, to ensure faithful wavelet reconstruction and boost signal-to-noise ratios (Dean et al., 2017). Although finer resolution could be achieved with a densely sampled geophone array, this would require much instrumentation and be cumbersome to deploy in a borehole. Geophones would also be vulnerable to loss in a deforming/freezing borehole, whereas the fiber-optic cable is a long-term and/or disposable installation.
A further limitation of DAS, compared to a three-component geophone, is that sensitivity to particle motion varies as a function of the incidence angle, θ, of seismic energy to the cable axis. For P and vertically polarized shear (SV) waves, this sensitivity varies as cos 2 θ and sin2θ, respectively (Mateeva et al., 2014), and there is no sensitivity to horizontally polarized shear (SH) waves. Nonetheless, DAS VSP data sets offer much potential for monitoring and characterizing the in situ seismic properties of glacier ice.
Here, we interpret data from the first glaciological borehole deployment of DAS technology, performed during the summer of 2019 on Store Glacier in West Greenland. We highlight the scope of the data set, determining (i) vertical P wave velocity and attenuation structure of the ice column, (ii) S wave velocities and the implied Poisson's ratio, and (iii) thickness of a subglacial sediment layer.

Field Setting and Deployment
Store Glacier (Greenlandic: Sermeq Kujalleq) is a marine-terminating outlet glacier that drains~34,000 km 2 of the Greenland Ice Sheet (Rignot et al., 2008). Previously, the subglacial environment of the glacier was accessed at Site S30 (inset, Figure 1) with boreholes drilled by pressurized hot water. Here, flow speeds of 600 m a −1 occur in response to high basal water pressure ) over a sedimentary bed .
In July 2019, a suite of vertical boreholes was drilled in the center of supraglacial lake Site L028 (70.56793°N, 50.08697°W; Figure 1), which drained via hydrofracture on 31 May 2019 (Chudley et al., 2019), 28 km upstream of the Store Glacier terminus. Borehole BH19c was drilled to the glacier bed, at~1,043 m depth, and instrumented with a Solifos BRUsens fiber-optic cable. This borehole, and two others located 70 m away, drained rapidly on breakthrough implying intersection with either the basal interface or a basal crevasse. With similar drainage observed previously at S30 , we exclude basal crevasses: it would be unlikely for each borehole to intersect these since they form less abundantly than surface crevasses. We therefore favor the simpler explanation that BH19c drained by instersecting the basal interface. The cable contained two single-mode fibers (for DAS) and four multimode fibers (for distributed temperature sensing, DTS, to be reported elsewhere), enclosed in a gel-filled stainless-steel capillary tube, with 4.8 mm-diameter polyamide sheathing and steel wire reinforcement. DAS VSPs were then recorded during a 3-day period, between 6 and 8 July 2019.
A Silixa iDAS™ interrogator (Figure 1b), recording at 4,000 Hz sampling frequency, established the seismic response every 1 m along the cable using a 10 m gauge length (~½ of the expected seismic wavelength). The vertical cable responds sensitively to vertically propagating P waves but is insensitive to those traveling horizontally; the converse holds for SV waves.
Seismic energy was generated with a 7 kg sledgehammer striking a polyethylene impact plate (dimensions 0.20 × 0.20 × 0.05 m). Shots were made (i)~1 m away from the borehole top (effectively zero offset; Figure 1c); (ii) along azimuths parallel, orthogonal, and at 45°to ice flow (~217°), up to 500 m offset from the borehole; and (iii) radially, at 30°intervals, at 300 m nominal offset from the borehole. For stacking, 30 hammer strikes were recorded at each location over a~4 min period. Recorded seismic wavelets were consistent: the typical correlation coefficient between successive traces in a shooting sequence is >0.96.
The fiber-optic cable was monitored in a passive mode. Zero-offset shots were extracted from the record by identifying the earliest arrival of seismic energy at the approximate time of shooting. For offset shooting, impact times were recorded with a GPS-synchronized surface geophone installed within 1 m of the impact place and used to extract data from the DAS record. These times are accurate to ±1 ms; residual static corrections, typically less than ±0.75 ms, were applied to synchronize shots prior to stacking.   (Figure 1a), hereafter termed "300 m offset." Prominent arrivals and initial velocity assessments are annotated for each record. The enlarged panels in Figure 2b show direct-wave first-break picks for each trace (red points) and their mean and standard deviation within each 10 m gauge length (blue symbols). VSP acquisition geometries are described in the supporting information.

Zero-Offset VSP
The direct P wave is prominent in the zero-offset record (Figure 2b), having linear moveout and no abrupt gradient changes, hence indicating no strong englacial P wave velocity (v P ) contrasts. However, linear trends fitted in shallow and deep intervals suggest a progressive v P increase with depth:~3,750 m s −1 in the upper 60 m, versus~4,000 m s −1 in the lowermost 100 m (insets i and ii). The direct-wave amplitude spectrum, stable throughout the VSP, has dominant frequency and bandwidth of~200 Hz (inset, Figure 2a). The absence of a basal reflection at 1,043 m depth is consistent with low-reflectivity bed conditions previously observed by Hofstede et al. (2018).
The low-amplitude events that are parallel to the direct P wave, but with hyperbolic moveout between 50-70 ms arrival time, are interpreted as P wave arrivals scattered from surface crevasses. The linear arrival following the direct wave has a velocity of~1,825 m s −1 , close to the typical S wave velocity (v S ) in ice. This is interpreted as a Stoneley wave (Cheng & Toksöz, 1981), a vertically traveling phase that propagates at a solid-solid interface, assumed to be the refrozen borehole wall.

The 300 m Offset VSP
The offset VSP features a prominent hyperbolic direct P wave (Figure 2c), detectable only below~70 m given the insensitivity of the cable to such near-perpendicular arrivals. The best fit hyperbola to the direct P wave traveltimes gives v P of 3,810 ± 12 m s −1 (Figure 2d). The later hyperbolic event (zero-depth arrival timẽ 160 ms) is the direct SV wave: its best fit hyperbola defines v S of 1,870 ± 2 m s −1 . SV wave amplitudes decay rapidly with depth due to combined propagation losses and cable sensitivity, but amplitude variations may also be complicated by the anisotropic radiation pattern of the shear component from the impact source (e.g., Hardage, 2011).
Consistent with the zero-offset record, neither arrival indicates any abrupt velocity changes. Backscattered energy is less prominent than in the zero-offset record although a strong scattered arrival is present, potentially reflected from a fracture~110 m south of the shot location ( Figure 1a). There is no apparent P wave basal reflection, but a low-amplitude reflection is visible in the lowermost 150 m (panel iii). This event does not converge on the direct wave but lags it by~44 ms at the deepest position (~1,000 m) that it is clearly visible. It is therefore interpreted as a subglacial reflection, propagating through a velocity:thickness model that defines a~44 ms two-way traveltime.

Vertical v p Model
Owing to the vertical cable and wave propagation, v P in any interval is the depth difference between two observation points, divided by the traveltime of energy between them. We evaluate traveltimes by cross-correlating traces from different depth intervals and recording the time lag of the cross-correlation peak. Thus, the velocity obtained is the average P wave interval velocity within the correlation gap.
Vertical v P trends are estimated for correlation gaps of 20, 50, and 100 m (gray, red, and blue, respectively, in Figure 3a). Curve widths represent the velocity uncertainty from misidentifying the peak cross-correlation lag by ±0.25 ms (one digital sample) and correspond to ±180, ±75, and ±40 m s −1 for the 20, 50, and 100 m gaps, respectively. The 20 m gap offers the highest vertical resolution but is noisy since uncertainties are a greater proportion of the interval traveltime. The longer gaps give smoother trends and suggest a constant v P of~3,700 m s −1 in the upper 800 m of the glacier; v P then rises steadily and exceeds 4,000 m s −1 between 880 and 950 m depth. Below 950 m, velocities progressively reduce to a minimum observed value of 3,730 m s −1 . These trends are evident regardless of the correlation gap, hence are considered to be genuine.

Vertical Attenuation Model
The amplitude, A x , of a propagating seismic wave decreases due to geometric spreading and anelastic attenuation losses, as where A 0 is initial amplitude, G(x) is geometric spreading losses, α is an attenuation coefficient, and x is the distance traveled. For the zero-offset VSP, x is the depth, z, of the observation point along the fiber-optic cable, hence G(x) = z −1 . The depth-averaged attenuation coefficient, α z , is therefore defined by rearranging Equation 1 as and the gradient in a plot of ln zA z A 0 against z gives −α −1 z .
The P wave attenuation coefficient is related, via v P , to the seismic quality factor, Q, by where f is wavelet frequency. Q has been related to englacial temperature (Peters et al., 2012), with warmer ice associated with lower Q values and therefore increased attenuation. Peters et al. (2012) measure Q from englacial reflections using spectral ratios (Dasgupta & Clark, 1998;Gusmeroli et al., 2010). In our data, the consistency of spectra (Figure 2a) suggests that Q is high and spectral ratios could be influenced more by noise than genuine attenuation signatures. We therefore measure α z and Q from the amplitude decay of the direct P wave.
Following corrections for spreading losses and assuming A 0 = 1, the variation of direct-wave amplitudes ( Figure 3b) has two linear sections. The nonphysical increase of amplitude to~50 m depth is attributed to interference with the Stoneley wave. Between 50 and 900 m, amplitude decay implies α z of [5.5 ± 0.1] × 10 −4 m −1 ; in the lowermost~100 m, α z is [118 ± 5] × 10 −4 m −1 . Assuming the v P derived previously and a dominant wavelet frequency of 200 Hz, Equation 3 suggests Q 200 of 311 ± 8 and 15 ± 2, respectively, in these intervals.

v S and Poisson's Ratio
Poisson's ratio describes how a material deforms when under stress, specifically as the ratio of transverse extension to axial compression. In ice, Poisson's ratio is typically 0.3-0.34 (Köhler et al., 2019) but is sensitive to the strength of anisotropy (Nanthikesan & Shyam Sunder, 1994) and increases weakly with temperature (Weeks & Assur, 1967). The dynamic Poisson's ratio, ν d , is defined using v P and v S as The variation of v P and v S (Figure 3c) is measured from the offset VSP following the same cross-correlation procedure applied to the zero-offset data. However, since direct waves in the offset VSP arrive at an oblique angle θ to the cable axis, path lengths are modified by cosθ. This neglects refraction effects (i.e., assumes straight rays), but these are likely small given the constant vertical velocity through the shallow ice column.
P wave velocities are consistent with those in Figure 3a and, although noisier, the same velocity variation is observed in the lowermost 200 m. Within the 100-500 m depth range that strong shear arrivals are recorded (Figure 2c), the mean v S is 1,815 ± 46 m s −1 . There is some tendency for v S to increase with depth, although this could be related to the straight-ray assumption, and a more sophisticated ray tracing algorithm will be applied in future analyses. Figure 3d shows two indicative measures of Poisson's ratio. The red curve is derived from the 50 m v P and v S trends and is therefore restricted to where the two have mutual coverage. This approach predicts the local ν d within the correlation gap around the borehole; however, the values it yields are lower than would be anticipated and are considered unreliable. Errors are compounded here since four straight-ray assumptions (two for each velocity estimate) are invoked for each ν d estimate. Instead, the orange curve is the path-averaged Poisson's ratio; neglecting refraction effects, direct P and SV waves follow the same path hence the ratio v P /v S is equivalent to the ratio of SV-to-P first-break traveltimes. The depth coverage extends beyond that of the correlation-based measure since direct-wave arrivals can be identified even where their signal-to-noise ratio is too weak for reliable correlation. This estimate of Poisson's ratio appears to be stable with depth: ν d = 0.338 ± 0.001.

Thickness Modeling of Subglacial Layering
Although the P wave velocity and thickness of layering beneath the borehole cannot be determined uniquely, permissible v P :thickness pairs satisfy a 44 ms two-way traveltime (Figure 2d). However, the deepest depth at which the reflection can reliably be perceived is 1,006 m, meaning 40 m of one-way propagation must be in basal ice. The traveltime contribution of the ice leg, t ice , is where h is the distance above the glacier bed, v ice is velocity through ice, and i is the angle of incidence to the vertical (Figure 4a). For the offset VSP geometry, i is~17°, and the basal v ice is 3,730 m s −1 (Figure 3). Therefore, Equation 5 shows that 22 ms of the 44 ms round trip comprises propagation through ice, with the remaining 22 ms being traveltime through the subglacial layer, t sgl . The thickness of this layer, h sgl , is where v sgl is the subglacial layer velocity and j is the refracted angle into the subglacial layer (from Snell's law, j = sin -1 (v sgl sin(i)/ v ice )). Equation 6 assumes a horizontal subglacial reflector, and negligible difference in incidence angles for the direct wave and the downgoing reflection leg (~0.6°for our geometry). The blue curve in Figure 4b shows pairs of v sgl and h sgl that satisfy t sgl = 22 ms. Hofstede et al. (2018) gave v P of 1; 839 þ1618 −94 m s −1 for subglacial material beneath Site S30; assuming the same velocity at L028, h sgl is 20 þ17 −2 m.

VSP Analyses
VSPs suggest little variation in either P or SV wave velocities within the majority of the L028 ice column, supported by the consistency of the depth-averaged Poisson's ratio. However, v P reaches 3,980 ± 120 m s −1 between 880 and 950 m depth. This velocity is consistent with that observed by Diez et al. (2015), for vertical P wave propagation through anisotropic fabric with a narrow cone opening angle. At Site S30, Hofstede et al. (2018) observed an anisotropic transition in the lowermost~80 m of Store Glacier (~84% ice thickness), consistent with enhanced deformation recorded in the lowermost 100 m of the ice column . This transition interpreted as the interface between Holocene-to Wisconsin-age ice (HWT). We therefore interpret the top of the fast-velocity layer as the HWT, here at 880 m depth. At~84% of ice thickness, this is consistent with a regional study of the HWT by Karlsson et al. (2013), who observed it at depths of 82-85%. Additionally, Horgan et al. (2008) identify pre-Holocene ice beyond~87% of ice thickness in Greenland's Jakobshavn Isbrae.
P wave attenuation rates are low throughout most of the ice column, implying efficient propagation through cold ice, but increase significantly below~900 m. Here, we measure Q 200 of 15 ± 2, similar to values measured by Peters et al. (2012) for basal temperate ice (approximately −1°C) in Jakobshavn Isbrae. We therefore predict similar conditions in the lowermost ice of Store Glacier. However, liquid water within temperate ice should also decrease v P (e.g., Endres et al., 2009), and a velocity reduction is only observed from~950 m depth. We therefore suggest that the cold-temperate transition occurs at the approximate depth of~950 m in BH19c. DTS observations made in the same borehole will verify this depth in future work.
Subglacial reflections were observed, originating beyond the borehole termination from the base of a layer that is 20 þ17 −2 m thick. Since the glacier bed has low reflectivity, the acoustic impedance of subglacial material must balance that of the basal ice (~3.4 × 10 6 kg m −2 s −1 ). With v sgl = 1,839 m s −1 , the impedance contrast would be removed if the density of subglacial material was~1,860 kg m −3 . Considering the elastic properties tabulated in Christianson et al. (2014), this density implies that Site L028 is underlain by consolidated, but neither deforming nor lithified, sediment.

DAS Applicability
DAS methods are clearly beneficial for characterizing the seismic properties of ice masses at high resolution. Nonetheless, vertically orientated fiber-optic cables are limited by their reduced sensitivity to oblique particle motion; a well-coupled three-component geophone would record all senses of particle motion, including SH waves that are undetectable with DAS. However, the ability of DAS to reconstruct independent seismic responses every~10 m throughout the ice column can outweigh these shortcomings, particularly if the goal of the survey is to benchmark P wave properties for conventional surface seismic analysis.
Our DAS cable has limited applicability as a horizontally orientated receiver in a conventional P wave reflection survey. Basal reflections would arrive with near-vertical incidence at the surface, hence with particle

10.1029/2020GL088148
Geophysical Research Letters motion orthogonal to the cable axis. Considering Figures 2c and 3c, in which P wave arrivals appear undetectable beyond~70°incidence, reflections from a 1 km-deep target would only be detected at offsets >600 m. However, the configuration has valuable applicability as an S wave receiver, recording, for example, source-generated S wave energy, reflected P-to-S conversions from the glacier bed , and the S wave components of natural glacier seismicity. The recent development of helically wound fiber-optic cables (Kuvshinov, 2015), sensitive to all directions of particle motion, could circumvent this limitation and thereby broaden DAS potential.
DAS interrogators are currently expensive to deploy, but their cost can be offset against enhanced survey efficiency. Once installed, the DAS cable allows field effort to be directed entirely toward source deployment; by contrast, conventional seismic surveys would require continual effort for the movement and operation of both source and receiver systems. This is particularly onerous for borehole surveying, where there is significant risk of instrumentation becoming stuck, thus requiring rapid surveying and/or coarser spatial sampling (e.g., Diez et al., 2015) than is possible with DAS. As a further advantage, DAS cables can also be monitored for colocated DTS, allowing the integrated interpretation of both seismic and thermal data.

Conclusions
We have demonstrated the feasibility and value of borehole DAS surveying in glaciology. Our VSP acquisitions sample the seismic response, at 10 m vertical resolution, throughout the 1,043 m thickness of a fast-flowing glacier in Greenland. These experiments allowed us to (i) measure the variation of P and SV wave properties; (ii) detect transitions relating to ice crystal fabric and temperature regimes with unique seismic properties, in a basal zone encompassing the lowermost 160 m of the glacier; and (iii) identify a subglacial layer of consolidated sediment up to 40 m thick.
Extended data analyses will consider the azimuthal variation of the seismic response and the passive seismic observations around Site L028, combining colocated radio echo-sounding and DTS measurements. We believe that DAS technology will play an increasingly important role in glacier seismology as the cost of DAS systems decreases and their potential is recognized.

Data Availability Statement
The two VSP data sets are available in .sgy format from the figshare repository Booth et al. (2020), under license agreement CC BY 4.0.