Not all Icequakes are Created Equal: Basal Icequakes Suggest Diverse Bed Deformation Mechanisms at Rutford Ice Stream, West Antarctica

Microseismicity, induced by the sliding of a glacier over its bed, can be used to characterize frictional properties of the ice‐bed interface, which are a key parameter controlling ice stream flow. We use naturally occurring seismicity to monitor spatiotemporally varying bed properties at Rutford Ice Stream, West Antarctica. We locate 230,000 micro‐earthquakes with local magnitudes from −2.0 to −0.3 using 90 days of recordings from a 35‐station seismic network located ∼40 km upstream of the grounding line. Events exclusively occur near the ice‐bed interface and indicate predominantly flow‐parallel stick‐slip. They mostly lie within a region of interpreted stiff till and along the likely stiffer part of mega‐scale glacial lineations. Within these regions, micro‐earthquakes occur in spatially (<100 m radius) and temporally (mostly 1–5 days activity) restricted event‐clusters (up to 4,000 events), which exhibit an increase, followed by a decrease, in event magnitude with time. This may indicate event triggering once activity is initiated. Although ocean tides modulate the surface ice flow velocity, we observe little periodic variation in overall event frequency over time and conclude that water content, bed topography and stiffness are the major factors controlling microseismicity. Based on variable rupture mechanisms and spatiotemporal characteristics, we suggest the event‐clusters relate to three end‐member types of bed deformation: (1) continuous creation and seismogenic destruction of small‐scale bed‐roughness, (2) ploughed clasts, and (3) flow‐oblique deformation during landform formation or along bedrock outcrops. This indicates that multiple processes, simultaneously active during glacial sliding, can accommodate stick‐slip behavior and that the bed continuously reorganizes.

Insights into basal conditions can be gained through the study of icequakes, if their hypocenters lie near the ice-bed interface (Röösli et al., 2016;Smith, Smith, et al., 2015). Basal icequakes have been detected widely at glaciers in Antarctica and elsewhere (Anandakrishnan & Bentley, 1993;Blankenship et al., 1987;Danesi et al., 2007;Helmstetter, Nicolas, et al., 2015;Röösli et al., 2016;Smith, Smith, et al., 2015;Walter, Deichmann, et al., 2008) and several reasons have been suggested for their occurrence. These include localized bed heterogeneities, water-pressure fluctuations or water-induced crack opening (Fischer, Clarke, et al., 1999;Smith, Smith, et al., 2015;Walter, Canassy, et al., 2013). There is also the possibility that a combination of these mechanisms may be at play simultaneously. However, many icequake studies suffer from short deployment times or heterogeneous network geometry. This makes it difficult to isolate the effect of spatiotemporally varying basal properties on icequake occurrence. Here, we use 90 days of passive seismic data to detect basal microseismicity at Rutford Ice Stream, West Antarctica ( Figure 1). The data were recorded by 35 seismometers with a nominal spacing of 1 km over a 10 × 10 km grid deployed ∼40 km upstream of the grounding line ( Figure 1b). Within this area, seismic surveys have shown that the bed consists of till, with varying water content, consolidation state, and degree of deformation (King, Pritchard, et al., 2016;Smith, 1997;Smith & Murray, 2009). Furthermore, it has been shown that surface ice flow is heavily modulated by a biweekly tidal signal (Gudmundsson, 2006;Minchew et al., 2017). Thus, our seismic network covers a region of diverse bed topography and rheology and captures several tidal cycles. This allows us to investigate basal slip in an ice stream with unprecedented spatiotemporal resolution. Here we locate icequakes, but also determine their source characteristics including event magnitude, source mechanisms and spatiotemporal clustering. Based on these results, we can better constrain the mechanisms for seismicity and how the icequakes reveal the basal properties of the Rutford Ice Stream (RIS). Figure 1a) drains ∼49,000 km 2 of the West Antarctic Ice Sheet into the Ronne Ice Shelf (Doake et al., 2013). To the west and east, RIS is bound by the Ellsworth Mountains and the Fletcher Promontory, respectively. At our study site, the ice flow velocity is ∼375 m a −1 (Adalgeirsdóttir et al., 2008) and the ice stream is around 2.2 km thick and grounded at 1.6-1.8 km below sea level (King, Hindmarsh, et al., 2009). RIS occupies a deep trough with a "w-shaped" cross-section (King, Pritchard, et al., 2016). The center of our network is deployed on the ice stream surface above a basal central high (∼1.8 km below sea level); the bed topography descends to the SW and NE into troughs on either site (Figure 1b). Slightly downstream (∼2 km) of our survey location, the bed topography is dominated by a prominent knoll, which also creates a surface expression. By contrast, the ice surface is flat within our seismic network. Nevertheless, the bed below the seismic network features a diverse morphology. This morphology can be emphasized through the "residual elevation" (Figure 1c), calculated by King, Pritchard, et al. (2016) by subtracting a filtered version of the bed DEM (2 × 2 km smoothing filter) from the original data. In the upstream part of our seismic network, the short-wavelength topography is formed of several elongated megascale glacial lineations (MSGLs), of which the central one is the most prominent. The topography in the downstream part of our network is more irregular and features multiple hummocks of nonuniform shape, orientation, and size (King, Pritchard, et al., 2016;Figure 1c).

RIS (
Based on radar and seismic surveys, it has been shown that the bed is composed of a basal till layer with spatially varying properties, resulting in different basal deformation regimes (King, Hindmarsh, et al., 2009;Schlegel et al., 2021;Smith, 1997;Smith & Murray, 2009). Upstream, the MSGLs are likely composed of water-saturated, deformable till. Seismic surveys over these landforms, repeated over timescales of a few years, reveal sediment transport and bedform erosion of up to 1 m a −1 at the downstream termination of the MSGLs (Smith & Murray, 2009;Smith, Murray, et al., 2007). The deformable till layer likely overlays a stiffer and more consolidated unit that outcrops locally and predominantly northwest of the central high (Schlegel et al., 2021;Smith & Murray, 2009). Downstream of the MSGLs, the till layer is generally stiffer and likely stiffest southeast of the central high where very consolidated till or possibly sedimentary rock are proposed to exist (Schlegel et al., 2021). This first-order discrimination of different bed domains was supported by drilling results ; see Figure 1c for drill locations). In the regions of stiffer till, basal sliding rather than bed deformation predominantly accommodates basal motion at the ice-bed interface. The two domains of dominantly bed deformation and basal sliding can be discriminated from each other based on their characteristic polarity and intensity reflection values from seismic surveys (Smith, 1997). Together with the different geomorphological appearance (MSGLs vs. hummocks), this allows for the definition of a "bed character boundary," separating the two domains (King, Pritchard, et al., 2016;G. Boulton-pers. communication in Smith, Smith, et al., 2015). This boundary is highlighted in Figures 1b and 1c and all subsequent map view figures. In the following, we will refer to the domain upstream of this boundary as a "soft sediment" bed and to the domain downstream of the boundary as a "stiff sediment" bed. This is based on the assumption that the bed character variability across the boundary mainly depends on porosity because the bed composition may be similar.
KUFNER ET AL.   Figure S1. (c) Zoom into the study region. Background color coding demarcates residual bed topography, which is calculated based on the difference between the short-wavelength topography and a long-wavelength trend surface (King, Pritchard, et al., 2016). Hummock locations and dashed bed character boundary are from King, Pritchard, et al. (2016), while the dotted pink-purple line represents an alternative bed character boundary defined by G. Boulton (pers. communication in Smith, Smith et al., 2015). Gray circles indicate the location of hot-water drill sites that where operated during the BEAMISH 2018/19 season (A. M. Smith, Anker, et al., 2020). Our study site has been the focus of previous passive seismic surveys. In a pioneer study, Smith (2006) deployed 10 geophones in a circular array with 3 km interstation spacing for a 11-day observation period. The recordings from this network were not sensitive enough to allow for precise event locations but showed six times more basal micro-earthquakes originated from regions of stiff sediments. Smith, Smith, et al. (2015) improved this understanding using another 10-station network, deployed in two subarrays with 1 km station spacing over a period of 35 days. Due to higher sensitivity instruments and different network configuration, 3,000 basal icequakes were precisely located, mostly in areas of stiff sediment bed. These confirmed findings from the earlier study and suggested that basal ice flow mechanics depend on basal conditions. Furthermore, they showed that seismicity generally featured low-angle faulting mechanisms, which indicates basal sliding in the flow direction as major source triggering seismicity.
In addition to icequake occurrence, the basal hydraulic system of RIS varies dependent on bed rheology . Based on radar and seismic surveys, it was shown that water channels or bodies exist in the region of soft sediments over a long distance along the landforms (at least 1 km long and 200 m wide, King, Woodward, et al., 2004;Murray et al., 2008;and Schlegel et al., 2021). Furthermore, water may be present on top of MSGL ridges (King, Hindmarsh, et al., 2009;Murray et al., 2008;Schlegel et al., 2021). Within the stiffer sediment region, free water appears in isolated spots and pools Schlegel et al., 2021).
Lastly, an ice stream flow velocity modulation, related to the spring-neap tidal cycle, has been measured at RIS (Adalgeirsdóttir et al., 2008;Gudmundsson, 2006;Murray et al., 2007). At the grounding line, the biweekly modulation of the surface ice stream flow velocity is up to 20%. This signal propagates, with decreasing amplitude, up to 60 km upstream. A linked hydrological and numerical modeling study suggests that only a combination of stress transmission through the ice and changes of basal water pressure can explain such modulations in surface ice flow velocity (Rosier et al., 2015). In addition, the model assumes a highly effective basal drainage system, low effective pressure and a nonlinear sliding law.

Network Description
We use 3 months (mid-November 2018-mid-February 2019) of continuous passive seismic recordings to generate a microseismic event catalog. This data set was collected as part of the BEAMISH project (A. M. Smith, Anker, et al., 2020) during the 2018/19 field season. The seismic network broadly forms a rectangle with ∼1 km station spacing. It overlays both bed domains ( Figure 1c). The geometry of the network was modified twice during the observation period. Initially, 19 stations in the northern and central part of the network were deployed in November 2018. The network was then extended with 14 stations to the east and north in December. In January 2019, two further stations were added and three stations from the westernmost corner of the array were redeployed in the central part of the network. In total, 38 sites were occupied, with 19-35 stations recording concurrently. During a strong storm in December, 15 of these stations were inactive for up to 5 days (see Figure S1 for details).
Each station consisted of a Reftek RT-130 data logger with a 4.5 Hz three-component geophone (either GS11-3D or L28-3D), which was buried to ∼1 m depth (see Figure S1 for details of the network). The sampling frequency was 1,000 Hz. Energy supply was ensured through a solar panel and battery. Timing was obtained from an attached GPS antenna.

Microseismic Event Catalog and Spatial Clustering
Thousands of microseismic events were recorded during the deployment period. These events tend to cluster closely spaced in time, are characterized by an impulsive P-wave onset, and two prominent S-wave arrivals ( Figure 2a). The detection of two independent shear waves is an indication of the anisotropic nature of the ice comprising the RIS (Smith, Baird, et al., 2017). Typical frequencies for P-waves are between 10 and 200 Hz. S-wave frequency is predominantly between 30 and 100 Hz. We use the QuakeMigrate software J. D. Smith, White, et al., 2020) to detect and locate events from the continuous seismic records. Instead of a classic station-by-station trigger, QuakeMigrate implements a detection scheme based on the coherency of seismic phase arrivals recorded at all seismic stations. This makes it an ideal detection tool if many, temporally overlapping, small earthquakes occur. Based on the P-and S-wave onset times and uncertainties derived in QuakeMigrate, the initial locations are refined using NonLinLoc (Lomax et al., 2000), which yields a more realistic location error estimate due to a probabilistic location approach and a weighting scheme for pick uncertainties. For QuakeMigrate, we use a homogeneous velocity model (vp = 3.841 km s −1 ; vs = 1.970 km s −1 ). The P-wave velocity (vp) corresponds to the ice velocity at RIS obtained from an earlier seismic survey (Smith, 1997). The S-wave velocity (vs) is derived from vp and the vp/vs ratio of 1.95 taken from a Wadati diagram using 1 day of data (∼36,000 P-and S-picks; Figure S2). In NonLinLoc, we further refine the velocity model and included a uniformly 100 m thick layer to represent firn (vp = 2.839 km s −1 ; vs = 1.456 km s −1 ; Smith, 1997) below the seismic stations and above the solid ice. Uncertainties in the velocity model (according to Smith, 1997, less than ±0.015 km s −1 ) are assessed through relocating sample events, while considering travel time dependent errors (see Section S1 for details). We do not include a velocity discontinuity below the ice, which would represent the glacial bed, as such a layer is likely to introduce artificial event clustering (Smith, Smith et al., 2015). Theoretical calculations based on till and bed properties expected at RIS (Smith, 1997) showed that the direct upgoing P-phase likely forms the first arrival for epicentral distances of up to 10 km. Less than 0.01‰ of our P-and S-picks have greater epicentral distances. Thus, including only the ice and firn layer in the velocity model does not lead to misidentified phases. Furthermore, we do not include the effect of anisotropy in ray tracing, but tune QuakeMigrate through the "detection threshold" parameter to pick the first possible S-wave onset. We assess the effect of anisotropy by using a sample event, which is located using the first and second peak in the S-waveform, respectively ( Figures S3 and S4). The discrepancy between the two resultant hypocenters is most significant in the vertical direction and can be neglected in horizontal directions, which we consider when discussing the results.
We apply quality restrictions to the automatically detected picks and events to ensure that no false picks and events are included in the final event catalog. We accept only events with a total root-mean-square (RMS) value of travel time residual of 0.02 s at the maximum likelihood hypocentral location, a maximum azimuthal gap of 280°, maximum 10% of picks with a P/S travel time residual (observed subtracted by predicted arrival time) larger than 0.02/0.2 s, and at least three P-picks and two S-picks. These selection criteria are obtained from the visual inspection of data sub-sets and reduce 295,785 potential events initially detected KUFNER ET AL.   Figure S1). Amplitude is in instrument counts. The windows used for M w derivation and the maximum amplitude used to calculate ML are highlighted. (b) Magnitude histogram for all events in 0.05 bins. "Cumulative values" refer to the all events greater or equal to a specific magnitude according to the Gutenberg-Richter law (Gutenberg & Richter, 1944). Solid lines represent regression lines based on the cumulative values. Sections with different log(ML) decay slopes ("b-values") are highlighted in red and blue, respectively. from QuakeMigrate to 227,029 events. Further details on location methodology and implementation to our data set can be found in the supplementary material (S1 and S2 and Table S1).
We account for the movement of the seismic stations relative to the bed due to ice flow by shifting each event location in the final catalog downstream. RIS moved ∼94 m downstream during the 90-day survey period, whereas our stations are specified at fixed locations during event location, clearly evidencing the necessity for such a shift. We perform this shift by calculating the stations' locations at the time of each individual event relative to the start of the deployment, using their GPS locations, and apply this lateral shift to that event hypocenter. This is repeated for all events in the catalog to compensate for ice flow.
Finally, we group the events into clusters as glacial microseismicity is known to occur in bursts of temporal and spatially focused activity (Smith, Smith, et al. 2015). We apply the DBSCAN ("Density-Based Spatial Clustering of Applications with Noise") cluster algorithms (Pedregosa et al., 2011) to search for spatial patterns in our microseismic event catalog. This SciKit python module is designed to find core samples of high spatial density and to extend clusters around them. Only events with magnitudes larger than the magnitude of completeness (including a 0.2 magnitude units buffer to account for uncertainty), which we determine as the maximum in the logarithmic magnitude plot of Figure 2b, are included in the cluster analysis. This event cut-off is implemented to avoid a bias in the output clusters due to spatially differing completeness magnitudes (see details on parameterization in Section S3).

Event Magnitudes
Magnitudes are calculated using a two-step approach. First, we determine the moment magnitude (M w ; Hanks & Kanamori, 1979) for a subset of our data (January 1, 2019, 1,520 events) from the far-field displacement of the P-wave (Shearer, 2009;implementation of Hudson, Brisbourne, et al., 2020) and assuming density (917 kg m −3 ) and seismic velocity (3.841 km s −1 ) at the source (Maurel et al., 2015;Smith, 1997). We then calibrated a local magnitude scale (ML), obtained from the maximum amplitudes in the waveforms, to the M w scale (see processing windows in Figure 2a and processing details in supplementary material S4). We choose this two-step approach as the M w calculation is most accurately conducted only if the focal mechanism is known and if a Brune model (Brune, 1970) can be fit to the displacement spectrum, whereas ML can be calculated for all events in our catalog.
The derived local magnitude scale is based on Smith, Smith, et al. (2015), whose ML for RIS is adapted from the well-established but empirically derived Richter scale (Richter, 1935) and follows the equation: where A is the maximum amplitude of either of the two horizontal components (in instrument counts; all instruments were corrected to a consistent "counts" scale). The distance term m accounts for the decay of amplitudes with increasing epicentral distance (d epi ) and t is a scaling parameter that bridges the offset between M w and ML. We derive ML or M w , for all stations of an event separately. The final magnitude of an event is then calculated as the median of all single-station measurements. The uncertainties are derived as the mean absolute deviation (MAD) of the single-station ML or M w , values from the median. Further processing steps for ML and M w are detailed in Section S4.
We obtain a 1:1 fit (Pearson correlation coefficient of 0.96, with 1 being a linear fit and 0 being no fit) between M w and ML for the data set used for scaling but also a high correlation (Pearson coefficient of 0.88) when considering additional events with rupture mechanism that were not used when initially deriving the M w -ML scaling ( Figure S5). This confirms that a linear M w -ML scaling is adequate to fit our data set (Butcher et al., 2020) and that the relatively simple approach of calculating ML is sufficient. This is likely because the total range of observed magnitudes spans only ∼1.7 magnitude units ( Figure 2b) and because picks from many different azimuths are available for each event.

Event Focal Mechanisms and Stress Inversion
We determine fault plane solutions from first motion polarities and P to S amplitude ratios using the HASH software ( et al., 2018). As the P-wave onset of RIS microseismicity is impulsive and the signal-to-noise ratio is high (Figures 2a and 3), an automated gradient-based polarity picker is implemented (see Section S5 for details on processing approach). Take-off angles are derived from the same velocity model used for the NonLinLoc relocations (a two-layer model of firn and ice). To account for errors in the polarity picks, 15% outliers (nonmatching polarities in the final solution) are allowed during the inversion. We further perform multiple inversions while perturbing take-off angles (standard deviation of 5°) to allow for uncertainties in the velocity model and the event location. The final set of good solutions is derived based on quality criteria, which are the stability of solution upon variations of input, the azimuthal gap of the final set of stations used (should be smaller than 180°) and the final number of input picks (should be larger than seven). Due to the clear waveforms, we derive stable solutions for events in the center of the network domain ( Figure 3a; 52% of all events with backazimuthal gap smaller than 90° have a rupture mechanism solution), but also for events at its extremities (Figures 3b and 3c; 28% of all events with backazimuthal gap between 90° and 180° have a solution).
In addition to single-event solutions, we calculate cluster-wise stress tensors using all individual focal mechanisms of a cluster as input data. The stress inversion is conducted using the software slick (Michael, 1987). Slick performs a linear inversion to minimize the number of rotations around an arbitrary axis necessary to rotate the input focal mechanisms to fit a uniform stress tensor. We assess the quality of the cluster-wise solutions via bootstrap tests. In these tests, the data are resampled 100 times while the fault and auxiliary plane are exchanged for 10% of all input mechanisms. The spread of the results obtained from bootstrap inversions provides a measure of inversion robustness. We only use clusters for which more than seven mechanisms are available. Events are generally well constrained with an average horizontal standard error (as defined by Lee & Lahr, 1972) derived from the pick uncertainties of 27 and 26 m in east-west and north-south, respectively. The mean vertical standard error is 48 m. In addition to these formal errors, we expect a perturbation of the hypocenters of ∼5 m due to errors in the velocity model when considering the uncertainty given by Smith (1997). In addition, a more severe hypocenter perturbation arises due to seismic anisotropy. This may introduce an error of ∼10-20 m horizontally and up to ∼100 m vertically. Last, laterally variable errors in the velocity model may be introduced as the firn-ice transition likely forms a gradual, rather than a sharp, velocity increase (for more details on uncertainty derivation see Sections S1 and S2 and Figure S6).

Spatial Icequake Distribution and Magnitudes
All events cluster near the ice-bed interface, the depth of which is derived from radar data (King, Pritchard, et al., 2016). On average, the events locate 16 m below the interface, which is within the average vertical location uncertainty derived here and the absolute location error of 10-20 m given by King, Pritchard, et al. (2016). The difference may result from velocity variations within the ice, which are not captured with the two-layer model used here or from uncertainties in the absolute reference frame used by King, Pritchard, et al. (2016), or both.
Despite their common depth location, the events show a discontinuous spatial distribution across the study region. Most events, including the largest in our data set, locate either at the boundary between the soft and stiff sediment regions or further downstream in the region of stiff sediment. Within the stiff sediments, there appears to be no correlation between event density and the location of hummocks identified from radar data (King, Pritchard, et al., 2016; see geographic labels in Figures 1b and 4b for orientation). However, more events occur southwest of the central high than northeast of it. In regions of very high seismic activity, KUFNER ET AL. events partly appear to arrange in distinct regions (∼300-500 m radius), which are seismically active at their rims and aseismic in their centers (e.g., Figure 6a). This configuration is robust, even when considering the event hypocenter uncertainties. Seismicity across the transition from soft to stiffer sediments correlates with a step-up in bed elevation across the boundary. If this step is large (e.g., 20 m residual elevation increase in Figure 5b), seismicity is most pronounced and large magnitude events occur, whereas negligible seismicity is associated with a transition without a change in residual bed elevation (e.g., Figures 5a and 5c). Events upstream of the bed character boundary tend to occur in the troughs separating MSGLs, while seismicity at the MSGL crests is absent (Figures 5d-5f).

Event Cluster Characteristics and Rupture Mechanisms
In addition to this large-scale icequake distribution, we observe a small-scale structure in the spatial distribution of most icequakes. Icequakes rarely occur as single events in space and time. Instead, seismicity is focused on spatially isolated spots of less than 100 m radius (highlighted in Figure 6a, zoom in one of these spots in Figure 6b), which produce many icequakes over a short timescale. We use DBSCAN to isolate these spots of focused activity, finding 828 spatial clusters with eight or more events. These clusters include KUFNER ET AL.
10.1029/2020JF006001 9 of 23  188,174 events with magnitudes larger than −1.5, which means that 93% of all events with magnitudes larger −1.5 are clustered. In the following discussion, the term "cluster" will be used to refer to discrete, spatially restricted, sites of icequake activity. By contrast, a "temporal sub-cluster" refers to a temporally limited period of high seismic activity at the cluster location. Within temporal sub-clusters, inter-event times are in the minute range, whereas the time of quiescence between two temporal sub-clusters is on average 3.7 days (see also Section S3 and Figures S7 and S8 for more detailed cluster characteristics).
Most clusters exhibit a common behavior regarding their rupture mechanisms and regarding their magnitude evolution with time: a) Events within one cluster feature highly similar rupture mechanisms (example in Figures 6b and 6c), resulting in well constrained cluster-averaged stress tensors (Figure 6d). We obtain stress tensors for KUFNER ET AL.
10.1029/2020JF006001 10 of 23 428 clusters, which comprise in total 70,023 individual mechanisms. The average spread value in these clusters is 4.6° (min./max.: 1.7°/17.9°). The spread is a measure of how well individual mechanisms match the resulting stress tensor. This small spread value indicates highly similar focal mechanisms within each cluster b) We observe a modulation of event magnitude within the clusters. Most ML-activity-time plots show a short-term increase and decrease of event magnitude with time ( Figure 7). On average, these activity cycles last from 2 to 6 hours (mean ∼3.5 h) and several of these activity cycles may occur in succession (e.g., Figure 7a). This relation is still valid when considering the uncertainty in magnitude (Figure S9). We observe these magnitude patterns for all clusters with a magnitude range larger than ∼0.4 ( Figure  associated to an icequake cluster. The mean P-axes azimuth of all stress inversions is 144 ± 12°. This is comparable to the surface ice flow direction (azimuth of 148°) measured with GPS, which suggests flow-parallel sliding at the base of the ice stream. This agrees with previous source mechanism observations at RIS Smith, Smith, et al., 2015). Here, the P-axes describe a gentle rotation (±11°) toward the ice stream margin on either site of the central high along with this large-scale trend. In addition, we observe a larger rotation (±36°) relative to flow for 5% (23 clusters) of all mechanisms (Figure 8a). These mechanisms primarily occur across the bed character boundary and indicate sliding along the base but at an oblique angle relative to ice flow b) Clusters show three distinct types of spatiotemporal behavior. Most clusters (81%) are active only once for few days (typically <5 days) while only a smaller percentage (19%) is active multiple times for few days. These repeatedly active clusters can be grouped into either "spatially stable" clusters (9%), which always occur at the same geographical position, albeit the ice moves above them, or in clusters where icequake hypocenters migrate downstream with ice stream flow velocity (10%; Figure 8b). The occurrence of this different spatiotemporal behavior is largely independent of their location with respect to the bed character boundary, the number of events in a cluster, or cluster duration (Figures 8c/8d). At the same time, there is no correlation between the number of events and duration of activity ( Figure S7).
Last, we note that the clusters appear to range in shapes and sizes in map view (e.g., circular or elongated). However, we do not consider these variations here. Determining the exact dimension and shape of these KUFNER ET AL.

10.1029/2020JF006001
12 of 23 features is at the edge of the resolution capacity of this icequake catalog. If double-difference relocation methods were used, the single clusters might collapse to even more concentrated features.

Temporal Evolution of Icequake Activity
Despite these pronounced spatial variations, the entire microseismic data set shows little overall systematic temporal variation in activity, nor a strong correlation with daily or biweekly periodicities in the tidal signal at the grounding line ( Figure 9). Instead, we observe the total number of detected events and the cluster KUFNER ET AL.

10.1029/2020JF006001
14 of 23  Table S2 lists all time windows used to derive these plots).
onset times to be dependent on the weather conditions at RIS (Figures 9a-9c, 9e). During periods of strong wind, noise levels are higher and therefore fewer events are detected. However, Figure 9 shows that ∼2 months (January/February 2019) of our data were acquired during stable weather conditions and with a consistent network geometry. Although the seismic network had been active during scientific drilling in the same area (A. M. Smith, Anker, et al., 2020), we do not observed a notable spatial or temporal correlation of icequakes in our catalog and the periods of drilling (5-8th January, 18-22nd January, 6-11th February; Anker et al., 2021).
During the stable weather period, we note that a weak correlation with biweekly tidal maxima might exist when considering only events from the larger magnitude population (Figure 9d; events with b = 10.9 in Figure 2b). For the period of stable weather conditions, the peaks in this histogram vaguely correspond to the temporal positions of the neap tides. The mean time difference of the events to the closest neap tide is 1.7 days (std of 1.0 days), whereas the mean time difference to the closest spring tide is 5.8 days (std of 1.7 days). However, apart from this weak correlation, Figure 9f illustrates the near chaotic temporal behavior of the event clusters. For instance, the four largest clusters in our data set behave completely differently with time. Whereas one of the clusters produces all events during ∼5 days of intense activity, the other clusters are split into several temporal sub-clusters with varying numbers of events and activity times. Neither of these clusters or sub-clusters correlates with daily or biweekly trends in the tidal signal. Finally, the largest clusters (those with the highest number of individual events) do not necessarily include the events with the largest magnitude.

Little Influence of Tidal Forcing on Icequakes
Tidally induced sea-level modulations at the grounding line have been shown to cause periodic coupling and decoupling of the ice and bed at other West Antarctic ice streams. For example, at Kamb Ice Stream (KIS), such modulations have been measured as far as 85 km upstream of the grounding line (Anandakrishnan & Alley, 1997). A likely tidally induced modulation of ice flow speed, yielding a temporally variable icequake rate has been observed at Whillans Ice Plain (WIP) (Barcheck, Schwartz, et al., 2020;Barcheck, Tulaczyk, et al., 2018;Winberry et al., 2013) as well. Also, at RIS, a biweekly modulation of ice flow velocity is observed at the surface (Adalgeirsdóttir et al., 2008;Gudmundsson, 2006;Murray et al., 2007). However, in our icequake data set, we observe only a weak correlation between the occurrence of the largest magnitude events and the neap tide at the grounding line ( Figure 9d). This is when the glacier flow velocity at the surface reduces (Gudmundsson, 2006). A similar weak correlation of larger magnitude icequakes with the tidal cycle has been observed before (Adalgeirsdóttir et al., 2008;Smith, Smith et al., 2015). Although this trend is weak, it has been observed in three different data sets, collected in different years (1997-1998, 2008-2009, and 2018-2019), so is likely a characteristic of basal microseismicity at RIS. Although magnitude variation can be caused by a number of other factors (e.g., fault size, variable loading velocity), a temporal modulation of icequake magnitudes could suggest that the pressure regime at the ice-bed interface changes temporally. The two different b-value trends, we observe for our data set (Figure 2b) also allude to this. Variable b-values can occur due to changes in the stress regime during the observation period (El-Isa & Eaton, 2014) or during the transition from tectonic to fluid assisted failure (Kettlety et al., 2019).
However, apart from the possibility of a gentle magnitude modulation with the tidal cycle, the bulk of basal seismicity does not show any clear biweekly trend. A correlation between icequake intensity and daily tidal height at the grounding line is not observed either. Thus, our observations contrast with the results at KIS or WIP, where variations of seismicity rates with one or two daily peaks, often but not always associated to the tidal cycle, have been observed (Anandakrishnan & Alley, 1997;Pratt et al., 2014). At KIS, basal icequakes are thought to accommodate significant parts of the basal ice stream motion (Anandakrishnan & Alley, 1997). At WIP, microseismic events likely indicate the nucleation phase for a tidally induced largescale movement of the ice stream (Winberry et al., 2013). Furthermore, Barcheck, Schwartz, et al. (2020) inferred an alignment of seismicity and MSGLs. This basal seismicity is periodic and influenced by glacial flow variations, likely produced by the tidal cycle.
The situation at RIS is clearly different and prompts two possible interpretations. On one hand, our observations could be explained by a scenario where basal sliding varies temporally, like the observed surface modulation, but is not reflected by stick-slip seismicity at the bed. Thus, icequakes would make up only a small proportion of the total motion and the tidal signal could be accommodated by aseismic deformation and movement at the bed (Smith, Smith et al., 2015). On the other hand, there could be intra-ice deformation at RIS, which modulates the deformation signal from the surface to the bed. Such tidally induced modulations in the vertical strain rate have recently been detected in ice sheets (Vankova & Nicholls, 2019). Furthermore, the ice at RIS is much thicker than at WIP (2,200 m at RIS, 650-800 m at WIP; Fretwell et al., 2013), which could explain why intra-ice deformation has a larger impact at RIS. In addition, at RIS, icequakes along the MSGLs occur with a similar spatial distribution to WIP but without a clear temporal pattern. This discrepancy would be an argument for tidal forcing at the base of RIS being less pronounced than at other ice streams. In addition, the bed topography of RIS is more extreme (see e.g., Figure 1c), compared to the relatively flat bed of WIP, which might hinder the upstream propagation of a tidally induced pressure change along the ice-bed interface. On the contrary, in modeling studies, tidal forcing has been suggested to periodically modulate the surface ice flow via friction at the bed (Rosier et al., 2015). This would require the tidal signal at the bed to be even more pronounced than at the surface and would be an argument for dominantly aseismic motion at the bed of RIS. Ultimately, measurements that monitor the strain or fabric modulations through the ice column might help to discriminated between these different scenarios. However, in either case, this study shows that the basal seismicity at RIS is not, to first order, controlled by the tidal cycle.

Network Wide Icequake Distribution: Role of Bed Topography, Bed Properties and Water in Triggering Icequakes
As the tidal influence in icequake distribution appears minor, other characteristics, such as bed topography and till properties at the bed, must have a greater impact on temporal and spatial icequake distribution. Soft till will accommodate ice flow by deformation whereas stiff till favors basal sliding (Smith, 1997). Accordingly, and in agreement with previous icequake studies at RIS (Smith, 2006;Smith, Smith et al., 2015), we observe more icequakes within the stiffer bed domain than in the soft sediment units (see Maps in Figures 1c and 4b for the geographic locations used in the following discussion). However, due to the superior network configuration and size compared to previous studies, we observe previously unresolved second order structures. We observe, for instance, fewer icequakes northeast of the central high than southwest of it within the broad domain of stiff sediments downstream of the bed character boundary. A lower radar reflectivity has been inferred for the latter, which suggests outcropping bedrock or very compressed sediments (Schlegel et al., 2021). Thus, larger regions of reduced seismicity within stiff till units could indicate the presence of compressed sediment or outcropping bedrock.
Furthermore, the icequake distribution highlights features that we suggest indicate variations in bed character over scales of hundreds of meters. This can best be illustrated within the soft sediment upstream of the bed character boundary, where we observe large-scale flow-parallel alignment of events within the valleys separating MSGLs (e.g., Figure 5d). This suggests that stiff till must be present in the valleys to favor seismogenic stick-slip behavior. Barcheck, Schwartz, et al. (2020) inferred similar alignment of seismicity and MSGLs at WIP. They also related these patterns to changes in frictional properties (soft sediments on top, stiffer at base). This alignment could be due to either the constructional or erosive creation of MSGLs. In both cases, soft sediments would be expected at the MSGLs crests. However, at RIS, seismic studies showed that MSGLs likely form when soft deformable sediments are accumulated . Thus, the similar character of seismicity at WIP and RIS might hint toward a constructional creation of MSGLs in general.
Another bed feature at the scale of a few hundreds of meters highlighted by our icequake catalog is the seismicity arranged in circular patterns with the centers depleted in seismicity within the broad domain of stiff sediments (e.g., Figure 6a). At least one of these central regions corresponds to an area where free water is proposed to exist at the glacier bed (Schlegel et al., 2021). However, we can rule out a direct role of fluid in creating icequakes through tensile crack opening as the RIS icequakes are likely caused by a double-couple source. Icequakes directly triggered through the hydraulic system at the glacial bed may manifest themselves through nondouble-couple tensile crack faulting (Walter, Canassy, et al., 2013). We infer the double-couple nature of the RIS icequakes from the station coverage that allows for many rupture mechanisms the coverage of the entire focal sphere (e.g., Figure 3a). If the icequake source would have a significant nondouble-couple component, a less clear separation of positive and negative polarities close to the nodal planes would be expected. In addition, results of full-waveform modeling for one icequake at RIS show a double-couple source to be more likely . Full-waveform modeling would allow for the resolution of different source types. Thus, both studies suggest that neither the direct role of fluids (e.g., through hydrofracturing) or other processes that require tensile forces (e.g., crevasse opening) seem to drive the icequakes at RIS. Instead, we propose that a weakening of the till resulting from the presence of fluid penetrating the till layer (Rathbun et al., 2008) or fluctuations in the hydraulic pressure caused by fluids (Röösli et al., 2016) may eventually result in a series of stick-slip events adjacent to regions of free water at the ice-bed interface. The role of fluids in promoting icequakes would also explain the temporal and spatial event clustering we observe (Fischer, Horálek, et al., 2014;Greenfield et al., 2019) and the large b-values (El-Isa & Eaton, 2014;Schlaphorst et al., 2017;Wilks et al., 2017). Large b-values are indicative for swarm-like earthquake behavior, which, in turn can be triggered from fluid induced pressure variations.

Zooming into Individual Icequake Clusters: Types of Subglacial Stick-Slip Deformation
Icequakes typically occur clustered in space and time, as observed at the bed of other glaciers in Antarctica and elsewhere (e.g., Barcheck, Schwartz, et al., 2020;Danesi et al., 2007;Helmstetter, Nicolas, et al., 2015). Despite cluster nature, size and repeat time being highly variable, clustered icequake activity is generally interpreted to be caused by sticky-spots (Barcheck, Schwartz, et al., 2020;Fischer, Clarke, et al., 1999;Jacobel et al., 2009;Röösli et al., 2016;Smith, Smith, et al., 2015), where basal resistance increases. Although descriptive, this interpretation does not necessarily imply a specific physical mechanism. Also, at RIS, we observe a large spatiotemporal variability of cluster nature. Due to the relatively long observation period and large network aperture compared to previous studies, as well as the detailed knowledge of bed properties from seismics, drilling and radar, we suggest that sticky-spots at RIS can be attributed to three different end-member types of stick-slip behavior at the glacier bed, which are schematically shown in Figure 10. We note, however, that there is likely no strict separation between these different end-member types and they KUFNER ET AL.

10.1029/2020JF006001
17 of 23 Figure 10. Schematic interpretation sketch on active basal processes. The loci of icequakes depend strongly on bed type, with most events occurring within the stiffer sediments. Different processes can trigger the icequakes. Among these, the continuous creation of sites of increased friction that develop due to sediment transport and/or due to temporal variable till properties within the bed (asperities) and their seismogenic destruction is most common.
can occur simultaneously or intermingled. As a whole they may be indicative of the deformation characteristics of subglacial till beds.
In the following interpretation, it is assumed that all icequakes occur very close to the ice-bed interface. This assumption is justified as the vertical location uncertainty ( Figure S6), including possible effects of model errors and anisotropy (Sections S1 and S2), and the uncertainty in the radar-constrained interface (King, Pritchard, et al., 2016), places all events at the interface. This agrees with full-waveform source inversions that suggest that such icequakes at RIS occur within meters of the ice-bed interface . Furthermore, we note that our event catalog does not allow us to draw detailed conclusions on the shape of the individual clusters. The event cluster size is generally in the 10-100 m range and with variable shape. However, based on the location errors derived from the pick uncertainties and additional uncertainty from unmodeled errors in the velocity model (isotropic and due to anisotropy), it is likely that the different shapes of individual clusters may be within the location uncertainty. Thus, all events in one temporal sub-cluster may originate from a single spatial location, that is, a single fault.

Type 1-Self-Destructive Asperities
Most of the icequake clusters (81%) are active for less than 5 days (mean 3.5 days, std 8.4 days; Figures 8d and S7a). During this time, the ice stream flows ∼3.7 m downstream. We detect icequakes with inter-event times in the 1-5 min range (mean 4.6 min, std 9.3 min; Figure S7d). These clusters are then inactive for the remainder of our observation period, which suggests that these spots are unlikely to be stationary obstacles in the ice stream bed. Stationary obstacles would likely produce repeating seismicity, for example, upon variations in basal water pressure (Fischer, Clarke, et al., 1999) or due to constant ice loading. Instead, we favor a concept of asperities within the subglacial till, which are randomly built by the glacial movement and subsequently destroyed through a sequence of stick-slip events. Such asperities may be envisaged as sites of increased friction that develop during continuous ice stream movement as sediment is transported and dilates and reorganizes (McBrearty et al., 2020;Thornsteinsson & Raymond, 2000;Van Der Meer et al., 2003). If glacial till is sheared, its pore volume is increased (Boulton & Hindmarsh, 1987). Till can then be weakened if pathways that permit water flow into the dilated material exist (Rathbun et al., 2008). This may lead to the formation of an asperity along which slip-deficit can build-up. Freezing-on of part of the bed could additionally contribute and would favor velocity weakening (Lipovsky et al., 2019), which is a requirement for stick-slip behavior. Once the shear resistance of the asperity is overcome, stress is released in a series of icequakes and the specific asperity is destroyed. The displacement per event at RIS is estimated to be in the range of 0.03-0.07 mm (Smith, Smith, et al., 2015), which is less than the glacial movement that would accumulate during typical interevent times (∼0.8 mm min −1 , assuming the same velocity at the base of the ice stream and the surface). Thus, it is unlikely that each new event in an icequake cluster is created by continuous loading. It rather suggests that a spot mostly deforms aseismically but slip-deficit can accumulate occasionally. As we rarely observe single events, but event clusters, it appears that glacial till does not support the accumulated slip-deficit to be released in a single large event (e.g., comparable to megathrust earthquakes in subduction zones), but rather in many small icequakes. If the asperities develop due to the ice stream movement and reorganization of till, different event counts per clusters, as observed here, can be envisaged. The sharp magnitude cut-off at larger magnitudes (b-value of 10.9) obtained here might also suggest that an upper magnitude threshold exists for the largest possible icequake. This magnitude threshold may be governed by the till properties and the maximum available normal stress.
In this concept of self-destructive asperities, the bed material must be strong enough to allow for the buildup of stress locally. This may explain why more icequake clusters occur in the stiff-sediment domain. Furthermore, the bed character boundary sections with a large step in residual topography may be favored for the occurrence of such clusters as they represent natural obstacles for flow. On the contrary, it seems that very stiff surfaces, like the stiff-sediment units northeast of the central high, are less favorable for asperity formation. This may be as sediment reorganization is expected to happen more slowly. Instead, they may give rise to polished surfaces, possibly overlain by a homogeneous water film, where aseismic glacial sliding is the dominant basal motion process.

Type 2-Ploughed Clasts
For some clusters (numbering 72; 9% of all clusters), we observe the downstream migration of the seismically active sites at the same speed as ice flow at the surface, which is ∼1.05 m day −1 (Figure 8b). This phenomenon occurs for ∼50% of all clusters which are active for a sufficient duration that the observed migration is larger than the single event location uncertainty. This observation suggests that an object, held within the ice, is being transported downstream and causing the icequakes. During this transportation process, the spot is periodically seismically active. Likely candidates for such a mobile object are clasts held in the basal ice and dragged through the glacial sediments or over harder materials (Zoet & Iverson, 2020). The presence of clasts embedded in the bed had been proposed based on scientific drilling at RIS (A. M. Smith, Anker, et al., 2020). If clast motion is hindered for some time, allowing slip-deficit to accumulate, the seismic activity could represent the moment in which the clast slips forward. An icequake with double-couple source would then be created by frictional sliding between the clast and/or the ice and till layer. Laboratory experiments showed that ploughing clasts can cause velocity weakening behavior (Iverson, 2011;Thomason & Iverson, 2008). The clast may eventually become lodged due to melt out or changes in the properties of the sediment. Such clasts will have variable shape, size and penetration depth, and so different numbers of events in the clusters appears logical. Our event catalog does not allow us to comment on the size or shape of such clasts, as we consider them to be within the horizontal resolution of the event locations. In contrast to icequakes originating from breaking asperities, bed deformation is expected in the case of ploughed clasts (Zoet & Iverson, 2020).
Apart from downstream migrating clusters, we observe some clusters (numbering 82; 10% of all clusters) that are active repeatedly at the same location ( Figure 8b). These could represent the presence of a more permanent obstacle to ice flow. Either basal drag could be too weak or the till matrix too strong to allow for the mobilization of a clast. Alternatively, these clusters could be related to bed asperities. Part of an asperity may remain locked after the initial cascade of icequakes and break at a later stage.

Type 3-Flow-Oblique Landforms as Obstacles
Our stress inversion data set contains 23 clusters (5% of all clusters with stress inversion results) that indicate flow-oblique deformation (Figure 8a). This rotation (±36°) is clearly supported by the data. For instance, seismic stations, crucial for constraining the rupture mechanism, show different polarities for either flow-parallel or flow-oblique mechanisms (Figures 3b and 3c). Furthermore, the rotated events occur close to mechanisms that are not rotated. Thus, their occurrence is unlikely to be an effect of network geometry. Such flow-oblique mechanisms have not been observed at RIS before and we suspect that it is the dense seismic network and the low-noise level that allows them to be resolved here. Based on our first-motion mechanisms, we infer that the main difference between these icequakes and those discussed above is the rotation of the strike of the rupture mechanism. A comparative analysis of the source characteristics of the different icequake populations might yield further discrimination.
Although scenarios can be envisaged where rotated rupture mechanisms originate from self-destructive asperities of ploughed clasts, it is striking that these mechanisms mainly occur along the bed character boundary and at the termination points of MSGLs. This suggests a causal relationship. These flow-oblique focal mechanisms may be related to intra-till deformation that occurs during the formation of subglacial landforms, either at the ice-bed interface or within the deforming till. This agrees with laboratory experiments conducted by Lipovsky et al. (2019), who concluded that shear seismicity may indicate geomorphological activity. Here, the flow-oblique mechanisms occur mostly along the bed character boundary. The bed character boundary is thought to be modified over time scales of a few years by sediment deformation (King, Pritchard, et al., 2016;Smith & Murray, 2009). Furthermore, the flow-oblique mechanisms tend to focus along the termination points of MSGLs, where active erosion and deposition has been interpreted from seismic data . Smith, Murray, et al. (2007) concluded that sedimentary processes may be the most likely explanation for this erosion. The flow-oblique focal mechanisms are likely the brittle manifestation of such sedimentary processes. Alternatively, or in addition, the flow-oblique focal mechanisms may originate at outcropping bedrock. Such bedrock units cannot be eroded and may form an obstacle that creates a local distortion of the stress regime.

Summary and Outlook
We present a microseismic event catalog for a 10 × 10 km region, ∼40 km upstream of the grounding line of RIS. The seismic network used to derive this catalog straddles a change in bed character properties (soft to stiff sediments) with consistent station spacing. Thus, we can identify seismic and aseismic regions within our network domain with high certainty.
All ∼230,000 micro-earthquakes (magnitudes between −2.0 and −0.3) detected in a 90-day observation period are located near the ice-bed interface. Most of these events indicate flow-parallel stick-slip. We propose that the interplay between the topography, bed character type, and the hydraulic system at the bed controls the spatiotemporal patterns in icequake occurrence. Icequakes focus at the transition from soft to stiff till and in defined spatial domains of stiffer till. The domains within stiffer till can be either large, coherent regions or more subtle structures, like the valleys separating MSGLs. Within the regions of stiffer till, fluids may modulate the strength of the till to promote seismicity. In contrast, tidally induced pressure fluctuations at the bed seem to be less pronounced or have little effect on icequake occurrence. This suggests that part of the tidally induced modulation is taken up by aseismic bed or intra-ice deformation.
On a smaller scale, most icequakes (93%) occur in clusters that are spatially and temporally restricted bursts of seismic activity. Accordingly we measure high b-values (between 3.3 and 10.9) in event number-magnitude plots, which are indicative for swarm-like behavior of earthquakes (El-Isa & Eaton, 2014). Modulations in b-values might be due to pressure fluctuations at the ice-bed interface or indicate an upper limit for the maximum icequake size. These clusters are generally less than ∼100 m in radius and are active for only a few days. Based on the calculated location uncertainty, we suspect that all events in a specific cluster could originate from the same spatial spot, that is, a single fault. These clusters show an increase and decrease of event magnitude with time while the events in a specific cluster feature highly similar rupture mechanisms. We further observe a gentle correlation of increasing inter-event times with increasing magnitudes. Similar icequake characteristics have been observed in very different glacial settings in the European Alps and in Greenland (Helmstetter, Lipovsky, et al., 2018;Helmstetter, Nicolas, et al., 2015;Röösli et al., 2016), although observation time spans in these studies are shorter than in this study. Thus, these common characteristics may provide insight to the rupture mechanism of icequakes, that is, the rupture of an asperity surrounded by aseismic slip, in general. Furthermore, such common characteristics hint toward a unique driving force within a cluster and suggests event triggering within the clusters once the activity period is initiated, possibly facilitated by frictional heating. A detailed investigation of the source mechanisms, the inter-event locations (e.g., through double-difference methods) and of the material properties surrounding the events might help to discriminate between such processes.
Apart from these common features, the clusters can be discriminated from each other based on distinct spatiotemporal evolution characteristics and the orientation of rupture mechanisms relative to ice flow. We attribute their distinct characteristics to different end-member deformation mechanisms that may act at the bed simultaneously. These are the dynamic creation and seismogenic destruction of spots of increased friction that develop due to sediment transport and/or due to temporal variable till properties ("asperities"), the ploughing of clasts through the underlying sediment, and flow-oblique deformation either associated with the erosion and formation of subglacial landforms or due to bedrock obstacles at the ice stream bed. Among these, the seismogenic destruction of asperities is the most common process. Taking these different processes together, we conclude that the bed of RIS can be envisaged as an actively and heterogeneously deforming subglacial bed mosaic (Piotrowski et al., 2004) with a variety of deformation processes active simultaneously. Our analysis suggests that the friction at the bed varies over a small scale and that the glacial bed is in a process of continuous reorganization. Both impact ice stream flow directly.