Distributed Acoustic Sensing (DAS) for Natural Microseismicity Studies: A Case Study From Antarctica

Icequakes, microseismic earthquakes at glaciers, offer insights into the dynamics of ice sheets. For the first time in the Antarctic, we explore the use of fiber optic cables as Distributed Acoustic Sensors (DAS) as a new approach for monitoring basal icequakes. We present the use of DAS for studying icequakes as a case study for the application of DAS to microseismic datasets in other geological settings. Fiber was deployed on the ice surface at Rutford Ice Stream in two different configurations. We compare the performance of DAS with a conventional geophone network for: microseismic detection and location; resolving source and noise spectra; source mechanism inversion; and measuring anisotropic shear‐wave splitting parameters. Both DAS array geometries detect fewer events than the geophone array. However, DAS is superior to geophones for recording the microseism signal, suggesting the applicability of DAS for ambient noise interferometry. We also present the first full‐waveform source mechanism inversions using DAS anywhere, successfully showing the horizontal stick‐slip nature of the icequakes. In addition, we develop an approach to use a 2D DAS array geometry as an effective multi‐component sensor capable of accurately characterizing shear‐wave splitting due to the anisotropic ice fabric. Although our observations originate from a glacial environment, the methodology and implications of this work are relevant for employing DAS in other microseismic environments.

. DAS therefore has great potential for seismology since it can provide dense, sub-wavelength sampling of a seismic wavefield over a range of possibly hundreds of kilometer. Such sampling could provide a step-change in our understanding and observing capability of seismic processes.
Here we present a study of naturally occurring microseismicity at Rutford Ice Stream, Antarctica (see Figure 1a), using DAS surface arrays to investigate the potential of DAS for natural microseismicity studies more broadly. The study site was chosen for the following reasons. Rutford Ice Stream flows at a speed of 100s of meters per year (Minchew et al., 2018), with icequakes generated as ice slides over the underlying bed (Kufner et al., 2021;Smith, 2006;Smith et al., 2015). A glacier data set is particularly suitable for such an investigation since icequake waveforms typically have high Signal to Noise Ratios (SNR) compared with other microseismic environments (Roeoesli et al., 2016;Smith et al., 2015). Furthermore, the velocity structure is approximately homogeneous, and therefore much simpler than other settings. The ice also exhibits a strongly anisotropic fabric (Harland et al., 2013;Smith et al., 2017), providing an ideal medium for attempting to measure seismic anisotropy using DAS. We first present a framework for initial detection and location of microseismicity using DAS, before demonstrating how DAS could be used for interrogating source physics and path effects. Source physics can help us understand basal sliding, while shear-wave splitting due to anisotropic path effects provides information associated with flow and deformation within the ice column. At each stage, we compare our results with conventional geophone data, quantifying the benefits, limitations, and factors to consider in future deployments of DAS for studying natural microseismicity.
HUDSON ET AL.

Data and Site Location
DAS systems measure the strain-rate along an optical fiber by sending a finite-duration pulse of light from a laser within a DAS interrogator along the fiber, as in Figure 1b. As photons travel along the cable, some undergo Rayleigh scattering from elastic collisions with particles in the fiber. As the cable deforms, for example, in the presence of a seismic wavefield, the position of the collisions relative to the end of the cable changes, resulting in a change in the two-way travel-time of the photons scattered back to the source (see Figure 1b, "In the fiber"). This is observed as a modulation in the phase of the returning light. If the change in length of a section of fiber through time can be measured by this phase modulation (see Figure 1b, "at the interrogator"), then the axial strain-rate can be calculated. This technique is called optical time-domain reflectometry (Masoudi & Newson, 2016;Zhan, 2019). Another important DAS concept is the gauge-length of the system, the length scale over which a change in strain is measured. The local change in strain is found by measuring the phase difference in the backscattered light from two closely separated points on the fiber. This measurement is proportional to the overall change in strain between these two points, the distance between which is referred to as the gauge length. It is controlled by the duration between two red pulses of light in Figure 1b, "At the interrogator," multiplied by the velocity of light in the fiber. Therefore, the gauge length limits the spatial resolution of the system, governing the response to different frequency signals (Dean et al., 2017). Here we use a Silixa iDAS™ system (Parker et al., 2014) with a channel spacing of 1 m and a gauge-length of 10 m (see Table S1 for a more complete interrogator specification). This effectively acts as a moving average measurement of strain over 10 m of fiber, in 1 m increments, with the maximum resolvable frequency signals for ice therefore being ∼200 and ∼90 Hz, for P-and S-waves with velocities of 3,841 and 1,970 m s −1 , respectively. These frequencies are greater than typical observed icequake corner frequencies at Rutford Ice Stream (Hudson et al., 2020;Kufner et al., 2021;Smith et al., 2015).
The seismic data were acquired in January 2020, during the austral summer. The deployment consisted of a Silixa iDAS interrogator with a 1 km fiberoptic cable (see Supplementary Information), as well as 16 4.5 Hz three-component geophones with Reftek RT130 dataloggers. Three of these geophones were colocated with the fiber. The interrogator was housed in a tent to keep it dry and within the operational temperature range. It was powered by a generator coincident with the tent. The full specification of the interrogator and fiber can be found in Tables S1 and S2. The fiber was deployed in several horizontal geometric arrangements, including a linear arrangement and a triangle, as shown in Figures 1 and 4. The line was 1 km in length and the triangle had sides of equal lengths of ∼330 m. The line data correspond to a 6-h period from 0100 to 0700 UTC on January 14, 2020, and the triangle to a 6-h period from 0100 to 0700 UTC on January 17, 2020. The interrogator was located at the NE end of the line. To investigate coupling, we also deployed the cable in buried and unburied configurations over separate periods of time. The fiber was buried by laying it in skidoo tracks and backfilling the track with snow at a depth of 5-10 cm (see Figure S1). The sampling rate of both the geophones and the DAS is 1,000 Hz.
The local magnitudes of icequakes in this study range from −1.9 to −0.9, calculated using the geophones. An example icequake arriving at a geophone and one channel of the DAS fiber is shown in Figure 2. This icequake is the same event as presented in the main source mechanism analysis later in this study. The SNR of the geophones is significantly higher than a single DAS channel for our data. However, the methods we present in this study generally make use of the dense-spatial sampling of the DAS data to attempt to minimize the effect of this potential limitation.
HUDSON ET AL.  Example icequake (04:01:42, January 14, 2020) arriving at a geophone coincident with the fiber (R104) and a single Distributed Acoustic Sensing channel 950 m along the fiber. Data here are bandpass filtered between 10 and 120 Hz. Phase arrivals for the geophone are indicated as red for the P-wave and blue for the S-wave. Note the shearwave splitting as indicated by the two S-wave phase picks on the horizontal components.

Microseismic Detection Method
To detect the icequakes we apply a waveform-based migration method called QuakeMigrate (Hudson et al., 2019;Smith et al., 2020), which has previously been successfully used to detect microseismicity in a range of settings, including at glaciers (Hudson et al., 2019). QuakeMigrate approximates the seismic energy associated with a seismic phase arrival at a particular receiver as an onset function, which in our case is defined as a continuous STA/LTA function through time. Onset functions for each receiver are stacked and back-migrated through time and space, in order to search for a coalescence of seismic energy corresponding to an event. One particular strength of this method is that although incoherent noise is back-migrated, it will not coalesce, therefore reducing the possibility of false detections. Furthermore, the initial back migration does not require phase picks to be made, which is often challenging given the SNR of DAS data (for example, see Figure 2). Phase-picking is only attempted for data where a peak in coalescence is detected. Another advantage of the method is that data from multiple instrument types can be combined once the waveform observations are approximated by onset functions, allowing for DAS and geophone time series data to potentially be used in combination in the detection and location algorithm. A further advantage of the QuakeMigrate algorithm is that it is open source, therefore allowing others to apply the methods demonstrated in this paper for other studies involving DAS, and hopefully improve upon our methods in the future.
Once the events have been detected and initially located, we refine the event locations using the non-linear earthquake relocation software, NonLinLoc (Lomax & Virieux, 2000). We refine event locations using NonLinLoc because it provides quantification of the statistical spatial and temporal uncertainty of the icequakes, allowing us to quantify the performance of DAS only versus geophone only versus combined network detection and location. The parameters used for the detection and location are given in Table S3.

DAS Source Mechanism Inversion Method
An earthquake source mechanism inversion involves using the observed radiation pattern of an earthquake to constrain the physical elastic failure mechanism that generates the earthquake. The earthquake source mechanism inversion used in this study is based upon the method described in Hudson et al. (2020) and available as the open source package SeisSrcInv (Hudson, 2020a(Hudson, , 2020b. Here we summarize the method, as applied to DAS data in this study, since a number of modifications are required. This specific DAS source inversion workflow is now implemented in SeisSrcInv. The source inversion method is a full-waveform Bayesian source mechanism inversion, randomly sampling the model space millions of times in order to obtain an estimate of the posterior probability distribution. We constrain the source model to be a Double-Couple (DC) model, which is appropriate for the predominantly stick-slip seismicity observed at Rutford Ice Stream (Hudson et al., 2020;Kufner et al., 2021;Smith et al., 2015). The DAS source inversion workflow is as follows:

Spatially Downsample the DAS Data
Due to the computational expense associated with calculating Green's functions and performing the source inversion, we only use every 10th channel, approximately equal to the gauge length (10 m), along the DAS cable to approximately sample the wavefield. This is found to be a sufficient resolution of spatial sampling for our case. Alternatively one could improve the SNR by stacking all these channels together, if they could be correctly aligned.

Filter the DAS Data
We first apply a rectangular fk-filter passing wavenumbers, k, of −0.04 to 1 0.04 m and frequencies between −150 to 150 Hz. To remove surface wave noise from a generator, we also apply notch filters centered at 33 and 66 Hz, with a bandwidth of 2.5 Hz. Alternatively, a more sophisticated fk filter informed by apparent surface wave velocities would likely be capable of removing the generator noise.

Generate Synthetic Modeled Waveforms
We generate synthetic modeled waveforms, for comparison to the observed data, for each moment tensor component ( ,  ,  , , , xx yy zz xy xz yz m m m m m m ), for each DAS channel, for an event. These are calculated using the program fk, which uses the Thompson-Haskell propagator matrix method (Zhu & Rivera, 2002). They are modeled for an isotropic, homogeneous ice medium overlaid with a 100 m firn layer (Smith, 1997) of decreasing velocity toward the surface.

Account for Anisotropy
Rutford Ice Stream has a strong elastically anisotropic fabric (Harland et al., 2013;Smith et al., 2017), causing Shear Wave Splitting (SWS) that has to be accounted for. This can be done either by applying a linearization and time shift correction to the observed data, or by simulating the effect of SWS on the isotropic modeled S wave phases to produce anisotropic synthetic waveforms. Correcting for such effects is valid as long as the anisotropy is a path effect, and the source region can be assumed to be isotropic. We implement the latter method, approximating the effect of SWS on the synthetic DAS data by applying an average anisotropic splitting angle and a fast-slow S-wave delay time. We assume that the S-waves arrive at approximately normal incidence (vertical) to the surface, due to the firn velocity structure. We can then approximate the fast and slow S-wave arrivals in the North and East axes from the LQT coordinate system of the synthetics using the equations, where  is the azimuthal angle from the source to the receiver and  is the average anisotropic splitting angle. These can then be combined into single N and E traces using, where  t is the fast-slow S-wave delay time.

Rotate Modeled Data Into DAS Axis
Once the simulated anisotropy has been applied to the synthetic modeled data, the North and East model components can then be rotated into the DAS axis, where  is the angle of the DAS fiber clockwise from North. It is important that this angle corresponds to the positive strain-rate direction along the fiber, as defined by the DAS interrogator.

Convert Data Into Same Units
The observed DAS data is in units of strain-rate, and the modeled DAS data is in units of particle velocity. In order to compare the modeled data to the observations, one therefore has to convert all the data either into strain-rates or particle velocities. We opt for converting the synthetic modeled data into strain-rate. The axial strain-rate,  x , the native measurement of DAS, by differentiating spatially, as given by, where x v is the velocity and x is the distance in the direction parallel to the cable. To make the modeled strain-rate consistent with the DAS, this differentiation should be applied over a length scale equal to the gauge length.

Perform Source Mechanism Inversion
Now that the observed and modeled data are in the same coordinate system and both in units of axial strain-rate, they can be compared to one another in a source mechanism inversion. We do this for a DC-constrained inversion with 1 × 10 6 samples.

DAS Shear-Wave Splitting Inversion Method
The strongly anisotropic nature of Rutford Ice Stream, as observed by shear-wave splitting (Harland et al., 2013;Smith et al., 2017), prompts the question of whether DAS can be used to observe seismic anisotropy. Shear-wave splitting provides a measure of the anisotropy along the raypath between the source and the receiver, which can be characterized by a delay time, dt, between the fast and slow S-wave arrivals, and the polarization, ɸ f , of the fast S-wave (e.g., Wuestefeld et al., 2010). SWS is typically estimated using three-component particle motion analysis on geophones. However, such an approach is not possible with linear DAS data because of its single-component nature and because it measures strain-rate rather than particle motion. Shear-wave splitting has previously been observed and analyzed using DAS in strongly anisotropic shales , where S-waves arrive as two distinct arrivals allowing dt to be measured. In that example, the anisotropy was known to have a relatively simple Vertical Transverse Isotropy (VTI) symmetry fabric such that determining the polarization was straightforward. However, at Rutford Ice Stream, the anisotropy has a more complicated orthorhombic symmetry (Smith et al., 2017), making determination of the polarization using DAS challenging.
Here we describe a methodology for estimating shear-wave splitting using the triangular DAS array. The motivation for this is to provide a proof of concept demonstrating that a 2D DAS geometry can be effectively used as a multi-component sensor capable of measuring shear-wave splitting. Using the triangular array partially alleviates for the inherent single component nature of DAS fiber because it records strain in a 2D plane rather than a 1D line, albeit with measurements at different orientations not at precisely the same location. However, if we assume that at the scale of the array the S-waves can be approximated as plane waves, we can approximate the triangular array as a point sensor by correcting for the spatial distribution, similar to the approach of Innanen et al. (2019). We then need to determine how any amplitude changes or polarity reversals recorded on the three sides of the array relate to the polarization of the S waves. The strain sensitivity pattern of an S-wave depends on both the orientation of the ray slowness vector (i.e., wavefront propagation direction) and of the polarization vector Benioff, 1935;Karrenbach et al., 2019). Thus, if we can estimate the orientation of the slowness vector, then we can invert for the polarization that best fits the observed data.
For the splitting analysis we use a simplified velocity model with a constant firn layer S velocity of 1,456 m s −1 (Smith et al., 2015), rather than a velocity gradient. We also assume an isotropic firn layer, such that the majority of the observed splitting is accrued as the wave passes through the ice column, as in Smith et al. (2017).
The anisotropy inversion processing steps to find the fast-slow S-wave delay-time and polarization are as follows: 1. The triangular array is divided into the three linear segments, before applying a slant-stacking processing technique to stack data over various possible linear velocities represented by apparent slowness values. These are normalized by the number of channels in the stack to preserve amplitudes. This single processing step: (a) provides an estimate of one of the components of the slowness vector, which is required to estimate the propagation direction in order to forward model the strain sensitivity; (b) reduces the varying travel times over the length of the array to a single point measurement at its midpoint; and (c) aids the identification and picking of the two S phases by separating them by their apparent slowness. 2. The fast and slow S-wave arrivals are then picked in the slowness-time domain, with the peak amplitude of the phase immediately after the pick taken as the recorded strain-rate amplitude of that arrival. The travel time difference between the fast and slow S-waves provides an estimate of the shear-wave splitting delay time dt. This analysis is then repeated for each side of the triangle, with the estimated dt at the midpoint of the array taken as the average of all those recorded on each side. 3. We then determine orientation of the slowness vectors (i.e., the wavefront propagation directions) of the two S waves. Since we have measured the apparent slowness of each S wave on the three sides of the array, we can apply a least squares inversion to solve for the best fitting horizontal slowness components. One should note that this slowness method provides an independent estimate for the azimuth of the source epicenter relative to the array from the hypocentral locations derived earlier in this study. This independence from our detection and location method is important because our original locations could be biased by the assumption of no anisotropy. In addition, we need to estimate the inclination of the propagation through the firn layer (to estimate emergence angle) and through the ice column (to define the path over which most of the splitting was accrued). We estimate this using Snell's law with the estimated horizontal slowness and a simplified two-layer velocity model. 4. Once we have determined the S-wave slowness vectors, we can then perform an inversion to obtain the fast and slow S-wave polarizations. The strain tensor, e S, associated with an S-wave propagating in the x 1 direction and polarized in the x 3 direction has the form, where the scalar factor e S determines the overall amplitude of the strain. Thus, to simulate the DAS response to an S-wave with arbitrary propagation and polarization directions, we rotate this tensor into the appropriate orientation, and project the resulting tensor onto the fiber geometry. To eliminate the need to solve for e S , we instead model the ratio of strain on different sides of the array, relative to a reference side (for example, A side 1 /A side 2 and A side 3 /A side 2 if side 2 is the reference). Since we have already estimated the azimuth and emergence inclination in step 3, this then leaves only the polarization to be determined, which is described by a rotation about the propagation axis. We then apply a least squares inversion to solve for the S polarization that minimizes the misfit between the modeled and observed amplitude ratios. We independently invert for the fast and slow S-waves.
In order to compare our splitting measurements with previously published splitting results at Rutford from Smith et al. (2017), we convert dt to a percentage difference in fast and slow S velocities dv S , by using the formula from (Wuestefeld et al., 2010), where r is the source-receiver distance and v Smean is the mean S-wave velocity along the full path.

Source and Noise Spectra
Earthquake spectra can provide insight into both the source physics and ambient noise levels that may hamper the detection of microseismicity. Figure 3a shows the spectra of a 1.5 s event window and two 10 s noise windows for the linear DAS configuration. Figure 3b shows the same time periods, but for a geophone collocated at one end of the DAS fiber.
The two noise periods correspond to times when the fiber-optic cable is unburied and buried a few cm deep, respectively (see Figure S1). The DAS spectra above 1 Hz in Figure 3a show that the noise in the unburied DAS data is more than an order of magnitude greater than when it is buried. Geophone data in Figure 3b, confirms that noise conditions are approximately constant over both periods. The increased noise present in the unburied DAS data is interpreted to be due to wind incident on the cable, which is a source of broadband noise and cable resonance with fundamental and higher order modes at 22 and 44 Hz, respectively. Burying the fiber therefore minimizes wind noise.
Paradoxically, the noise for the buried cable is greater below 1 Hz. Burying the cable increases the DAS response over this range by more than 15 dB. This is because the buried cable is coupled to the ice better, and is therefore more sensitive to the primary and secondary microseisms (Cessaro, 1994). Theoretically, the low-frequency DAS instrument response is only limited by the recording duration.
A further feature to note is the pair of minima in the DAS spectra at 33 and 66 Hz. These are a result of 2.5 Hz bandwidth notch filters applied to the data to remove the response of the DAS to surface waves produced by a generator used to power the DAS. This would obviously also affect a geophone collocated next to a generator too, although a generator is not required to power a geophone. DAS is sensitive to these anthropogenic surface waves, and although mitigating for such signals could be challenging, they could likely be removed with sophisticated fk-filters or other filters. Indeed, when working in remote areas, a power source for experimental equipment such as the DAS interrogator itself may be required. Whilst surface waves from the generator are detrimental to our study, since they mask icequake phase arrivals, these observations positively emphasise the potential of DAS for surface wave studies.
Figure 3 also includes the spectra of an event when the DAS cable is buried. No events are detected when the cable was unburied, and the data in Figure 3a demonstrate why, with the event energy well below the noise threshold of the unburied DAS. The event spectrum observed by the DAS has a higher SNR than the spectrum observed by the geophone. Although a dominant frequency of ∼70 Hz is visible in the geophone data, the signal is significantly higher above the noise level in the DAS data, even with a notch filter applied at 66 Hz. This higher spectral SNR is due to the DAS data being stacked over ∼1,000 channels. It is worth noting that our comparison between DAS and geophones is not strictly fair as stacking hundreds or thousands of geophones over a 1 km length would likely yield a higher SNR spectra than DAS. However, an inherent benefit of DAS is exactly this, in that it allows for thousands of channels to be measured simultaneously and at lower cost.
HUDSON ET AL.
10.1029/2020JB021493 8 of 19 The resolution, bandwidth and SNR of DAS spectra when buried suggests that DAS is promising for studying the spectra of microseismic sources. The frequency content of a microseismic source allows for limited interpretation of the event source-time function. However, more involved analyses, such as spectral measurement of magnitude (Butcher et al., 2020;Stork et al., 2014), require knowledge of the absolute amplitude of the spectra, in addition to the frequency response of the fiber. This is commonly referred to as the transfer function. Understanding the transfer function of both the interrogator and the fiber-medium coupling should therefore be a priority for studies utilizing DAS, given the potential of DAS for providing higher resolution, bandwidth, and SNR spectra than conventional geophones or broadband seismometers . Indeed, although the transfer function of the DAS is uncalibrated in this study, the optical properties of the fiber can be used to calculate the absolute strain-rate , and if the coupling of fiber to the medium were also known, then the fiber-medium coupling transfer function could then be determined.
The advantages of DAS for studying the lower frequency limit of the spectrum also suggest its applicability for ambient noise studies in glaciated environments, such as those undertaken using DAS in other settings Dou et al., 2017;Spica et al., 2020;Zeng et al., 2017). Furthermore, DAS may provide a means of observing other environmental processes with low frequency signals, such as gravity waves in ice shelves (Lipovsky, 2018).
The quality of the DAS spectra also prompts the question of whether a spectrum-based earthquake detection method (Hudson et al., 2019) would prove more effective than the QuakeMigrate time-series method. However, this would require averaging the spectrum over many DAS channels, which would reduce the spatial sampling, or increasing the sensitivity of the DAS system to facilitate higher SNR spectra from individual channels.

Detection and Location of Natural Microseismicity Using DAS
An important question for the applicability of DAS for studying natural microseismicity is whether one can detect and locate seismicity using DAS alone. To address this question, we compare results of icequakes detected using DAS and geophones separately.
Unlike the three-component geophone records, P-wave arrivals are not visible in our DAS data because only the horizonal component of strain is recorded along the fiber. A near-surface firn layer, of substantially lower seismic velocity, refracts P-waves toward vertical incidence. Surface DAS recordings in areas without a firn layer would yield P-wave observations . However, fast and slow shear-waves are visible in the DAS data. This is a diagnostic of seismic anisotropy in the ice column, which has been previously documented at Rutford Ice Stream (Harland et al., 2013;Smith et al., 2017) and Korff Ice Rise . For icequake detection and location, we assume an isotropic fabric, as QuakeMigrate can only pick a single S-wave arrival, typically the fast arrival. This will result in some uncertainty in hypocentral location if the slow S-wave arrival is picked instead.
Figures 4 and 5 show the icequake detection results. Only the hypocentral distance and azimuth are resolved adequately for our DAS geometries. The depths of the DAS-only results are therefore artificially constrained prior to detection and location to between 1,700 and 2,100 m bsl, based on previous Rutford icequake observations (Hudson et al., 2019(Hudson et al., , 2020Kufner et al., 2021;Smith et al., 2015). The geophone-only depths are not artificially constrained. The lateral spatial clustering observed in Figure 4 for all network configurations, interpreted to be sticky patches of the bed, is expected from icequake datasets (Hudson et al., 2019;Roeoesli et al., 2016;Smith et al., 2015;Winberry et al., 2009). It shows that both the DAS line and triangle geometries are able to resolve spatial clustering in the data, if depths are independently constrained. The addition of a vertical section of fiber (e.g., Lellouch et al., 2020Lellouch et al., , 2021 could also provide depth constraint in the vertical plane, in the same way that the triangle array breaks the symmetry, therefore providing better epicentral constraint in the horizontal plane. A second observation is that significantly fewer events are detected in data from the DAS-only configuration compared to the geophone-only configuration. Only 499 events are detected using the DAS triangle, compared to 1,321 geophone event detections, and only 139 events are detected using the DAS line, compared to 1,270 geophone event detections (see Figures 5a-5d). The line is therefore a less effective configuration than the triangle for event detection, likely partly due to the two-dimensional nature of the triangle configuration resulting in a higher peak coalescence of energy at a single location, rather than being split by the geometric ambiguity of the line, but also because it may record more S-wave energy due to its effectively quasi two-component measurement. The addition of a vertical section of fiber could also provide depth constraint in the vertical plane, in the same way that the triangle array breaks the symmetry, therefore providing better epicentral constraint in the horizontal plane. Such a vertical section would also have reduced surface wave noise compared to a horizontal section. A spatial detection bias is also apparent in the data, with icequakes from clusters beyond the SW end of the linear DAS configuration detected, but icequakes at closer offsets but perpendicular to the cable not detected. This sensitivity is likely a result of the single component nature of the DAS making it insensitive to certain S-phase polarizations, for the basal icequake radiation pattern orientations observed at Rutford Ice Stream (Hudson et al., 2020;Kufner et al., 2021;Smith et al., 2015).
We also compare the spatial and temporal uncertainty of detected events. Figures 5a and 5c show the epicentral uncertainty of events detected for the line and the triangle, respectively. We compare epicentral uncertainty, rather than hypocentral uncertainty, since the depths of events for the DAS-only detection are artificially constrained. While for some events the line and the triangle DAS-only epicentral uncertainties are similar to geophone-only detection uncertainties, there are a number of DAS-only events that have spatial and temporal uncertainties greater than the plotted range (i.e., epicentral uncertainty >0.5 km and origin-time uncertainty >0.1 s). In summary, the geophones provide better spatial constraint of the icequake epicenters and can also constrain icequake depth, while both the DAS and geophones appear to provide similar constraint on event origin times. There are insufficient events detected by the line to quantify whether the triangle or line has better spatial constraint. Figures 5b and 5d show the origin time uncertainty of the events. For both DAS geometries, the geophone-only and DAS-only data suggest similar constraint on the origin time uncertainty. Although the epicentral uncertainties using DAS presented here are large, DAS could still be particularly useful for monitoring applications where rates of seismicity, rather than accurate locations, are of primary importance.
These results suggest that DAS has limitations for microseismic detection. First, for this experiment arrangement with sources ∼2 km below surface, geophones are significantly better than DAS for detecting microseismicity. This is because: (a) the spatial extent of the geophone network is much greater than the DAS (see Figure 4), a limitation specific to our DAS deployment that could be overcome by deploying more fiber; (b) the geophones measure three components of ground motion, and so are sensitive to P-, SV-, and SH-phase arrivals; and (c) the geophones have a much higher SNR than single DAS channels. Although the firn velocity structure causes the DAS to be sensitive to only S-phases in this study, horizontally deployed DAS on ice without firn would be sensitive to both P-and S-phases, providing better event depth constraint . A second limitation is the complexity of combining DAS and geophone data together for detection and location. Theoretically, DAS and geophone data could be combined to reduce spatial and temporal uncertainty. However, there could be a trade-off between the gain of additional observations and detrimental additional noise. Weighting the combined data to mitigate for this is complex, and likely siteand network-geometry specific. We therefore do not investigate weighting methods in this study, instead HUDSON ET AL. , respectively, except for the triangle fiber deployment. Note that epicentral uncertainties are clipped at an upper limit of 0.5 km for clarity. There are events with uncertainties greater than this but these are not plotted.
suggesting this as an area for future work. Ideally one would also convert the time-series for each type of instrument into mutually consistent units, so as to avoid any arrival-time uncertainty due to possible phase shifts caused by the difference between a strain-rate and velocity measurement. Overall, the poor performance of our deployment for detection and location is likely dominated by the spatial extent of the DAS compared to the geophones. We suggest that if the horizontal spatial extent of the DAS is comparable to or greater the spatial extent than the depth of seismicity, then better performance may be achieved, such as in Walter et al. (2020).

Source Mechanism Inversion
Microseismic source mechanisms are often studied using conventional seismic networks. However, although theory and observations show that DAS has promise for source mechanism analysis Karrenbach & Cole, 2020;Vera Rodriguez & Wuestefeld, 2020), a full, quantitative source mechanism inversion has not previously been attempted for such data. Here, we show an example of a source mechanism inversion using DAS, and compare the result with that of an inversion using conventional seismic instrumentation for the same event. We use the full-waveform method described in Hudson et al. (2020), with further DAS specific considerations described in the Supplementary Information.
The source mechanism inversion results for the icequake highlighted by the yellow star in Figure 4 are shown in Figure 6. This icequake occurred when the DAS was configured in a linear arrangement. Source mechanism results for additional icequakes are provided in the Supplementary Information. Figure 6a shows the event arrival at every channel along the DAS cable. The fast shear wave is dominant at smaller epicentral distances, with the slow shear wave prevalent at greater distances. The north component of a collocated geophone at the SW end of the linear cable is shown in black, at its approximate location. Figure 6b shows the DAS source mechanism inversion result. It should be noted that although we use S-waves in the inversion, we plot the P-wave radiation pattern, for consistency with other studies. The double-couple (DC) upper hemisphere focal mechanism is as one might expect (Hudson et al., 2020;Kufner et al., 2021;Smith et al., 2015), showing the observed waveforms used in the inversion, the modeled result, and the difference, colored by normalized strain-rate. The blue scatter points on the focal mechanism plot indicate the fiber and the red solid arrow and dashed lines indicate the slip vector and its associated uncertainty, respectively. (c) DC constrained full-waveform source mechanism inversion result using geophone observations. The Z components include the P phase arrival only, and the R and T components include the S phase arrival only. Note that the Z and RT components are not plotted temporally aligned with one another, so as to clearly see P-and S-phase arrivals on the Z and RT components, respectively. Both focal mechanism radiation patterns are for an upper hemisphere projection for P-wave radiation, for consistency with other similar studies. Black quadrants correspond to compressional first arrivals. Both inversions used 10 6 samples of the moment tensor space.
(see Figure 4). Uncertainty in the slip vector, shown by the red dashed lines, is of the order of   20 . It should be noted that a P-wave radiation pattern is plotted for consistency with other studies, although the fiber, and hence the inversion, is actually only sensitive to S-wave arrivals. Due to the  45 rotation of the S-wave radiation pattern relative to the P radiation pattern, the DAS cable lies near an S-wave radiation pattern maximum. The fit between the model and the observations is quantified by a variance reduction value of 0.65. The 2D observation, model, and difference fields show that the model provides a good fit, with negligible discernible coherent energy shown in the difference field for at least the majority of the slow S-wave. The scatter in the best fitting model result is due to the automatic alignment algorithm not correctly aligning every individual channel to the data. However, the general trend and the small proportion of scatter relative to fitted signal provides us with confidence that the inversion is successful. To our knowledge this is the first microseismic source mechanism inversion performed using DAS observations.
For comparison, we also show DC-constrained moment tensor source inversion results using geophones instead of DAS. All P-phase polarities are correctly fitted by the geophones, as are the majority of S-phase polarities, although there is more uncertainty for the S-wave matches due to some of the SWS not being entirely compensated for. Unfortunately, geophone coverage of the focal sphere for the network configuration is not as well configured for source mechanism analysis as previous studies (Hudson et al., 2020;Kufner et al., 2021;Smith et al., 2015), leaving the most likely DC source poorly constrained compared to previous observations. While the slip vector is approximately the direction of ice flow, and therefore in agreement with the DAS source inversion, the uncertainty, denoted by the dashed lines, indicates   270 azimuthal uncertainty. This is significantly greater than the   45 azimuthal uncertainty for the DAS source inversion. These results show that DAS provides better constraint of the source mechanism than the geophone network configuration, at least for this icequake.
Although the DAS source mechanism inversion outperforms the geophone inversion here, there are a number of challenges and limitations of the inversion. Indeed, some other DAS-constrained source mechanism inversions of other icequakes in this data set (see Figure S2) have large slip-vector uncertainties compared to the icequake in Figure 6. The first limitation is that due to the low velocity firn layer, body wave phases arrive approximately vertically. This means that DAS is not sensitive to P-wave arrivals that have particle motions perpendicular to the fiber, unlike studies with no firn layer , and so the source mechanism can only be constrained using S-waves. The SNR of the DAS strain-rate data also limits the performance of our source inversion method. This is particularly evident for the icequakes inversion results in Figure S2. Comparing DAS strain-rate observations to model outputs can also be challenging. While some models output strain-rate, the results from the wave propagation code used here, fk (Zhu & Rivera, 2002), have to be converted from particle velocities. Even if the model used did output strain-rate directly, one still needs to simulate gauge-length effects. Another challenge is accounting for any anisotropy at the study site. Our results suggest that we have adequately accounted for anisotropy, but our method only holds for anisotropy path effects that can be approximated as homogenous on the spatial scale of the cable. If anisotropy were present at the source, or varied and was unknown across the length of the cable, then our method would be invalid. However, one could instead view the clear anisotropy observations in Figure 6a as an advantage of DAS, as an inversion for ice fabric anisotropy is theoretically possible. There are also further limitations due to the inversion method itself. Due to uncertainty in the source location and velocity model, the modeled waveforms may not perfectly align with the observed phase arrivals. We therefore have to use an automatic alignment algorithm, described in Section 2.3. This does not always align with the observed phase arrival but instead on noise at a different time. This behavior is likely to be primarily due to the low SNR of individual DAS channels, but could also be partly caused by inhomogeneous attenuation of the medium and/or assumptions related to the frequency content of the source-time function. Given these limitations, we therefore present these source mechanism inversion results as a foundation for future analysis of higher SNR data associated with more earthquakes, possibly in combination with more advanced inversion methods yet to be developed.
Here, we have treated the DAS and geophone data separately for the source mechanism inversions to provide a clear comparison between DAS and geophones. However, it is theoretically possible to combine DAS and geophone observations in a joint inversion. While this would not provide substantial gains for the event discussed here due to little increase in spatial constraint over the focal sphere, it might in other situations.
One could imagine constraining the inversion better by deploying DAS within the center of a network and geophones located more sparsely and at greater distances, for example, Performing such a joint inversion would provide its own challenges, such as how to weight DAS and geophone observations, since SNR and the single-component versus three-components may introduce bias into the solution.

Seismic Shear-Wave Splitting Inversion
We demonstrate how DAS observations can be used to quantify shear-wave splitting, and hence seismic anisotropy of the ice column, using an icequake observed with the triangular array geometry. These results are shown in Figures 7 and 8 Figure 7 shows how the anisotropy-inversion methodology is applied to measure the slowness of the icequake S-wave arrivals. The icequake clearly shows the differing moveout observed on each side of the array (see Figure 7a). The slant-stacking result is shown in the time-slowness domain in Figure 7c. The fast and slow S-wave arrivals can then be identified in this domain, as well as the peak strain-rate amplitude. To find the fast and slow S-wave polarizations, we choose side two of the array as our reference (see Figure 7b) and normalize the recorded strain amplitudes on each side by the amplitude on side two. These amplitude ratios allow us to determine the polarization, and are plotted in Figures 8a and 8c for the fast and slow S-waves, respectively.
The S-wave polarization inversion results for the fast and slow S-wave for this example icequake are shown in Figure 8, along with the modeled strain sensitivity pattern of the best fitting inversion results. The best fit solution for the fast S-wave polarization, ɸ f , is −74°, where ɸ f is defined in the ray frame according to the convention of Wuestefeld et al. (2010) (i.e., ɸ f is 0° for SV and ±90° for SH, where positive angles are measured clockwise looking along the ray toward the source). The best-fit solution for the slow S-wave polarization is 15°. The observation that the fast and slow S-wave polarizations, ɸ f , are approximately perpendicular gives us confidence that this method is correctly measuring SWS.
HUDSON ET AL. Overall, this example event is found to have a shear-wave splitting delay time of 0.0093 s, based on averaging all three vertices of the triangular array to estimate the splitting at the center of the array. We find a best fitting propagation azimuth of 72°, as indicated by the arrow in Figure 7b, with horizontal slownesses of 0.28 and 0.26 s/km for the fast and slow waves, respectively. Furthermore, we estimate a ray emergence inclination of ∼22°-25° based on an estimate of a mean firn layer S-wave velocity of 1,456 m s −1 (Smith et al., 2015). However, we expect that the true inclination may be steeper than this, due to a nonlinear gradient in the firn layer velocity structure. This is evidenced by the lack of P wave sensitivity. We find that reducing the inclination modulates the magnitude of the strain recorded but has only a minor effect on the relative strain sensitivity patterns observed, so we do not believe that uncertainty in emergence angle is a significant source of error.
We compare this splitting measurement with that of a previous study at Rutford (Smith et al., 2017) by converting our splitting fast-slow delay time into a S-wave velocity difference, dv s . From the estimated source location for the event and the approximate path-length between source and receiver, we estimate a dv s of 0.78% with an average inclination of 29°. Although Smith et al. (2017) found some splitting as high as 4%-5%, they also found that dv S and ɸ f varied substantially for different azimuths and inclinations, due to the orthorhombic symmetry of the anisotropy. If we restrict comparison of our results to only splitting parameters from Smith et al. (2017) with ray paths within azimuthal and inclination ranges of ±15° of our measurement, they have a mean dv S and ɸ f of 0.89% and −75°, respectively. These are remarkably similar to our result.
The consistency between our result and previously published results gives us confidence in our approach, suggesting that smaller 2D DAS arrays can be effectively used as multicomponent sensors to detect and quantify SWS. Practically, however, this is unlikely to become the preferred method of characterizing SWS. Deploying a DAS system in a small, azimuthally varying 2D array removes one of the major benefits of DAS systems: their large aperture, instead effectively reducing the array to a point sensor. In our case, the source-receiver geometry is such that the raypath sampled a relatively weak splitting axis, producing a splitting measurement not very representative of the broader anisotropic fabric. In contrast, a less dense but broader aperture geophone array could record multiple splitting measurements from a single event, covering a wide range of azimuths and inclinations, which would provide a more comprehensive picture of the anisotropy. Nevertheless, we suggest there could be benefits of including small azimuthally varying 2D segments of DAS arrays as part of a larger aperture deployment to introduce some multi-component sensitivity.

Considerations for Future DAS Deployments
Current limitations of DAS that constrained our experiment include practical limitations on the spatial extent of cables, the single component nature of the strain measurement, and the limited amplitude sensitivity leading to low SNR on individual channels. We suggest the following considerations for future surface DAS experiments based on lessons learned from this experiment. Many of these suggestions are not new concepts but are nonetheless useful for planning any DAS deployment. (a) Consider hybrid geophone and DAS networks for microseismic detection. (b) Use cable lengths comparable to or greater than the source depths to be studied. This might be logistically challenging for remote sites. (c) Orient the cable in multiple directions, or use multiple cables with time-synchronized interrogators. Breaking the symmetry of a DAS system can improve detection at different azimuths. (d) Isolate the cable from wind-shear and other surface noise. (e) Consider the wavelength of seismic sources being studied, ensuring that the gauge-length is less than half the wavelength, for adequate Nyquist sampling. For example, our study, basal icequake S waves typically have wavelengths of ∼25 m, which is greater than twice the gauge-length. (6) Consider deploying in boreholes (e.g., Lellouch et al., 2020Lellouch et al., , 2021 in combination with surface arrays. Advantages of this would include being sensitive to P waves and increased SNR. However, deployment in a single, vertical borehole would result in complete azimuthal ambiguity in source location. Ideally, one would deploy DAS horizontally and vertically for increased sensitivity and hypocentral constraint. As well as the aforementioned considerations for utilizing current DAS technology, there are also technological developments that would improve DAS for natural microseismic study applications. First, the gauge-length imposes a limit on the ability to observe high frequency, near source observations. The ability to vary gauge-length during an experiment, or use of a chirp-style signal to instantaneously measure multiple gauge-lengths may prove useful for certain investigations. A more significant improvement might be the development of readily available helically wound fiber to measure strain in three dimensions (Baird, 2020;Lim Chen Ning & Sava, 2018). This may allow for greater sensitivity of P-wave observations for various cable orientations or shallow velocity structures. However, helically wound fiber would conversely affect the sensitivity to S-waves, possibly negating any advantages of enhanced P-wave sensitivity. Practically, DAS interrogators currently consume considerable power and the interrogator is expensive compared to geophones or standard seismometers. Reducing power consumption and unit cost would allow DAS to be deployed in more remote locations and for longer periods of time than in this study. However, even without such improvements to DAS interrogators, the fiber itself is relatively inexpensive, allowing for intermittent long-term monitoring of a site.

Conclusions
DAS provides dense spatial sampling of seismic wavefields and therefore has significant potential for studying natural microseismicity. We investigate the potential of DAS deployed on the Earth's surface. We show that DAS alone is relatively poor at detecting and locating seismicity compared to geophone observations, when deployed with an aperture less than the depth of the earthquake hypocenters being studied. However, DAS can outperform geophones in other analysis of natural microseismic sources. The SNR and bandwidth of the spectrum measured by stacking a number of DAS channels is significantly improved over a single geophone, and likely seismometers, especially for low, near quasi-static frequencies. Obtaining the DAS transfer function would allow the use of spectral observations to constrain moment magnitude of microseismicity. We also demonstrate that the dense spatial sampling that DAS provides can constrain source mechanism inversions better than conventional geophone networks, at least for our specific source-receiver geometry. Finally, we show that for a 2D fiber geometry, it is possible to quantify the shear-wave-splitting exhibited due to anisotropy in the ice column. In summary, whilst DAS has significant potential for natural microseismicity studies, it also has limitations. Any future study should consider utilizing the strengths of DAS to address the underlying science question, but also include conventional seismometers to compensate for the limitations of DAS, in a complementary, hybrid deployment.

Data Availability Statement
The seismic data will be made available through IRIS and the UK NERC Polar Data Centre. The data are also made available through Zenodo (https://doi.org/10.5281/zenodo.4778368). The microseismic detection software used is QuakeMigrate, which is open source software that is deposited on Zenodo (Winder et al., 2021). The moment tensor source inversion python software used is SeisSrcInv, which is deposited on Zenodo (Hudson, 2020a). The moment magnitude calculation python code used in this study, SeisSrcMoment, is also deposited on Zenodo (Hudson, 2020b). GMT was used in the production of the map figures (Wessel et al., 2019). 299622), which is part of Accelerating CCS Technologies (ACT) program. We thank Silixa for the loan of an iDAS interrogator. Finally, we thank the reviewers A. Lellouch and E. Smith, as well as the editor N. Nakata, for such extensive feedback and comments that has undoubtedly improved the manuscript.