Fine Structure of Microseismic Glacial Stick‐Slip

Frictional instabilities exist in many geological settings, including glaciers and tectonic plate boundaries. However, investigations of suggested analogies between stick‐slip “icequakes” and earthquake faulting have been hampered by the noisy, melt‐prone and inaccessible nature of glacial environments. Here, we reveal details of stick‐slip events beneath an Alpine glacier using seismic sensors within a few meters of a seismically active bed region. We present evidence that widely detected stick‐slip events, which are measurable at the ice surface, are in fact dynamic ruptures over many smaller asperities, whose individual seismic failures are usually too small to be recorded at the surface. Characteristic recurrence times of such multi‐asperity ruptures and their sizes suggest an analogy to Parkfield earthquakes on the San Andreas Fault, questioning traditional glacier sliding theories. Although several trillion times smaller, glacial seismic sources presented here may therefore be ideal for studying earthquake faulting due to much higher event rates.

stick-slip motion not evenly distributed across the glacier bed but concentrated at irregularly spaced asperities? (b) Do asperities' rupture planes locate at bi-material interfaces such as ice-till, or are they located within homogeneous materials such as ice, till or bedrock? (c) What is the spatial extent of asperities? (d) What are the co-seismic and aseismic displacements during and in between microseismic stick-slip events? (e) Within what radius of influence does an active stick-slip asperity affect basal sliding? Answering these questions is necessary for parameterizing frictional sliding in large scale ice flow models, thus ultimately for predicting future sea level rise.
New observations have established links between seismogenic stick-slip activity and the material characteristics of glacial beds (Kufner et al., 2021;Lipovsky et al., 2019). Other studies have estimated slip distances and fault sizes (Anandakrishnan & Alley, 1994;Helmstetter et al., 2015). However, these studies all rely on seismic data recorded at the ice surface. Seismic waves emitted by stick-slip icequakes are therefore prone to significant path effects, and their surface records may contain little or no information about dynamic properties of microseismic ruptures. Accordingly, the questions above are difficult to address with seismic surface records and existing findings have uncertainties, which are unlikely resolvable with further surface measurements.
Here, we use seismic borehole instrumentation within a few meters of an active microseismic stick-slip asperity at an Alpine glacier bed in order to study stick-slip icequakes. Our recordings show the substructure of a basal stick-slip asperity, quantify its extent, the co-seismic slip distance and highlight its relevance for basal sliding.

Study Site and Experimental Setup
Our study site is located in the ablation zone of Rhonegletscher, Switzerland (Figure 1a), where the ice is ∼200 m thick and at the pressure melting point (Rutishauser et al., 2016). Previous surface studies in the same region revealed seismogenic stick-slip motion clustering at distinct asperities at the ice-bed interface (Gräff & Walter, 2021;Hudson et al., 2020;Walter et al., 2020).
Between 2 July and 26 August 2020, we deployed eight Lennartz 3D/BHs seismometers sampled at 500 Hz and placed in shallow boreholes at depths of 4 m. The seismometers formed a concentric array with 400 m diameter, extended by one seismometer 400 m north of the array center (Figure 1b, Table S1 in Supporting Information S1). Similar to previous studies (Gräff & Walter, 2021;Hudson et al., 2020;Walter et al., 2020), we (a) used this array for stick-slip icequake detection by searching the continuous data for high frequency events with distinct P-and S-waves and no Rayleigh wave component (for details see Gräff & Walter, 2021), and (b) located stick-slip asperities at which stick-slip events cluster with the probabilistic nonlinear hypocenter location algorithm NonLinLoc (Lomax et al., 2000, Table S3 in Supporting Information S1). In contrast to previous works, for the present study we searched the data for active basal stick-slip asperities already during the field campaign in order to guide subsequent hot water drilling.
Based on data from the first 2 weeks of the near-surface array deployment, we chose one active asperity (henceforth "main asperity") for further detailed studies, including borehole monitoring. At depths of ∼145 m (∼35 m above the glacier bed), we installed five ASIR AG4.5 deep-borehole geophones, sampled at 2 kHz, forming a circle of 30 m radius centered on the main asperity's hypocenter. The array's location close to the glacier bed reduces seismic wave attenuation from long travel distances, especially within the uppermost crevassed and weathered ice layers (Mikesell et al., 2012;Walter, 2009) and thus overcomes noise-induced detection limitations of surface records (Roeoesli et al., 2016;Walter et al., 2008).
With the waveform data of the near-surface seismometers and deep-borehole geophones, we re-located the basal stick-slip asperity with NonLinLoc to an accuracy of 5 m. This enabled us to guide the hot water drill into the asperity that we expect to cover several square meters (Gräff & Walter, 2021;Walter et al., 2020). On 23 July, we drilled one borehole 1 m south of the best estimate epicenter through 180 m ice until we reached the glacier bed and installed two pressure sensors, from which we could unfortunately only retrieve corrupted measurements due to incorrect electronic wiring. However, this borehole indicated the local ice thickness and so we drilled two additional boreholes, located at the estimated asperity epicenter and 5 m south of it, respectively, to a depth of (175 ± 1) m. These boreholes were each equipped with one IMS tri-axial 25 kHz accelerometer sampled at 48 kHz (IMS1 & IMS2) and located ∼5 m above the glacier bed ( Figure 1c). We expected the sensors to freeze to the ice for a solid coupling, which we further enhanced by embedding the sensors in sand. After retrieving all sensors end of August we used the boreholes for inclinometry.

Near-Surface Measurements
Two stick-slip asperities remaining active for more than 1 week were found during the first 2 weeks of the near-surface seismometer deployment. We decided to focus on the one asperity exhibiting unambiguously dilatational P-wave first motions at seismometers up-glacier and compressive down-glacier (Figures 1b, 1d and 1e), revealing a double-couple focal mechanism that excludes the possibility of a basal tensile crevassing icequake (Walter et al., 2013). Events from the main asperity have a typical recurrence time of 55 +28 −15 min ( Figure 2b) and a seismic moment of (1.0 ± 0.5) × 10 5 Nm determined by Sh-wave integration (Methods section of Gräff & Walter, 2021). Attenuation and undersampling of basal signals recorded at the glacier surface do not allow for calculation of source characteristics (Allstadt & Malone, 2014;Helmstetter et al., 2015).

Deep-Borehole Measurements
Recordings from the deep-borehole geophone array predominantly showed signals from the main asperity, located in the center of the array. Even though we sampled the geophones with 2 kHz, signals from the main asperity were still undersampled ( Figure 1e), pointing toward a frequency content above the Nyquist frequency of 1 kHz. Compared to the spectral peak of the near-surface sensors below 100 Hz, the spectral peak measured above 100 Hz emphasizes the high-frequency attenuating effect of the ice (Figures 1f  and 1g). From combined near-surface and deep-borehole recordings, we infer a double-couple fault mechanism with a slip pointing south to southeast. Using seismic moments calculated with the deep borehole geophones by Sh-wave integration, we determine the source orientation that minimizes the moment discrepancy between stations: Strike ∼97 • aligned with bed contours, dip ∼33 • (ca. 15° steeper than the bed gradient from radar-echo sounding [Rutishauser et al., 2016]), rake ∼−55° pointing toward south-east (beach ball in Figure 2a).
The deep-borehole geophone recordings show weak basal seismic activity near the main asperity, which we investigate with the earthquake detection and clustering algorithm REDPy (Hotovec-Ellis & Jeffries, 2016). REDPy detects repeating signals in continuous seismic data with an STA/LTA trigger (Trnkoczy, 2009) and clusters them based on waveform similarity determined by cross-correlation. We confined this analysis to only 3 days (24-26 July) in order to avoid interfering tremor-like signals present on most other days and data gaps. In addition to the main asperity detected with the near-surface seismometer array, we detect 28 "weak clusters" of basal seismicity that locate within ∼10 m of the main asperity and further basal seismicity in up to 150 m distance. In the near-surface records, signals of these weak clusters are missing or weak enough to evade detection ( Figure S1 in Supporting Information S1). Systematic differences of ∼1 ms in P-arrivals and few milliseconds in S-wave arrivals, respectively, indicate that weak clusters themselves locate at distinct individual asperities. We relocate these "weak asperities" relative to the main asperity with the double-difference earthquake location algorithm HypoDD (Waldhauser & Ellsworth, 2000). They locate south-east of  Figure 1b) with focal mechanism as beach ball. The locations of the two accelerometers (red triangles) and the rupture initiation of the main asperity (orange ellipse) are given. Blue regions mark weak asperities with their location uncertainties. (b) Recurrence time to seismic moment scaling of main asperity (orange) and weak asperities (blue). Faint scatter points indicate single events. Solid markers with error bars indicate asperities' median recurrence time and seismic moment within the 16th and 84th percentile. The marker size represents the number of events from an asperity. Inset shows the normalized cumulative number of events greater than the magnitude on the horizontal axis.
Blue bars indicate events from weak asperities, orange from the main asperity.
the main asperity's point of rupture initiation (Figure 2a), whereas the surroundings of these weak asperities do not produce basal seismicity.
Seismic moments of events from weak asperities are about 60× smaller than for the main asperity. Their median seismic moment is 2.0 +3.8 −0.8 × 10 3 Nm, whereas within the same time span it is 1.2 +0.5 −0.3 × 10 5 Nm for the main asperity. The event moment magnitudes from weak asperities follow a Gutenberg-Richter law (Gutenberg & Richter, 1955) with a b-value of b = 1.73 ± 0.05, whereas there is a moment gap between events from weak asperities and the main asperity, breaking the scale invariance of the Gutenberg-Richter law (Figure 2b and Figure S2 in Supporting Information S1). Within the 3-day time window, recurrence times of weak asperities are widely spread at 15 +70 −12 min and thus overlap with those of the main asperity (65 +35 −12 min) within 1σ (Figure 2b).

In Situ Measurements
The two accelerometers located 5 m above the glacier bed recorded stick-slip events from the main asperity for 6 days until the asperity's activity stopped (31 July) and one additional event 8 days afterward. This latter event had a seismic moment 20 times smaller than previous events. The sensor in the southern borehole (IMS2) recorded 112 stick-slip events from the main asperity, whereas we only captured two events with the northern sensor (IMS1) as the digitizer unexpectedly overwrote almost all recorded events. One of those two events recorded on IMS1 occurred only a few minutes after installation, that is, before the sensor was embedded in sand. The other event is the last and weak event.

Acceleration Waveforms
In the acceleration waveforms (Figure 3a), we interpret the first clear arrival to be the direct P-wave and a second lower frequency phase arrival around 2 ms later to be the S-wave. The event's signal arrives first at sensor IMS2 and (0.52 ± 0.02) ms later at IMS1, revealing that IMS2 must be closer to the point of rupture initiation. Diminishing oscillations in event waveforms when the sensors were embedded in sand and probably frozen to the ice provide evidence that oscillation in the signal coda are caused by a weak sensor coupling ( Figure S3 in Supporting Information S1).
We locate the rupture onset of stick-slip events at the main asperity relative to the two accelerometers with NonLinLoc by picking the first motion of the P-and S-waves and constrain it by the travel direction of the P-wave. This results in an unambiguous location at the glacier bed (i.e., not above the sensors). Figure 2a shows the relative location of the rupture initiation and the two accelerometers. A relative distance of (8.0 ± 1.0) m for IMS1 and of (5.5 ± 1.0) m for IMS2 indicates that for frequencies in the upper range of what we measure ( 10 kHz), our sensors are in the intermediate radiation field-even for source extents of only ∼1 m (Equation 10.12 in Aki & Richards, 2002).

Spectral Content
In view of our sensors' locations in the intermediate field, we compare the acceleration amplitude spectra of the first stick-slip event, which was recorded on both accelerometers (IMS1 and IMS2). The strongest spectral peak which resides at approximately 700 Hz matches for both sensors and represents the oscillations in the coda. Recordings of IMS1 (north of the rupture initiation and at greater distances to it), contain suppressed high frequencies above 2,000-3,000 Hz compared to IMS2 located almost directly above the rupture initiation point (Figure 3b). Since the near-field affects low frequencies (Madariaga et al., 2019), and since attenuation, which is stronger for higher frequencies cannot not explain the contrast between the concave high frequency signal at IMS2 and the convex signal at IMS1, we suggest that a Doppler effect is responsible for this discrepancy (Bakun et al., 1978;Pacor et al., 2016). Higher frequency content for IMS2 is expected to arise for a rupture front propagating south or south-east, away from accelerometer IMS1 whilst laterally or even slightly toward IMS2 ( Figure S4 in Supporting Information S1). Strong rupture directivity in slip direction is predicted for faults with high bi-material contrasts (Shi & Ben-Zion, 2006). In a glaciological context, such a contrast is given at the ice-till interface (Iverson, 1999), which is known to favor stick-slip for sediments entrained in the ice sole . This suggests, that the fault plane of the main asperity is located on such a bi-material surface, rather than within homogeneous material such as ice, till or bedrock which is supported by borehole imaging ( Figure S13 in Supporting Information S1) and in agreement with previous studies Walter et al., 2020).

Residual Displacement
The intermediate radiation field records of our accelerometers inhibit conventional calculation of source properties from corner frequencies (Brune, 1970;Haskell, 1964;Savage, 1972). However, they allow measurements of a residual displacements due to the co-seismic slip of the glacier over its bed. From the acceleration recordings, we determine the residual displacement using an integration method minimizing sensor drift (Boore et al., 2002). On all components, displacement waveforms show a clear offset from zero-displacement after the stick-slip event (Figure 3c). The strongest residual displacement resides on the north components. It has a negative sign, in accordance with a roughly southward slip direction. Considering all 112 stick-slip events recorded at accelerometer IMS2, we observe a systematic correlation between residual displacement and seismic moment determined independently with the deep-borehole geophone array (Figure 3d). In contrast, we find only a weak dependence between inter-event time and residual displacement ( Figure S6 in Supporting Information S1).
During stick-slip events, the shear stress built up at asperities is partially removed by the seismic rupture, leading to a local stress drop. Integrated over the ruptured area, the stress-drop sums up to the resistive shear force acting on the glacier ice and building up elastic deformation due to gravitational loading between stick-slip events. In order to determine this resistive shear force, we model the elastic response of the glacier ice above the main stick-slip asperity. We perform a Bayesian inversion (Foreman-Mackey et al., 2013) of the residual displacement measurements considering two cases: first for a resistive point force acting within a full space (Kelvin problem), and second for a resistive point force acting tangentially on a half space (Cerruti problem) with the elastic properties of ice (Ling et al., 2002, and Figure S8 in Supporting Information S1). The Kelvin and Cerruti Problems reflect the end member cases of (a) a strong coupling at the ice-bed interface (e.g., a frozen bed) and (b) a completely free interface surrounding the point where the force acts on the ice (e.g., weak subglacial till or basal cavities), respectively. Because the resistive point force models simplify the real seismic rupture at the bed and thus do not represent three-dimensional displacements adequately, we only fit the absolute residual displacements measured at the two accelerometers. The inverted resistive point force is F = 0.06-0.12 MN, twice as high for the isotropic space as for the half space (the isotropic space can be approximately regarded as two coupled half spaces; Figure S7 in Supporting Information S1). From this resistive force and the seismic moment that we determined by Sh-wave integration of the deep-borehole geophone recordings, we calculate the fault area to (67 ± 35) m 2 (Supporting Information S1) by applying the Cerruti problem and therefore assuming (a) a circular fault, (b) significantly lower shear resistance surrounding the asperity, and (c) a weak bed at the fault extent. We favor this setting as the Doppler effect in the accelerometer data testifies to a rupture propagating along an interface of high material contrast, that is, a weak bed. However, the calculated fault area only constitutes a rough estimate, due to the simplicity of the point force model. With the calculated area, and assuming an equally weighted sum of ice-side and till-side partial moments (Ampuero & Dahlen, 2005;Wu & Chen, 2003), the ice-side slip D i ≤ (8 ± 3) μm (Supporting Information S1). For an ice-till interface, the total slip can be multiple orders of magnitude larger due to small shear moduli of weak unconsolidated subglacial till resulting in more elastic deformation within the till compared to the ice (Lipovsky et al., 2019, and references therein). For a weak (till) bed with shear moduli in the range of several tens of MPa , the total co-seismic slip may reach millimeters per stick-slip event, summing up to centimeters per day. This magnitude becomes relevant for local basal sliding, and can even dominate local basal motion over viscous creep and regelation, indicated by borehole inclinometry ( Figure S12 in Supporting Information S1).
For a more realistic source model, we also applied a double couple source distributed over a fault of finite extent ("Okada model," Okada, 1992). However, because of the sparse coverage of the three dimensional displacement field measured with only two sensors, and the many degrees of freedom (e.g., fault size, fault orientation, location of rupture initiation), an ambiguity between the fault extension and the co-seismic slip couples both parameters making them unresolvable (Figures S9-S11 in Supporting Information S1).

Sub-Structure of the Stick-Slip Asperity
Our measurements close to the glacier bed reveal the substructure of a microseismic stick-slip asperity at the ice-bed interface. Weak asperities cluster downglacier of the main asperity's rupture initiation point, and produce events with approximately 60 times smaller moments than the main asperity. Assuming similar material properties and loading within our basal seismicity region, we conclude that the fault sizes of the weak asperities must be smaller than the main asperity. We calculate the average fault size of the weak asperities to 0.8 +1.8 −0.4 2 , assuming non-overlapping circular asperities with radii of half the distance to their closest neighbor. Accordingly, the total area covered by all 28 weak asperities sums up to ∼50 m 2 , that is, 60× the area of a single weak asperity. The main asperity's seismic moment is also ∼60× larger than the seismic moment of individual weak asperities. This linear dependence on fault area is expected for ruptures that include multiple subfaults (Madariaga, 1979, Supporting Information S1). Thus, we suggest that what we call the "main asperity" is a dynamic sequence of all weak asperities rupturing together as a multi-asperity rupture. These multi-asperity ruptures initiate upglacier of the weak asperities and propagate approximately downglacier by rupturing all weak asperities. This is in line with the moment gap in Figure 2b and agrees with Gräff and Walter (2021), who assumed that seismic moment variations at an asperity detectable at the ice surface are primarily modulated by stress drop and thus by the slip distance (Beeler et al., 2001;Deichmann, 2018, and Figure 3d). Both quantities are influenced by varying subglacial water pressures and the related loading at the bed.
Apart from the main asperity ruptures, weak asperities rupture mostly on their own and follow a Gutenberg-Richter law (Gutenberg & Richter, 1955) in which multi-asperity ruptures do not fit ( Figure S2 in Supporting Information S1). For super-critical asperity densities (Dublanchet, 2019), such multi-asperity ruptures that break through mechanical boundaries have been predicted and also observed at the Parkfield section of the San Andreas Fault: A typical "Parkfield earthquake" ruptures the fault section's microseismic asperities (Schorlemmer & Wiemer, 2005) to produce an M6 earthquake approximately once every 25 years (Tormann et al., 2012). This results in a gap in the magnitude distribution between these "characteristic" earthquakes and the microseismicity (Figure 1 in Jackson & Kagan, 2006). In this view, our "main asperity" is not a purely rate-weakening frictional patch surrounded by rate-strengthening material (T. Chen & Lapusta, 2009;Johnson & Nadeau, 2002), but a superposition of nearly simultaneous ruptures at many smaller asperities with a spatially constant location of rupture initiation. Some hypotheses explaining such behavior may be, for example, locally higher normal stress (Cattania & Segall, 2021) or regions that exhibit lower levels of fracture surface energy (Ide & Aochi, 2005;Noda et al., 2013). This is supported by one event at the end of the main asperity's active period that shows a magnitude 20 times weaker than typical events, indicating that it may have only ruptured a subset of weak asperities. Indeed, this would agree with the missing Doppler effect at the two accelerometers, that we attribute to a smaller rupture area ( Figure S5 in Supporting Information S1).

Seismic Versus Aseismic Glacier Sliding
The stress drop during a basal stick-slip event modulates the co-seismic slip. We calculate stress drops in the kPa-range for the multi-asperity ruptures. This is less or equal to 1%-10% of the 0-1.5 MPa inferred basal shear stress from gravitational loading (effective normal stress is unknown), and indicates that the asperity does not focus basal shear stresses from large surroundings, as the stress drops at the asperity would be orders of magnitude larger in such a case. Assuming that we investigated a typical glacier bed region hosting shear resistance, we conclude that the spatial separation of stick-slip asperities detectable from the surface of Rhonegletscher (Gräff & Walter, 2021) is much larger than the radius over which they influence basal sliding. This means that the effect of these microseismic stick-slip asperities on the total basal motion of Rhonegletscher must be small compared to aseismic sliding as described in traditional theories. However, we emphasize that Gräff and Walter (2021) omitted discrete basal seismicity with dominant isotropic radiation patterns, although they may be a manifestation of stick-slip sliding at the glacier bed with a stronger volumetric component caused by dynamic effects. Together with other forms of basal seismicity that were potentially overlooked so far, the influence of frictional sliding at the glacier bed may become relevant after all.
In contrast to Rhonegletscher, where the glacier bed mainly consists of hard bedrock with rocky humps and till-filled troughs, the situation on soft-bedded glaciers and ice-streams could be dramatically different. Kufner et al. (2021) found regions at the bed of Ruthford Ice Stream, Antarctica, that are covered with stickslip micro-seismicity. Their dense spatial distribution suggests that these stick-slip events may combine to crucially influence basal sliding in such a setting.

Summary and Conclusion
Regarding the questions posed in the introduction, our study provided the following new insights: (a) The concentration of meter-scale frictional asperities allows for their successive seismogenic rupture, eventually resulting in multi-asperity ruptures detectable at the glacier surface. (b) Rupture directivity suggests that the fault planes are located at bi-material interfaces, in agreement with ice resting locally on till. (c) Assuming a circular shape, the fault area of the asperity whose seismic ruptures are detectable at the surface is about multiple 10 square meters. (d) Stick-slip displacement over this asperity may reach millimeters for single events based on assumptions of elastic properties of the basal till that change with pore pressure. (e) For the investigated data, the stick-slip asperity is unlikely to affect basal shear stresses over significant distances.
The presented experiment, which included borehole instrumentation within a few meters of a frictional asperity, sets a new standard for monitoring microseismic stick-slip of glacier flow, especially for Alpine glaciers where conditions are less challenging than on polar ice streams. Event detection sensitivity and location precision are unprecedented, revealing hitherto unknown details of the structure and dynamics of a glacial stick-slip asperity, which behaves surprisingly similar to earthquakes on the Parkfield section of the San Andreas fault. Therefore, the responsible frictional mechanisms may be relevant over vast ranges of spatial and temporal scales. Relatively short recurrence times of glacial stick-slip, as well as hydraulic conditions and loading velocities varying over hourly to daily scales, provide more constraints to models of seismogenic sliding than earthquake observations once every few decades. We therefore suggest that glacier borehole seismology can provide a natural laboratory that could be used to elucidate various aspects of frictional sliding, including earthquakes.

Data Availability Statement
Seismometer data from near-surface stations RA81-RA88 and deep-borehole stations RA91-RA93 of the 4D local glacier seismology network (https://doi.org/10.12686/sed/networks/4d) are archived at the Swiss Seismological Service and can be requested via its web interface http://arclink.ethz.ch/webinterface/. Mini-SEED data of the IMS sensors, a catalog of the basal seismicity around the main asperity, and seismic station information are available in ETH Zurich's Research Collection (https://doi.org/10.3929/ethz-b-000511274).