Array processing in cryoseismology: a comparison to network-based approaches at an Antarctic ice stream

Seismicity at glaciers, ice sheets, and ice shelves provides observational constraint on a number of glaciological processes. Detecting and locating this seismicity, specifically icequakes, is a necessary first step in studying processes such as basal slip, crevassing, imaging ice fabric, and iceberg calving, for example. Most glacier deployments to date use conventional seismic networks, comprised of seismometers distributed over the entire area of interest. However, smaller-aperture seismic arrays can also be used, which are typically sensitive to seismicity distal from the array footprint and require a smaller number of instruments. Here, we investigate the potential of arrays and array-processing methods to detect and locate subsurface microseismicity at glaciers, benchmarking performance against conventional seismic-network-based methods for an example at an Antarctic ice stream. We also provide an array-processing recipe for body-wave cryoseismology applications. Results from an array and a network deployed at Rutford Ice Stream, Antarctica, show that arrays and networks both have strengths and weaknesses. Arrays can detect icequakes from further distances, whereas networks outperform arrays in more comprehensive studies of a particular process due to greater hypocentral constraint within the network extent. We also gain new insights into seismic behaviour at the Rutford Ice Stream. The array detects basal icequakes in what was previously interpreted to be an aseismic region of the bed, as well as new icequake observations downstream and at the ice stream shear margins, where it would be challenging to deploy instruments. Finally, we make some practical recommendations for future array deployments at glaciers.

1 Introduction

Cryoseismology is a rapidly emerging field that shows promise for studying glaciological processes (Podolskiy and Walter, 2016; Aster and Winberry, 2017). For example, icequakes associated with basal slip (Smith et al., 2015; Roeoesli et al., 2016; Kufner et al., 2021) can provide observational constraint on frictional processes in ice dynamics models (Gräff and Walter, 2021; Hudson et al., 2023, 2020; Köpfli et al., 2022; Lipovsky et al., 2019; Zoet et al., 2013). However, glaciers are often challenging to access and logistically expensive to operate within, and potential seismically active areas are typically highly variable temporally and spatially. To date, detecting subsurface icequakes has generally been performed using conventional network-based detection and location methods. Conventional seismic networks typically require receiver spacings similar to the distance of the events from the receiver and are predominantly sensitive to icequakes within the spatial extent of the network. Therefore, adequately sampling a specific region may require tens to hundreds of receivers. Conversely, seismic arrays are predominantly sensitive to icequakes outside the array aperture, with individual arrays requiring significantly fewer instruments than a conventional network. Seismic arrays could therefore facilitate smaller seismic deployments than conventional networks while enabling event detection at greater distances, in a similar way to the gain provided by arrays for nuclear test ban monitoring (Bowers and Selby, 2009), for example. Here, we address the following question: “to what extent are array deployments useful for subsurface microseismic icequake studies?” The scope of this work is to investigate the advantages, challenges, and practical application of array-processing methods for icequake cryoseismology studies. Specifically, we assess the sensitivity of arrays compared to networks for icequake detection and icequake location. We also assess whether arrays can detect event types typically missed by standard network-based detection algorithms and discuss the limitations of arrays compared to networks. We also present an implementation of a frequency-domain array-processing method that is made available to the community to accelerate the application of array processing in cryoseismology.

Array processing has existed as a method for decades. The core component of array processing, beamforming, involves combining plane-wave arrivals from all receivers in an array to find the direction of the origin of the wave from the array (Rost and Thomas, 2002). In this study, we focus on plane-wave beamforming rather than on other similar methods such as matched field processing (Sergeant et al., 2020; Nanni et al., 2022), which are generally also viable alternatives to traditional network-based detection and location approaches but are beyond the scope of this study. Beamforming was originally developed for nuclear test ban monitoring (Bowers and Selby, 2009) but has since been applied to other topics, including studying the structure of the earth (Wolf et al., 2023; Wang and Vidale, 2022; Thomas et al., 2002), monitoring offshore seismicity (Jerkins et al., 2023), ambient seismic noise source analysis (Bowden et al., 2021; Löer et al., 2018), and detection of seismicity using distributed acoustic sensing instrumentation (Van Den Ende and Ampuero, 2021; Näsholm et al., 2022; Klaasen et al., 2021; Lellouch et al., 2020; Hudson et al., 2021). However, the use of array processing within cryoseismology is limited. Beamforming has been used to locate large slip events at Whillans Ice Stream, Antarctica (Pratt et al., 2014). Multiple regional arrays have been used to locate glacier earthquakes, likely caused by significant calving events (Ekström et al., 2003; Ekström, 2006; Tsai and Ekström, 2007), and to study calving processes locally (Köhler et al., 2015, 2016, 2019, 2022; Podolskiy et al., 2017). Array processing has also proven useful for locating glacial tremors (Lindner et al., 2020; Umlauft et al., 2021; McBrearty et al., 2020), locating near-surface icequakes associated with crevassing (Lindner et al., 2019), locating seismic sources on ice shelves (Hammer et al., 2015), and even measuring the thickness of sea ice (Serripierri et al., 2022). The closest applications to that investigated in this study, specifically targeted at basal microseismic icequake detection, are Cooley et al. (2019), who use beamforming event location but not detection, and Lindner et al. (2019), who detect and locate surface crevassing using beamforming. Here, we present results of icequake detection and location solely using array processing and comparing it with a conventional network-based approach, using a dataset from Rutford Ice Stream (RIS), Antarctica (see Fig. 1). The array-processing detection uses an array of 10 sensors shown in Fig. 1d, with the network-based comparison performed using a network of 16 sensors (see Fig. 1).

Figure 1 The seismic array deployment. (a) Location of Rutford Ice Stream (RIS) with respect to Antarctica (basal topography is from Fretwell et al., 2013). (b) Overview of the deployment within RIS. (c) Plot of the seismic network and array, as well as the previously assumed seismic and aseismic regions of the bed (Smith et al., 2015). (d) Plot of array in detail.

2 Data

The dataset used in this study is from receivers deployed on the surface of RIS, Antarctica. The instruments deployed consist of 16 4.5 Hz geophones connected to Reftek RT130 data loggers, sampling at 1000 Hz, in the geometry shown in Fig. 1. A total of 10 of these instruments were deployed in a 90 m aperture array, buried under several metres of snow. They were deployed in the austral summer of 2019 to 2020 (20 December 2019 to 23 January 2020). RIS flows at ∼400 m yr⁻¹, with the bed known to be seismically active from previous studies (Kufner et al., 2021; Smith et al., 2015; Smith, 2006; Hudson et al., 2021). RIS is therefore an ideal site at which to test array-based methods on glaciers.

3 Methods

3.1 Array processing

For the purposes of this study, we define a seismic array as a collection of seismic receivers configured in a geometry that allows the data to be processed using beamforming methods. In comparison, we define a network as a collection of seismic receivers that are used to perform conventional seismic detection and location, spaced so as to be most sensitive to events inside the network.

3.1.1 Detection and location

There are several different methods for detecting and locating icequakes using arrays. Here, we describe the specific method chosen for this study, a frequency-domain, frequency-wavenumber (fk) beamforming method, chosen because it is more computationally efficient than time-domain methods (Rost and Thomas, 2002). Similar methods are routinely deployed by various seismic observatories globally (Schweitzer et al., 2009).

The method comprises four overall steps (see Fig. 2):

1. Compute the beam power, slowness, and back-azimuth time series (P(t), S(t), Θ(t)).

2. Detect and associate icequake phase arrivals.

3. Locate each icequake (via either full-3D or fixed-depth methods).

4. Apply earthquake slowness and beam power filters in an attempt to remove falsely triggered events.

P(t), S(t), and Θ(t) are calculated using the frequency-wavenumber method as follows. First, the discrete frequencies on which to perform beamforming are specified. Here, we use 20 values between 10 and 150 Hz, representing the bandwidth within which icequakes have previously been observed at RIS (Smith et al., 2015; Kufner et al., 2021). The slowness limits for the beamforming calculation are then defined. The maximum slowness for this study is defined as 1 s km⁻¹, which is chosen based on an extreme estimate of the minimum apparent velocity (1 km s⁻¹) that one might expect to detect for an icequake at RIS. We also make the assumption about an incident plane wave at the array. This is valid for an approximately 90 m aperture array with a minimum hypocentral distance of 2.2 km (ice thickness). The velocity time-series data at each receiver are then windowed in time in order to calculate the beam power, slowness, and back azimuth for each time window. The time window used here is 0.2 s. This value is greater than the slowest time it would take for a seismic wave to travel across the array and provides frequency resolution up to 5 Hz. To provide sufficient phase arrival time resolution, the window time step is 0.01 s; hence consecutive windows overlap. The beamforming is performed as follows. The total beam energy E(u_i,θ) for a given slowness u_i and back azimuth θ_i across a given time window t_win is calculated as follows (Rost and Thomas, 2002; Bowden et al., 2021):

$$\begin{array}{}\text{(1)}& E({\mathit{u}}_i,{\mathit{\theta }}_i)={\displaystyle \frac{\mathrm{1}}{N}}\underset{-\mathrm{\infty }}{\overset{-\mathrm{\infty }}{\int }}\sum _{n=\mathrm{1}}^N{s}_n(t+\mathrm{\Delta }{t}_i{)}^{\mathrm{2}}\mathrm{d}t,\end{array}$$

where s_n(t)² is the power spectral density of the time series at a single receiver n, and Δt_i is the time shift required to align receiver n with the centre of the array as follows:

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }{t}_i,n({\mathit{u}}_i,{\mathit{\theta }}_i)=\left(\begin{array}{c}{u}_{x,i}\sin \left({\mathit{\theta }}_i\right)\\ u_{y,i}\cos \left({\mathit{\theta }}_i\right)\end{array}\right).{\mathit{r}}_n,\end{array}$$

for an array with negligible surface topography variation and where r_n is the distance of receiver n to the centre of the array. Here, we undertake beamforming in the frequency-wavenumber (fk) domain, where Eq. (1) can be rewritten as follows:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}E({\mathit{u}}_i,{\mathit{\theta }}_i)& ={\displaystyle \frac{\mathrm{1}}{N}}\underset{-\mathrm{\infty }}{\overset{-\mathrm{\infty }}{\int }}\left|\sum _{n=\mathrm{1}}^N{S}_n\right(\mathit{\omega })e^{\mathrm{2}\mathit{\pi }\mathit{\omega }\mathrm{\Delta }t}{|}^{\mathrm{2}}\mathrm{d}t\\ & ={\displaystyle \frac{\mathrm{1}}{N}}\underset{-\mathrm{\infty }}{\overset{-\mathrm{\infty }}{\int }}\left|\sum _{n=\mathrm{1}}^N{S}_n\right(\mathit{\omega })e^{\mathrm{2}\mathit{\pi }\mathit{\omega }{\mathit{u}}_i{\mathit{r}}_{\mathit{i}}}{|}^{\mathrm{2}}\mathrm{d}t,\end{array}\end{array}$$

where S_n(ω)² is the power spectral density of the time series in the Fourier domain, and k_i=ωu_i is the wavenumber associated with the given slowness and angular frequency. Here, we linearly stack over all discrete frequencies, so the integral in Eq. (3) simply becomes a linear sum. The power for a particular slowness and back azimuth over a particular time window $P({\mathit{u}}_i,{\mathit{\theta }}_i,t)$ is hence given by

$$\begin{array}{}\text{(4)}& P({\mathit{u}}_i,{\mathit{\theta }}_i)={\displaystyle \frac{E({\mathit{u}}_i,{\mathit{\theta }}_i)}{t_{\mathrm{win}}}}.\end{array}$$

One now has the beam power for all points in the 2D slowness space for each time window (see Fig. 2c). From this 2D slowness space, the peak beam power P(t) and its associated slowness S(t) and back azimuth Θ(t) can be calculated for each window. This process is performed for both the vertical and horizontal components of the seismic data separately.

Figure 2 (a, b) Plane wave arriving at array in the horizontal and vertical planes, respectively. All arrivals are time-shifted to align arrivals (diagram inspired by Bowden et al., 2021). (c) Heatmap of beam power in slowness back-azimuth space for one time window. An event is shown by the peak in beam power. (d) Schematic example of power time series, with vertical peaks associated with P arrivals and horizontal peaks with S arrivals. (e) Array-based 3D location algorithm, involving stacking PDFs of P–S travel time, P-wave takeoff angle, and S-wave takeoff angle. (f) Fixed-depth array-based location using back azimuth and P–S travel time only.

P(t), S(t), and Θ(t) can then be used to detect icequakes. For this study, the vertical component beam power time series P_V(t) is used to detect potential P-wave arrivals, and the two horizontal beam power time series are added together to give P_H(t), which is used to detect potential S-wave arrivals. Potential P- and S-wave arrivals are detected by searching for peaks in P_V(t) and P_H(t), respectively (see Fig. 2d). Here, we assume that a steep velocity gradient due to a near-surface firn layer causes approximately all P-wave energy to be incident on the vertical component and approximately all S-wave energy to be incident on the horizontal components. Surface waves may have energy on both vertical and horizontal components, but such phases are hopefully removed from the catalogue by the slowness-ratio filter. Peaks in P_V(t) and P_H(t) are identified using a median absolute deviation (MAD) threshold, with any peak in P(t) having a MAD multiplier value of 2 (optimized for this dataset) and a minimum event time separation >0.25 s defined as a potential phase arrival (based on the minimum P–S separation time for a basal icequake). Potential P and S phases then need to be associated. We only trigger an icequake detection if we can associate both P- and S-wave arrivals for an event. The phase association algorithm is as follows, with an example of the result shown in Fig. 3a. First, only potential P and S arrivals with a maximum separation of $\mathrm{\Delta }{t}_{\mathrm{P}-\mathrm{S},\max }=\mathrm{10}$ s are considered, to minimize the possibility of incorrectly associating arrivals. The value of Δt_P−S limits the hypocentral distance from which we can detect earthquakes (see Eq. 6). In order to account for possible overlapping event arrivals, only P and S arrivals with back azimuths in agreement within 15^∘ are paired. A third criterion is as follows: if multiple P and S arrivals meet the first two criteria, the highest beam power P arrival is associated with the highest beam power S arrival that has a back-azimuth difference of <15^∘ and a time difference of $<\mathrm{\Delta }{t}_{\mathrm{P}-\mathrm{S}}$, and so it continues for consecutively lower beam power arrivals. In practice, the maximum beam power criterion is rarely required, as two events are normally distinguishable by the $\mathrm{\Delta }{t}_{\mathrm{P}-\mathrm{S}}=\mathrm{10}$ s window using back azimuth alone.

The phase-associated icequake arrivals can then be used to locate events. Array-based icequake locations differ from network-based methods. Instead of using many P- and S-phase arrivals to locate an event, one has single P and S arrival observations. However, there is additional information on slowness and back azimuth associated with these two phase arrivals. To use an array to locate an icequake in 3D, one has to calculate a takeoff angle ϕ_EQ and a radial distance d_EQ from the event to the receiver, in addition to using the back azimuth. In this study, we compare two methods of finding icequake hypocentres: a 3D location method and a fixed-depth method (see Fig. 2e, f). The motivations for applying these two methods are discussed in Sect. 4.2. The takeoff angle for phase i (P or S) ϕ_EQ,i is calculated from the apparent velocity and slowness $v_{\mathrm{app},i}=\mathrm{1}/S_{\mathrm{app},i}$ and the approximate seismic velocity at the receiver v_rec,i, using the following equation:

$$\begin{array}{}\text{(5)}& {\mathit{\varphi }}_{\mathrm{EQ},i}={\sin }^{-\mathrm{1}}\left({\displaystyle \frac{v_{\mathrm{app},i}}{v_{\mathrm{rec},i}}}\right).\end{array}$$

In both cases, d_EQ is calculated using the P–S travel time delay Δt_P−S using the following equation:

$$\begin{array}{}\text{(6)}& d_{\mathrm{EQ}}={\displaystyle \frac{v_{\mathrm{P}}{v}_{\mathrm{S}}}{v_{\mathrm{P}}-v_{\mathrm{S}}}}\mathrm{\Delta }{t}_{\mathrm{P}-\mathrm{S}},\end{array}$$

where v_P and v_S are the path-average P- and S-wave velocities, respectively. These are approximated here to be the bulk ice velocities (v_P=3841 m s⁻¹ v_S=1970 m s⁻¹), obtained using a relationship between seismic velocity and temperature (Smith, 1997; Smith et al., 2015; Kohnen, 1974). In the 3D location method, for a vertically homogeneous ice structure, the ray path from the icequake source to the receiver is linear, so the takeoff angle, back azimuth, and hypocentral distance (d_EQ) can be used to directly back-project the event location in 3D for a vertically homogenous ice structure. However, at RIS there is a ∼100 m thick firn layer directly below the surface, which has increasingly slow seismic velocities towards the surface (Smith et al., 2020; Zhou et al., 2022) (see Sect. 3.1.2 for details). This causes the rays to dip steeply vertical as they near the surface. This effect has to be accounted for when locating icequakes, since the ray path cannot be assumed to be linear near the surface, with the takeoff angle representing the final ray path angle rather than the initial takeoff angle at the source. To account for a firn layer, we specify a vertically varying but laterally homogeneous velocity model and perform ray tracing on every grid point to a 2D depth vs. horizontal-radial-distance grid from the array (at 10 m grid resolution) to create a 2D spatial takeoff angle probability density function (PDF) (see Fig. 2e). Similarly, we create a second 2D spatial PDF for d_EQ, using Δt_P−S for a given icequake (see Eq. 6) and assuming an appropriate uncertainty δd_EQ again in depth and horizontal distance from the array for d_EQ. The 2D takeoff angle and d_EQ PDFs are then stacked and normalized to create an overall misfit space (see Fig. 2e), where the most likely icequake location within the 2D plane corresponds to the peak in this combined PDF space (see Fig. 2e). The 3D icequake location is then calculated by including the back-azimuth information. Conversely, in the fixed-depth method, the icequake location is calculated without using the takeoff angle, instead simply projecting the distance d_EQ onto a fixed-depth horizontal plane at the measured back azimuth (see Fig. 2f). The fixed-depth horizontal plane is equal to the ice thickness, explicitly assuming that all events originate at the bed.

Figure 3 Example of array phase detection and association performance. (a) Plot of beam power and peaks associated with phase arrivals for two events. Inset diagram shows possible focal mechanisms that could cause the observed P and S amplitudes. (d–f) Waveforms arriving at centre array station. (g–i) Mean of stacked and time-shifted waveforms for the entire array, using the slowness and back azimuth of event 2 in panels (b) and (c).

The initial icequake catalogue will likely contain both real events and false triggers. In order to discriminate between real events and false triggers, we apply two filters. The first filter is a slowness-ratio filter. In this work, we are interested in subsurface icequakes, where the array-processing correctly identified P and S body wave phases. P-wave, S-wave, and surface wave phases all have different propagation velocities and so will have different slownesses (since slowness is the inverse of velocity). However, the beamforming algorithm used here calculates apparent slowness rather than actual slowness, which is dependent upon the angle of incidence of the wave with respect to the surface. Apparent slowness therefore cannot be used to discriminate between phase arrivals. Instead, we use the slowness ratio ($S_{\mathrm{S}}/S_{\mathrm{P}}=v_{\mathrm{P}}/v_{\mathrm{S}}$), which is independent of apparent slowness effects if both the P- and S-waves have approximately the same ray path. A histogram of $v_{\mathrm{P}}/v_{\mathrm{S}}(=S_{\mathrm{S}}/S_{\mathrm{P}})$ ratio for the initial icequake catalogue is plotted in Fig. 4a. Ice at RIS has a $v_{\mathrm{P}}/v_{\mathrm{S}}$ ratio of 1.95 (Smith et al., 2015), approximately corresponding to the peak in Fig. 4a. We remove any events with $S_{\mathrm{S}}/S_{\mathrm{P}}<\mathrm{1.8}$ and $S_{\mathrm{S}}/S_{\mathrm{P}}>\mathrm{2.1}$. This filters our catalogue from ∼700 000 events to ∼250 000 events. The second filter we apply is a beam power filter, only keeping events with a combined beam power from both the P and S arrivals greater than 2×10⁸ counts² (1 count = $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{9}}$ m s⁻¹). This filter will likely remove some real events as well as false triggers, since smaller events that are further away from the array will inevitably have smaller beam power phase arrivals. Although the beam power filter is likely less sensitive to whether an event is real or a false trigger than the slowness filter, we still apply it in an attempt to remove more false triggers, so as to obtain a catalogue with hopefully a higher real event to false trigger detection rate. We use a beam power filter of 2×10⁸ counts² based on the assumption that englacial events directly beneath the array (with depths $<\sim \mathrm{2}$ km) at RIS are likely false triggers (see Fig. 4b). The beam power filter reduces our catalogue from ∼250 000 events to ∼25 000 events, which is herein referred to as the array-based icequake catalogue.

Figure 4 Summary of how the initial array-processing-derived icequake catalogue is filtered. (a) Histogram of slowness ratio for all earthquakes in the initial catalogue (note that $v_{\mathrm{P}}/v_{\mathrm{S}}=S_{\mathrm{S}}/S_{\mathrm{P}}$). Red dashed lines indicate filter values used to produce the final icequake catalogue. (b) Histogram of combined beam power (P_P+P_S) vs. depth below surface for icequakes with $\mathrm{1.8}<v_{\mathrm{P}}/v_{\mathrm{S}}<\mathrm{2.1}$, with depth derived using the 3D location method. Ice thickness is ∼2.2 km. The red dashed line indicates the minimum beam power filter used to obtain the final icequake catalogue.

3.1.2 Velocity model

At Rutford Ice Stream, there is a ∼100 m thick firn layer, overlying an assumed isotropic bulk ice column with a thickness of several kilometres. Bulk ice velocities are assumed to be v_P=3841 m s⁻¹ and v_S=1970 m s⁻¹ for the P- and S-waves, respectively. These are obtained using a relationship between seismic velocity and temperature, calibrated using seismic refraction methods (Smith, 1997; Smith et al., 2015; Kohnen, 1974). The 100 m firn-layer velocity structure used in this paper is based on the model derived from refraction seismic methods (Smith et al., 2020), which has been approximately independently verified using ambient seismic noise methods (Zhou et al., 2022). The full velocity model is included in a data repository (see ”Code and data availability” for details). The velocity model used for 3D location increases vertically from approximately v_P=660 m s⁻¹ and v_S=338 m s⁻¹ exponentially at the surface to v_P=3841 m s⁻¹ and v_S=1970 m s⁻¹ at 110 m below the surface and beyond, as in Smith et al. (2020). For the fixed-depth method, a homogeneous velocity model using the bulk ice velocities mentioned earlier is used. We assume that uncertainty in the velocity model is negligible compared to uncertainty in phase arrival times.

3.1.3 Uncertainty estimation

There are various sources of uncertainty in icequake phase arrival times and hypocentres. Sources include uncertainty in the velocity model, the GPS-derived receiver locations, and potentially anisotropic effects. Here, we estimate the overall uncertainty in event phase arrival times and location as follows. Phase arrival time uncertainties are defined as the full-width half maximum of the peak in the beam power time series associated with that particular phase arrival. Spatial uncertainties are defined by uncertainty in slowness and back azimuth. The uncertainty in slowness is defined as the full-width half maximum of the peak in beam power in the radial direction of the 2D slowness space (Fig. 2c). Similarly, the uncertainty in back azimuth is the full-width half maximum of the peak in beam power azimuthally in the 2D slowness space of Fig. 2c. In each case, the full-width half maximum is used, rather than the standard deviation of a Gaussian fit, to improve computational efficiency. We assume that there is negligible uncertainty in the velocity model. However, in reality, uncertainty in the velocity model will likely contribute to uncertainty in event locations that we have not accounted for in the array-processing-derived icequake locations.

3.1.4 Array sensitivity

Seismic arrays are sensitive to icequakes with a particular bandwidth, which is governed by the minimum and maximum spacing between individual receivers. The minimum optimal frequency that an array is sensitive to is proportional to the maximum receiver spacing and vice versa for the maximum optimal frequency. The equation that describes this is given by (Bowden et al., 2021)

$$\begin{array}{}\text{(7)}& f={\displaystyle \frac{v}{x_{\mathrm{rec}}}},\end{array}$$

where f is frequency, v is the seismic velocity, and x_rec is the receiver spacing. For the array in this study, the minimum and maximum receiver spacings are 20 and 90 m, respectively. The corresponding optimal bandwidth for the array is 40 to 200 Hz for P-waves and 20 to 100 Hz for S-waves, using the bulk ice velocity for P- and S-waves, respectively. Depending upon the array geometry and level of radial symmetry, the receiver spacing and therefore sensitivity could vary with azimuth. The array in this study is designed to be approximately radially and azimuthally symmetric, so the azimuthal sensitivity is approximately constant.

3.1.5 A note on stacking

Stacking the phase-shifted waveforms within the beamforming process suppresses incoherent noise. In this study, we linearly stack the data. However, one could also use nth root stacking (Rost and Thomas, 2002) or phase-weighted stacking (Schimmel and Paulssen, 1997). One could also consider stacking different frequency ranges separately in order to further improve sensitivity (Gal et al., 2014). We find negligible differences in performance, using different stacking techniques for the array setup in this study. This is likely because the array has a small (100 m) aperture, and minimal scattering occurs within the shallow firn layer at RIS, therefore resulting in insignificant incoherent noise compared to local icequake signals, regional earthquakes, and coherent ambient seismic noise sources. We would recommend reconsidering the choice of stacking method for larger-aperture arrays or sites with significant near-surface heterogeneity, either of which would increase incoherent noise levels.

3.2 Network-based icequake detection and location

The array-based detection and location method presented in this study is benchmarked against a current network-based migration method for icequake detection, QuakeMigrate (Hudson et al., 2019). This method approximates the energy from an icequake arriving at each receiver with a Gaussian onset function, which is then back-migrated through time and space to search for a coalescence of energy that corresponds to an event. If a sufficiently high coalescence of energy is found over a particular time window, it triggers an event detection. Full details on this method can be found in Hudson et al. (2019). QuakeMigrate provides an estimate of a hypocentral location for each event that equates with the 3D grid node, corresponding to the maximum coalescence of energy. We therefore relocate all events using the non-linear, probabilistic earthquake location algorithm NonLinLoc (Lomax and Virieux, 2000) in order to obtain a precise location and physically meaningful estimates of hypocentral uncertainty. The QuakeMigrate method shares similarities with the beamforming method (see Sect. 3.1), with both seeking to identify coherent arrivals of energy within a particular frequency bandwidth. The fundamental difference between the two methods is that for plane-wave array processing, earthquake sources need to be sufficiently far away that the plane-wave assumption holds, whereas the network-based method is optimized for sources within the network extent. Another key difference between the two is that the array-based method searches a pre-defined slowness azimuth space, which in this case is undertaken in the fk domain, while the network-based approach searches over a pre-defined 3D spatial grid through time.

3.3 Moment magnitude

One can use the icequake magnitude distribution to quantify the sensitivity of the array compared to the network for smaller vs. larger events. We choose to use moment magnitude M_w (Hanks and Kanamori, 1979), as it provides an absolute measure of the actual moment release of an icequake, rather than relying on an empirical relationship. Another benefit is that unlike local magnitude scales, it does not exhibit a break in the scaling relationship at low magnitudes (M_w<3) (Deichmann, 2017; Hudson et al., 2022). Moment magnitudes are calculated using SeisSrcMoment (Hudson, 2020) (see Hudson et al., 2022, for details on applying the method). This involves fitting a Brune model (Brune, 1970) to the frequency spectrum of the icequake.

4 Results and discussion

4.1 Icequake detection using arrays vs. networks

Detected icequake hypocentres are shown in Fig. 5, with a summary of the magnitude distribution of these icequakes shown in Fig. 6. The magnitude distributions for three detection setups are shown in Fig. 6: (1) the array-based detection method applied to the array shown in Fig. 1d (Fig. 6b), (2) the network-based detection method applied to the same array shown in Fig. 1d (green points, Fig. 6a), and (3) the network-based detection method applied to the entire network shown in Fig. 1c (red points, Fig. 6a). The results allow for the comparison of an array deployment with a network deployment, as well as array-based with network-based detection and location algorithms more generally.

The array-based detection method outperforms the network-based detection method in several areas. First, the array detects more icequakes than the network (see Fig. 6). The additional icequake detections in the array data originate predominantly at greater distances from the array/network centre than icequakes detected using the network-based approach (see Fig. 5). Icequakes continue to be detected at distances of 10 km or more, albeit with an increasing average magnitude with distance, whereas the network has a sharp detection limit at ∼7.5 km from the network centre (see Fig. 6d). The sharp limit for the network-derived catalogue is the consequence of the boundary of the network-based search grid, with the search area significantly exceeding the physical microseismic detection limits of the network-based algorithm. Additionally, the array detects ∼1000 icequakes with magnitudes >0, whereas the network only detects a negligible number of these larger icequakes, many of which originate far outside the spatial extent of the network (see Fig. 6). The array-based method also detects a higher proportion of the smaller icequakes, with a magnitude of completeness $M_{\mathrm{c},\mathrm{array}}=-\mathrm{1.4}$ compared to $M_{\mathrm{c},\mathrm{network}}=-\mathrm{0.81}$ for the network-based method.

However, the network-based detection method outperforms the array-based detection method for discriminating false triggers. This is evidenced by the clustering of seismicity at sticky spots clearly exhibited in the network-based event locations (Fig. 5c). The array-based data exhibit some clustering near the array (Fig. 5a), but as hypocentral distances from the array increase, the events become more scattered, likely a combination of both false triggers and poorer hypocentral constraint.

Overall, the array-based method is more sensitive than the network-based method, detecting more icequakes across the magnitude range. This result is likely received for two reasons. First, the array-based method outperforms the network-based approach in phase association, with multiple P- and S-wave phase arrivals possible to associate within a given time window. This is possible due to the accuracy of back-azimuth measurements, allowing arrivals with back azimuths differing by <15^∘ to be paired, even if phase arrivals overlap. This is particularly powerful when accounting for radiation pattern effects, as shown in Fig. 3. For the example in Fig. 3, two events close in time originating from different back azimuths have inverse P/S amplitude ratios, due to radiation pattern effects, yet both are detected within the same window. Theoretically, tens of events could overlap within each time window, limited only by the back-azimuth tolerance and distribution of event back azimuths. This is in contrast to the network-based approach, where only one event association would be allowed within a given time window, so as to minimize the risk of incorrect phase associations. However, the greater number of events detected by the array-based method could also be due to a lack of metrics to filter the catalogue by. Although our array-processing method provides uncertainty estimates, these have a particularly coarse temporal resolution, limiting their use for filtering the data to remove false event detections. Conversely, the network-based approach measures uncertainty with a higher resolution, allowing both temporal and spatial uncertainty filters to be used to remove false detections. However, given our strict (<15^∘) back-azimuth phase association criterion and slowness-ratio filter, we are confident that the difference between the array-based and network-based methods cannot be attributed solely to false event triggering.

Figure 5 Comparison of array-derived to network-derived icequake catalogue. (a, b) Array-derived icequake catalogue hypocentres. Red points use the full 3D array location method. Grey points use the fixed-depth method (transparency corresponds to icequake amplitude). A^′ and B^′ denote the extent of the plotting in panel (c) and (e), respectively. (c, d) Network-derived icequake catalogue hypocentres, coloured by icequake cluster (see text for details). (e, f) Enlarged plot showing only the array-derived icequake hypocentres obtained using the full 3D array location method. As in panels (a) and (b), transparency corresponds to icequake amplitude.

Figure 6 (a) Moment magnitude distribution for entire network and array stations only, with events detected using migration-based method. Lower scatterpoints indicate data for each individual bin, while upper scatterpoints indicate cumulative moment magnitude distribution. (b) Same as panel (a) but for arrays using the array-based detection method (using the fixed-depth hypocentres; see Fig. 5) and the mean of the linearly stacked waveform data from the array. (c) Plot of histogram-binned data for individual spatial icequake clusters from the entire network data from panel (a). Data are coloured to corresponding clusters in Fig. 5c. (d) Plot of mean M_w with distance from the array/network centre for the entire network data from panel (a) and the array data from panel (b).

A note on icequake magnitude distributions

Tectonic earthquake magnitudes typically follow a logarithmic scaling relationship (Gutenberg and Richter, 1936, 1944). The array-based icequake detection results shown in Fig. 6b exhibit a similar relationship, with the tail-off at magnitudes below the magnitude of completeness M_c caused by icequake signal-to-noise ratios (SNRs) falling below the noise level, leading to not all events being detected. However, the network-based icequake magnitude distributions do not exhibit a clear linear trend (see Fig. 6a). This effect is not caused by S-wave anisotropy or assumptions about the source mechanism orientation, with P-wave M_w and average moment tensor M_w distributions exhibiting similar peaks and troughs in the binned data. We instead find that these peaks and troughs in the M_w distribution are caused by spatially distinct icequake clusters with their own narrow magnitude distributions (see Fig. 6c). The limited extent of magnitude variation for each cluster is presumably governed by bed properties, whether that be the extent of slip, the rupture velocity of the ice–bed interface, or other similar effects (Gräff and Walter, 2021; Hudson et al., 2023; Zoet et al., 2012). This clustered distribution of icequakes also likely plays a role in mean M_w with distance (see Fig. 6d), where the network is more sensitive to small icequakes in clusters within the network extent. The array-based detection results do not exhibit this cluster-dominated behaviour since it can detect events at greater distances, therefore sampling a greater distribution of icequake clusters and potentially icequake sources.

4.2 Icequake location using arrays vs. networks

The network outperforms the array in icequake location. This is evident from the clearly discernible icequake clusters in the network-based icequake catalogue shown in Fig. 5c, expected at RIS based on findings from the same area of bed when a much denser seismic network was deployed (Kufner et al., 2021). Similar observations are obtained when locating the icequakes only using the 10 inner array stations with the network-based location method. Icequake clustering is generally indiscernible in the array-based icequake location results in Fig. 5a. One reason for this is the filtering of high-quality icequake phase arrivals in the network data using temporal and spatial uncertainty measurements, as described in Hudson et al. (2019). Noisy, poorly constrained yet real icequakes are likely filtered out by definition in the network-based results, whereas these noisier, low-SNR icequakes are more likely kept in the array-based catalogue compared to the network-based catalogue. This could also affect the behaviour of M_w with distance (see Fig. 6d).

However, there is also a more fundamental limitation on the location results of the array-based method: the presence of a near-surface, low-velocity firn layer (Smith et al., 2020; Zhou et al., 2022) that causes seismic waves to steeply dip towards the surface (see Fig. 2). Glacier settings with a thinner or no firn layer would result in better constraint on event location. Uncertainty in the velocity structure of the firn layer, especially at P-wave wavelengths (<10 m), limits the measurement of the takeoff angle from apparent slowness used in the array-based method's 3D icequake location procedure. This is what causes the icequakes located using the array-based method to be incorrectly located, approximately directly beneath the array (red scatterpoints, Fig. 5e, f). The firn-layer effect on the 3D location method is most pronounced in Fig. 5f, where events are projected below the ice–bed interface. To mitigate this issue, we are forced to neglect firn-layer effects and project icequake epicentres onto an artificial horizontal plane at approximately the depth of the ice–bed interface using the fixed-depth location method (see Fig. 2f). Obviously, this is an approximation, with both neglecting the firn layer and differences between the average ice–bed interface and the true icequake depth resulting in greater uncertainty in the icequake epicentres, likely making any icequake clusters challenging to discern. High-resolution (<10 m) 3D imaging (lateral in addition to vertical) of the firn velocity structure beneath the array could allow one to apply an array transfer function to minimize these effects, although without access to such data, we cannot test this hypothesis. We therefore include the 3D-method location results to emphasize the importance of understanding the near-surface velocity structure when performing array processing.

4.3 New insights into Rutford Ice Stream

The array deployment provides new insight into seismicity at RIS. Previous studies of RIS suggest that bed properties vary upstream of the deployment vs. downstream (Smith, 1997; Smith and Murray, 2009). The upstream bed is thought to be comprised of unconsolidated sediment that fails aseismically, while the bed downstream can fail seismically at sticky spots (Smith et al., 2015) (see Fig. 1). Kufner et al. (2021) use a larger network of 35 receivers to find icequakes upstream of the seismic–aseismic boundary (blue points and purple line, respectively, Fig. 7a). These are interpreted to occur in the depressions of mega-scale glacial lineations (MSGLs) in the previously aseismic region. Array-processing enables event detection at greater hypocentral distances, allowing us to confirm the findings of Kufner et al. (2021), with seismicity extending further upstream again, likely only limited by array sensitivity, the maximum $\mathrm{\Delta }{t}_{\mathrm{P}-\mathrm{S},\max }$ time that we impose on P–S phase association, and our false-trigger beam power filter. Figure 7b shows an example of a previously undetected upstream event. The icequake has high-SNR P- and S-wave arrivals. The network-based method would likely fail to detect the event, even if the search grid were sufficiently large, since there are multiple phase arrivals in the window from other overlapping events. Overall, observing seismicity in the previously inferred aseismic region has implications for bed complexity and icequake nucleation, potentially supporting ideas such as basal water pressures modulating bed friction and seismicity (Hudson et al., 2023; Gräff et al., 2021).

Figure 7 Examples of previously unobserved icequakes at Rutford Ice Stream. (a) Map of overall seismicity detected in this study within the context of RIS more widely. Grey points represent the icequake catalogue from array processing in this study, blue points are icequakes from Kufner et al. (2021), the yellow dashed lines indicate the shear margins, and the purple solid–dashed line indicates the boundary between the previously inferred seismic–aseismic region (Smith, 1997; Smith and Murray, 2009). (b) Example icequake located upstream in the “aseismic” region from the $\mathrm{\Delta }{t}_{\mathrm{P}-\mathrm{S},\max }=\mathrm{2.5}$ s icequake catalogue. (c) Downstream icequake from a previously unstudied region. (d) Shear-margin icequake, again from a previously unresolved region.

The array-derived icequake catalogue also contains earthquakes further downstream than those found in Kufner et al. (2021) (see Fig. 7a). It is expected that events might originate from this region. However, what is surprising is the number of seismic signals in the time series (see Fig. 7c). Typical icequake repeat times for individual spatial clusters are of the order of hundreds to thousands of seconds at RIS (Hudson et al., 2023), and the waveforms do not look similar, so the other potential events within this window likely originate from various locations. The high number of signals within the 6 s window shown in Fig. 7c would likely be challenging to separate using the network-based detection algorithm, but the array-based algorithm can associate P and S phases based on the slowness measurements. This emphasizes the phase-association benefits provided by the array-based method compared to the network-based detection method.

The final observation in the array results that we emphasize here is the occurrence of icequakes at the shear margins of RIS. Again, these events are observed due to the greater detection distances of array-based methods (up to 20 km) compared to networks (<7.5 km). An example icequake from the shear margin is shown in Fig. 7d. This event also likely originates at or near the glacier bed, since the slowness-ratio false-trigger filter applied to the icequake catalogue should remove any surface wave detections. The large-amplitude S-wave compared to the P-wave is likely a combination of the position on the icequake focal sphere and perhaps also the higher shear rates near the ice stream shear margins. The S-wave also appears to potentially exhibit shear-wave splitting associated with seismic velocity anisotropy. Detecting such icequakes could provide information on shear-margin dynamics, such as how important damage to the ice fabric is for impeding or accelerating ice stream flow.

4.4 Lessons learnt and recommendations for using arrays in the cryosphere

1. Seismic networks are more sensitive than arrays for studying smaller areas in more detail. The network-based approach has a smaller magnitude of completeness than the array-based approach while also providing greater spatial constraint on icequake hypocentres that allow clusters of events to be identified. For studying a particular glaciological process in as much detail and as comprehensively as possible, we recommend deploying a seismic network rather than an array. This is especially true for sites with a thick (greater than the seismic wavelength) low-velocity firn layer.

2. Arrays can detect icequakes at greater distances than networks. Arrays can outperform networks in detecting events, especially when P- and S-wave phase arrivals overlap in time. Theoretically, the array-processing method here can detect many events within a given time window, compared to a single event using the network method. Furthermore, phases from an event that has highly differing SNRs can still be readily associated as they are from the same event. In this work, we implement a standard fk-domain method, stacking all frequency ranges equally, but the sensitivity of our method for simultaneous event detection could be further improved by stacking different frequency ranges separately (Gal et al., 2014). Together, these properties of the array-processing method enable more events to be detected from greater distances and during noisier time periods. We therefore suggest that arrays are useful for the initial scoping of a field site before potentially deploying a more comprehensive network to study a particular process.

3. Arrays typically require fewer instruments than a network. The array in this study is comprised of only 10 receivers, whereas a seismic network would ideally comprise tens to hundreds or more receivers. Arrays therefore provide an efficient means to investigate seismic activity, at least initially.

4. Array and network geometries limit detection performance in different ways. As described earlier, arrays and networks have different advantages and compromises. Arrays require receiver spacings optimized to the spectral content of the earthquakes to be detected, which may not be known in advance. In this study, that spacing is 20 to 90 m. In contrast, networks generally perform best when the receivers are evenly spaced, when icequake depths are not greater than the maximum receiver spacing. The optimal receiver spacing for a network to study basal icequakes at RIS (∼2 km) is therefore much greater than the maximum optimal array receiver spacing (∼100 m). To summarize, it is challenging to design a deployment optimized for both array and network processing.

5. Consider networks of sub-arrays to capitalize on the advantages of both networks and arrays. If one has a sufficient number of receivers, then it may be possible to deploy evenly spaced sub-arrays within an overall network. This could facilitate a hybrid detection approach, taking advantage of both network and array benefits.

6. Consider generating a more complete catalogue of events based on the initial catalogue using other methods. Once a catalogue of initial events has been obtained, one could use other methods such as template matching to increase the number of events detected (Helmstetter, 2022; Gimbert et al., 2021; Gibbons and Ringdal, 2006).

5 Conclusions

Here, we focus on how useful array processing is for deployments used to study icequakes at glaciers. The motivation for using arrays rather than networks is that the cryosphere is often challenging to access, logistically expensive to operate within, and potential seismically active areas are typically highly variable temporally and spatially. Seismic arrays could facilitate smaller seismic deployments than conventional networks while enabling event detection at greater distances. We find that arrays can detect icequakes over a greater spatial extent than seismic networks but provide poorer spatial constraint on seismicity within a network, where networks have the potential to elucidate glaciological processes in greater detail. Array slowness ratios play an important role for discrimination of real events from false triggers at greater distances. At Rutford Ice Stream, events are detected in what was previously thought to be an aseismic region with different bed properties, downstream where icequakes have never previously been observed, and at the otherwise inaccessible ice stream shear margins. These results emphasize the value of array-based icequake detection for more comprehensively studying glacier dynamics over larger spatial footprints than otherwise possible. We make a number of recommendations based on the findings from this study, especially that arrays might be particularly useful for initial scoping field deployments, where one only has access to an insufficient number of instruments to deploy a suitable network, or for regions that are potentially hazardous to operate in, such as crevassed shear margins.