Icequake Source Mechanisms for Studying Glacial Sliding

Improving our understanding of glacial sliding is crucial for constraining basal drag in ice dynamics models. We use icequakes, sudden releases of seismic energy as the ice slides over the bed, to provide geophysical observations that can be used to aid understanding of the physics of glacial sliding and constrain ice dynamics models. These icequakes are located at the bed of an alpine glacier in Switzerland and the Rutford Ice Stream, West Antarctica, two extremes of glacial settings and spatial scales. We investigate a number of possible icequake source mechanisms by performing full waveform inversions to constrain the fundamental physics and stress release during an icequake stick‐slip event. Results show that double‐couple mechanisms best describe the source for the events from both glacial settings and the icequakes originate at or very near the ice‐bed interface. We also present an exploratory method for attempting to measure the till shear modulus, if indirect reflected icequake radiation is observed. The results of this study increase our understanding of how icequakes are associated with basal drag while also providing the foundation for a method of remotely measuring bed shear strength.


Introduction
Understanding how glaciers slide over the underlying bed is an important process that is not yet fully understood. Glacial sliding is important because it is the dominant process controlling how solid ice moves off the land and into the oceans, contributing to sea-level rise (Ritz et al., 2015). However, "basal drag is a fundamental control on ice stream dynamics that remains poorly understood or constrained by observations" (Morlighem et al., 2010). Here, we use passive glacial seismicity observations, that is, icequakes, to study the basal drag of glaciers.
Icequakes are sudden releases of seismic energy due to the movement of ice. Icequakes originating at or near the bed of a glacier, associated with glacial sliding, can be used to investigate a number of physical properties and processes at or near the ice-bed interface (Podolskiy & Walter, 2016). Icequakes cannot completely elucidate glacier sliding processes, since ice flow is also accommodated aseismically through creep and viscous deformation. However, they do provide brief snapshots that provide insight into the physics of glacier sliding.
In this study, we use two icequakes associated with different glacial extremes to explore the following questions: (1) What icequake source mechanism fits the seismic data best? (2) To what extent can icequake source mechanisms be unified over two extremes of glacial settings and spatial scales? (3) Do the icequakes originate from the ice-bed interface, and if so, what can we learn about ice-bed mechanical coupling? (4) What fundamental properties of the bed can be remotely measured, such as the shear modulus of the till? (5) What are the fundamental limitations of using icequakes to investigate glacial sliding? The two particularly pertinent questions relevant for understanding basal drag better, and therefore the most significant results of our work, are as follows: how the ice is mechanically coupled to the bed, and whether it is possible to measure the shear modulus of the bed material.
The shear modulus of the till is an important parameter for ice dynamics modeling, since it is a measure of the elastic stiffness of the till. If slip of the ice is governed by failure at the ice-till interface or in the till, then the strength of the till controls the point of failure and therefore slip at the glacier bed. The shear modulus of the till is dependent upon till properties such as the density, porosity, and water content (Leeman et al., 2016). Measurements of the till shear modulus can therefore be used to obtain estimates of these till properties, which in combination with laboratory studies (Leeman et al., 2016;Tulaczyk et al., 2000) could be used to calculate till shear strength. Although such calculations are beyond the scope of this study, we present a novel method of remotely estimating the till shear modulus.
To explore these questions, we analyze icequakes from two glaciers that represent the extremes of different spatial scales (see Figure 1). The first location is an alpine glacier in the Swiss Alps and the second is an ice stream in West Antarctica. We present a detailed analysis of one icequake from each location. Each icequake is from a cluster of similar icequakes and so represents repeatedly observed behavior near the bed of each respective glacier. The icequake hypocenters are approximately at the ice-bed interface and are likely to represent the extremes of different glacial settings for which glacial sliding of ice over a bed occurs. While the icequakes analyzed here are thought to be representative of stick-slip seismicity at these locations, it is worth noting that we only present results for two icequakes, each only representative of a single cluster location geographically, and so these results should be treated primarily as exploratory findings that lay the foundations for implementation on larger data sets. Figure 1 shows the seismometer network geometries used to locate the icequakes and derive the most likely icequake source mechanisms. A source mechanism is a physical model of the most likely mode or modes of failure of a material subjected to an external stress, as well as the orientation of that failure. These source mechanisms, combined with their associated seismic radiation patterns and seismic moment of the energy released during failure, can be used to learn about the dynamic behavior of the slip of ice over the bed and the material properties of the surrounding media.
Icequakes originating at or near the ice-bed interface have previously been observed in glacial settings including: Antarctic outlet glaciers and ice streams (Anandakrishnan & Alley, 1994; 10.1029/2020JF005627 Bentley, 1993;Barcheck et al., 2018;Blankenship et al., 1987;Danesi et al., 2007;Smith, 2006, Smith et al., 2015Zoet et al., 2012); Greenland outlet glaciers (Roeoesli et al., 2016); and alpine glaciers (Allstadt & Malone, 2014;Dalban Canassy et al., 2016;Deichmann et al., 2000;Helmstetter et al., 2015;Walter et al., 2008Walter et al., , 2010Weaver & Malone, 1979). Much of this observed seismicity is interpreted to be associated with glacial sliding, specifically stick-slip behavior. Stick-slip seismicity occurs where patches of the bed, or ice-bed interface, are interpreted to have a higher shear strength, where basal drag is sufficient to inhibit flow until either the stress increases, or shear strength decreases, sufficiently to allow slip. Basal icequakes associated with tensile faulting have also been observed (e.g., Dalban Canassy et al., 2016;Walter et al., 2010). Although a significant number of studies have been undertaken on basal icequakes associated with glacial sliding, few have analyzed the icequake source mechanisms (Anandakrishnan & Bentley, 1993;Helmstetter et al., 2015;Roeoesli et al., 2016;Smith et al., 2015;Walter et al., 2010). To date, it has often been assumed that stick-slip seismicity should exhibit double-couple source mechanisms. This mechanism represents two coupled moment release pairs acting against one another to conserve angular momentum. One common example of this is when an earthquake is generated during slip between two tectonic plates. Here, we test this assumption by investigating all known types of fundamental earthquake source mechanisms, as well as two coupled mechanisms. The majority of previous studies have only inverted for first motion P wave polarities. Here, we perform source mechanism inversions using the full waveform for P, SV, and SH phases. This allows us to gain more information from the basal icequakes and allows us to explore the aforementioned questions in more detail than would otherwise be possible.

Methods
Source mechanisms for the two icequakes shown in Figure 1 are used to study the process of slip of ice over the bed. To derive the icequake source mechanisms, we constrain potential source models using the full waveform arrivals of P and S phases at seismometers near the glacier surface.

Data Processing
The icequake data presented in this study were collected by the networks shown in Figure 1. The network at Rhonegletscher, Switzerland, was comprised of three 3-component 1 Hz Lennartz borehole seismometers sampling at 500 Hz connected to Nanometrics Centaur digitizers and four 3-component 4.5 Hz geophones each connected to a Digos Data-Cube3 digitalizer sampling at 400 Hz. The Rhonegletscher data used in this study were collected in February 2018, corresponding to alpine winter conditions. The network at the Rutford Ice Stream, West Antarctica, was comprised of ten 3-component 4.5 Hz geophones connected to Reftek RT130 digitalizers sampling at 1000 Hz. These data were collected in January 2009, during the austral summer. The icequakes were detected using QuakeMigrate and a spectrum-based method, as discussed in Smith et al. (2015) and Hudson et al. (2019). This provides us with a catalog of icequakes from which we can select icequakes located near the glacier bed. Below we detail how specific icequakes are processed and why certain processing-related decisions are made.
In order to reduce the noise present for each phase arrival, we filter the data using the parameters shown in supporting information Table S1. We filter between 5 and 100 Hz for the Rhonegletscher icequake and 10 and 200 Hz using a four-corner causal Butterworth filter. Different filter parameters are used for the different glacial settings based on the different spectra of noise sources, the dominant source frequency of the basal icequakes, and the sampling rate of the data. The source of the higher frequency noise filtered out of the data could be due to natural sources such as surface winds or perhaps more likely instrument noise; hence the bandpass rather than highpass filter is applied. The icequakes' energy observed at receivers generally lies between 5 and 200 Hz. The phases are then separated, with the length of the waveforms passed to the full waveform inversion method specified in Table S1. Phases are rotated into the vertical (Z), radial (R), and transverse (T) components so as to approximately isolate the P, SV, and SH phases.
The icequakes are located by picking the P and S phase arrivals manually and then using the nonlinear location algorithm, NonLinLoc (Lomax & Virieux, 2000). Information regarding the phase picks are provided in Table S4. The ice velocity models used in the location procedure are given in Figure S1. The origin times and hypocentral locations are given in Table 1. In each case, the icequake depths correspond to the depth of the bed of the respective glacier found using ground penetrating radar (Church et al., 2018;King et al., 2016). Although the depth uncertainty is high, at ∼10% of the total icequake depth in both cases, this does not significantly affect the full waveform modeling, since the phase arrivals are manually aligned and the locations of the various layers and interfaces are all relative to the source location rather than the absolute geometry of the real glaciers.
Although we only analyze one icequake at each glacier in detail, these icequakes are representative of an entire cluster of icequakes observed at each location. Icequakes clustered both spatially and temporally are commonly observed at glacier beds and are thought to be caused by sticky spots, where the failure mechanisms are approximately identical when the sticky spot is seismically active (Roeoesli et al., 2016;Smith et al., 2015;Winberry et al., 2009). This is commonly referred to as stick-slip motion. The similarity of each icequake to its associated cluster is evidenced in Figure 2. The single Rhonegletscher icequake arrivals (red) and other icequakes in the associated cluster are shown in Figure 2a, and the single Rutford icequake and other icequakes in that associated cluster are shown in Figure 2b. In both cases the icequake that we study in detail is almost identical to all the other icequakes in the cluster. This repeatability is particularly remarkable for the Rutford icequake cluster. These observations provide us with confidence that the icequakes that we study here are representative of the behavior of basal icequakes at least for an individual cluster, and likely basal activity more generally, at each glacier. We are therefore confident that despite presenting the analysis of single events within this manuscript, the events used represent well the basal seismicity in that location.
Examples of the icequake arrivals at one station are shown in Figure 3a for the Rhonegletscher icequake and Figure 3b for the Rutford Ice Stream icequake. The seismograms for all the stations for each icequake can be The filters applied are specified in Table S1.  Table S1. Seismograms for all the stations for each icequake used in this study can be found in Figure S2. found in Figure S2. All P and S phase arrivals are clearly impulsive. The manually picked P and S arrivals are shown in red and blue, respectively. The P phase arrivals can clearly be seen on the vertical (Z) components and the S phase arrivals can be seen on the horizontal channels, as expected. The P-S delay times are much greater for the Rutford icequake because the source is ∼2 km below the glacier surface, compared to ∼200 m below the surface for the Rhonegletscher icequake. There are no surface wave phases observed, which in combination with the hypocentral locations gives us high confidence that the icequakes originate from near the glacier bed .
Significant shear wave splitting is observed in the Rutford Ice Stream icequake data, as observed in the same data set in Smith et al. (2017), probably because of the strongly anisotropic ice fabric (Harland et al., 2013;Smith et al., 2017) combined with ray paths of lengths greater than 2 km. We correct for this shear wave splitting using the method of Wuestefeld et al. (2010), implemented using the software SHEBA. This is based on rotating and shifting the seismograms in time (Silver & Chan, 1991) to find the most robust solution. SHEBA also implements the multiwindow clustering analysis method of Teanby et al. (2004) to minimize the impact of the choice of S-wave window used in the automated shear wave splitting analysis (Wuestefeld et al., 2010). The parameters found by this method and applied to the data to remove the splitting effects are given in Table S3.

Full Waveform Source Mechanism Inversion
The icequake source mechanisms presented in this study are found by using a Bayesian inversion method similar to that detailed in Pugh et al. (2016), but instead using the full waveform of various phases. We use a Monte Carlo based technique to randomly sample potential source models, ensuring no bias within the n-dimensional space (where n is the number of dimensions of the source model). For such a source model, we can calculate the observed displacement, u n , at a seismometer (Walter et al., 2009), where n denotes a particular seismometer, M is a vector composed of the source model parameters, for example, of length six for a full moment tensor model, and G n (⃗ x, t) is a two-dimensional matrix containing the Green's functions associated with each model component. The Green's functions account for path effects due to the medium.
We investigate the following source mechanisms in this study: a Double-Couple (DC) source mechanism (three free parameters); a Single-Force (SF) source mechanism (three free parameters); an unconstrained Moment Tensor (MT) source mechanism (six free parameters); a DC-crack coupled mechanism (seven free parameters); and a SF-crack coupled mechanism (seven free parameters). Examples of the physical manifestation of these source mechanisms are shown in Figure 4. First motion radiation patterns for each source model are shown, to indicate an instantaneous component of the overall waveform that we simulate. The DC and MT models implicitly suggest that away from the source, the ice is mechanically coupled to the bed, while the SF sources suggest that the ice and bed are mechanically decoupled away from the source (Dahlen, 1993). We use the term mechanically coupled to refer to regions distal to the fault behaving such that the ice-bedrock interface is static with no slip occurring. This latter source is typically used to describe landslide source mechanisms (Allstadt, 2013;Dahlen, 1993;Kawakatsu, 1989). A single-force source suggests mechanical decoupling of the ice from the bed because it describes one body accelerating over another, which can only occur if the two bodies are decoupled. This is in contrast to the DC and MT models, where even at a bimaterial interface, the moment release is constrained to a finite length fault plane and the moment tensor only describes deformation at the source (Vavryčuk, 2013). Beyond the finite spatial limits of Source-time function the source, the material is required to be mechanically coupled, even for a bimaterial interface, for example, in the model presented in Shi and Ben-Zion (2006).
The Green's functions used in Equation 1 are generated using the software fk (Haskell, 1964;Wang & Herrmann, 1980;Zhu & Rivera, 2002). The program takes a one-dimensional layered velocity model, a source-time function, and the epicientral distance and azimuth of receivers from the source, with the parameters used for each icequake case given in Table 2. We do not invert for the source-time function but used a fixed time duration as specified in Table 2.
The displacement at a seismometer can be calculated from Equation 1, once the Green's functions have been generated for a particular randomly sampled source mechanism model. This modeled displacement can then be compared to the real, observed displacement. There are a number of methods for quantifying the misfit. We use the variance reduction method (Templeton & Dreger, 2006;Walter et al., 2009), where the variance reduction value is given by, where is the misfit, v n, data (t) is the observed velocity at seismometer n over time, and v n, model (t) is the modeled velocity for seismometer n over time, calculated by differentiating Equation 1 with respect to time. The probability of the data fitting the model, P(data|model), assuming Gaussian statistics, is then defined by the likelihood function,  (Bodin & Sambridge, 2009), The probability of the randomly sampled source mechanism model fitting the data can then be found using Bayes' theorem (Bayes & Price, 1763), where the posterior probability, P(model|data), is given by, All the sampled models are assumed to have identical initial prior probabilities, therefore P(model) is given by, where N is the number of samples used in the inversion, typically 1 × 10 6 . Evidence that this is sufficient is provided in Figure S3, which shows both the equal area sampling of the spatial orientation of source mechanism and equal angle sampling of the full moment tensor space. These distributions are representative of the spatial sampling for all source model types. However, obtaining P(data) is more challenging. We find Figure 5. Schematic diagram and example synthetic and observed data demonstrating the method used to estimate till properties in this study. The direct waves travel directly from the source to the receiver (gold triangle), passing through ice only. The indirect waves travel downward first, reflecting off a lower interface below a till layer, before traveling up toward the receiver. On the right are the Z, R, and T components predicted for a near surface seismometer offset laterally by ∼90 m from the source for the Rhonegletscher icequake. The waveforms show the direct, indirect, combined (75% direct, 25% indirect) synthetic energy arriving at the seismometer, as well as the observed energy. The till layer thickness associated with this inversion is 1 m. P(data) by using Bayesian marginalization (Tarantola & Valette, 1982), where P(data) can then be defined by, Using a Monte Carlo based approach to sample a large number of models, typically of the order of 10 6 , provides us with an estimation of the full posterior probability distribution (pdf) for a particular type of source mechanism model. The most likely source mechanism model can then be found, along with an estimate of its associated uncertainty, taken to be the standard deviation of the pdf.
The different source mechanism models shown in Figure 4 have different numbers of free parameters. In order to account for the complexity of a particular model when comparing the various model types, we use the Bayesian Information Criterion (BIC) (Schwarz, 1978). The BIC allows us to assess whether a model with more free parameters overfits the data relative to one with fewer parameters. It is given by, where k is the number of free parameters for the model and n is the number of samples, or data points, used in the inversion. The difference in BIC value between two models i and j, ΔBIC i, j , can be used to compare the relative fit of one model to the other. If ΔBIC i, j < 3.2, then there is insufficient evidence to suggest that model i is better than model j, whereas if 3.2 < ΔBIC i, j < 10 then there is substantial evidence to suggest that model i is more appropriate than model j, and if ΔBIC i, j > 10, then there is strong evidence that model i should be favored over model j (Kass & Raftery, 1995).
The full waveform inversion method described allows us to find both the most likely model for a specific type of source mechanism and to intercompare different types of source mechanism, while rigorously accounting for uncertainty in the results.

Subglacial Till Properties From Full Waveform Icequake Source Mechanism Inversions
If an icequake source has both direct and indirect arrivals, that is arrivals traveling straight from the source to the receiver and arrivals that have reflected off or refracted at some interface below the source, respectively (see Figure 5), then one can learn something about the medium beneath the icequake source. If the icequake is located at the ice-till interface, with a reflective till-bedrock interface below the till, then one can approximately measure the bulk and shear moduli, K till and till , of the till. A full derivation of this method can be found in the Text S5, with a summary provided here.

10.1029/2020JF005627
The observed velocity, v obs, i (t), at the receiver is given by, where i denotes the seismic phase (P or S). v dir (t) is associated with the energy propagating directly from the source to the receiver (see Figure 5). v indir, i (t) is associated with energy that is radiated downward, before reflecting off an interface that we define as the till-bed interface. For our model scenario, we approximate this path using a source within ice, at a variable distance vertically above a bedrock interface, with this distance representing the simulated till thickness. We do this because the method for generating the Green's functions that we use, fk (Haskell, 1964;Wang & Herrmann, 1980;Zhu & Rivera, 2002), does not allow us to place a source at an interface between two media. We therefore approximate a source at an interface between ice and till by separating direct and indirect arrivals using a homogeneous ice velocity model and an ice with a granite bed velocity model.  i is defined as the additional proportion of indirect waves observed at the receiver, ranging from 0 to 1.
An example of v dir (t) and v indir (t), the time derivatives of displacement, are shown in Figure 5. The combined modeled velocity and real observed velocity at the example seismometer are also plotted. One can see from the waveforms in Figure 5 that the differences between different arrivals are subtle, with small changes in relative amplitudes between the different modeled phases. It is therefore necessary to have sufficiently high resolution observations in order to perform this analysis.
Theoretically, the value  i is defined by, where R till − bed, i is the reflectivity coefficient for seismic phase i at the till-bed interface and T till − ice, i is the transmissivity coefficient for seismic phase i at the till-ice interface. If we make the assumptions (1) that the radiation is approximately at normal incidence to each bimaterial interface and (2) that any P to S and S to P conversions are approximately compensated for with one another, then from the Zoeppritz equations (Aki & Richards, 2002;Zoeppritz, 1919) we can obtain the following simplified relations for  P and  S , where Z p, ice, till, bed and Z s, ice, till, bed are the P and S phase impedance for the ice, till, and bed. Z p, ice , Z p, bed , Z s, ice , and Z s, bed are known, or can at least be approximately assumed. If we can obtain values of  P and  S then we can use these equations to solve for Z p, till and Z s, till , which in turn can be used to give us the bulk and shear moduli, K till and till , of the till in terms of the density of the till, till (Dvorkin et al., 1999), To solve Equations 12 and 13 to find K till and till , we therefore need to obtain values for  P and  S . This can be done by performing a full waveform source mechanism inversion as described previously but also inverting for the proportion of indirect P and S waves observed at receivers approximately directly above the source. To perform this inversion, we rewrite Equation 8 as, where v homo ice, i (t) is the modeled velocity signal for a medium comprised of only ice without material interfaces and v rock bed, i (t) is the modeled velocity signal for a medium with a faster velocity reflecting bed at a depth below the source equal to the thickness of the till layer.

10.1029/2020JF005627
Since we have models to calculate the velocity signals v homo ice, i (t) and v rock bed, i (t), we can therefore perform the full waveform inversion with an additional two parameters,  P and  S , which we vary randomly and uniformly between 0 and 1. The best fitting model can then be used to determine the best values for  P and  S . From this we can then calculate K till and till from Equations 12 and 13.
One assumption we make is that at the glacier bed, there are three layers with distinct velocity contrasts: an ice layer; overlying a till layer; overlying a bedrock layer. A justification of this assumption is given in section 3.1. A further important assumption we make is that the indirect radiation from an icequake (see Figure 5) is approximately at normal incidence to the ice-till and till-bed interfaces, in order to simplify the Zoeppritz equations (Zoeppritz, 1919). To make this assumption, in the till properties inversion we only use stations close to the icequake epicenter, with maximum angles less than 24 • from normal incidence. At an angle of incidence of 24 • , the reflectivity coefficients at the interfaces are predicted to vary from approximately 10% to 25% for P waves, depending upon the materials comprising the interface (Booth et al., 2016).
These are approximately accounted for at the reflecting interface, albeit for an ice-bedrock interface rather than a till-bedrock interface. Ideally one would also account for such variation at the transmitting interface between ice and till, although for simplicity we do not correct for angle of incidence effects at that interface in this exploratory study. A final note of relevance to our method is that we do not have to account for thin bed effects (Widess, 1973) since we are simulating the source at the upper interface of the thin bed, so there is no upper reflection that would otherwise interfere with reflections off the lower interface of the thin bed. In any case, a strength of our general full waveform source mechanism inversions undertaken in our entire study is that we generate all reflections in our modeled seismic waveforms, and so account for thin bed effects in our other inversion results presented in sections 3.1 and 3.2.

Variation of Icequake Source Depth, Source Mechanisms, and Bed Structure
Icequake source depth, source mechanism, and bed structure are varied for each glacial setting. The results are plotted for Rhonegletscher in Figure 6a and the Rutford Ice Stream in Figure 6b. Each point on the plots indicates the most likely result of one full waveform source mechanism inversion. The composition of the various bed structures with depth are shown Figure 6, below their associated plots. Both glacial settings generally indicate that the closer the source is to the ice-bed interface, the higher the similarity value and therefore the better the model fits the data. In the Rhonegletscher case, the highest likelihood model is for ice with bedrock but no overlying till layer. In the Rutford Ice Stream case, the highest likelihood model result is for an ice half space with no bed. The highest likelihood models are invariably those of greater complexity, with more free parameters (the full MT, DC-crack, and single-force-crack models).
The highest likelihood and therefore best fitting model for the Rhonegletscher icequake is a singleforce-crack source mechanism, with the icequake ∼5 m above an ice-rock interface. The best fitting model for the Rutford Ice Stream icequake is a full moment tensor source mechanism with no apparent bed below the source. However, these models have more free parameters than the DC or single-force models. Accounting for this additional complexity when comparing different types of source mechanism model can be undertaken by using the Bayesian information criterion (see Equation 7). Table 3 gives the ΔBIC complex − simple values for Rhonegletscher and the Rutford Ice Stream, with the high, positive values (>9) in Table 3 suggesting that the simpler, DC or single-force source model should be favored over the more complex models, for the highest likelihood models given in Figure 6. After accounting for this complexity, the most likely source mechanism is either the DC or single-force source mechanisms for the Rhonegletscher icequake and is the DC source mechanism with an ice-only half space for the Rutford icequake. Although the single-force mechanism for the Rhonegletscher icequake has a slightly higher similarity value at 0.43 as opposed to the DC similarity value of 0.42 (see Table 3), there is no statistically significant difference between the two, with ΔBIC DC − SF ≈ 0. We therefore cannot make a distinction between the two. However, the SF source provides a much poorer fit than the DC source for the simpler homogeneous ice velocity model for the Rhonegletcher icequake. We therefore infer that the DC source provides a universally better fit overall, and so we present the DC source model as the best overall fit for both the Rhonegletscher and Rutford icequakes. These source mechanisms are discussed in more detail in section 3.2.
One potential source of bias associated with the results in Figure 6 is polarity reversal of the P and S waves as the source depth is varied, with polarity reversals occurring at half the dominant wavelength of the relevant  Roethlisberger (1972), bedrock velocity (taken to be that of granite) from Walter et al. (2010), and till velocity is based on Antarctic till . Note that since the ice only case has no interfaces at depth, the inversion is performed for one depth only (0.2005 km below surface for Rhonegletscher, 2.0375 km below surface for Rutford) and extrapolated for comparison with the other inversion results.
phase. Such a polarity reversal could cause an anticorrelation between the modeled and observed signal, potentially resulting in a misleadingly low similarity value. The length scale over which the polarity of a phase would reverse is the order of 12 to 18 m for the P waves we observe and 24 to 36 m for the S waves. However, we manually align the P phase arrivals of the modeled greens functions with the observed seismic signals, and check that the first arrival polarities are consistent. This minimizes any polarity reversal bias for the P wave, but theoretically the S wave could still observe polarity reversals that are not compensated for. However, the peaks in the similarity values near the ice-bed interface are significantly narrower than the depth difference over which a P or S wave could reverse polarity, being approximately 1 to 4 m wide. We are therefore confident that our findings in Figure 6 are not biased by P and S phase polarity reversal caused by varying source depth.
The best fitting velocity models, the ice-only velocity model for Rutford and the ice-rock velocity model for Rhonegletscher, both indicate that either there is no till layer present at the glacier bed or that the seismic signals do not exhibit reflections off an ice-till impedance contrast. From experiments drilling to the bed (Gräff & Walter, 2019) and seismic studies, at alpine and Antarctic glaciers (Iken et al., 1996;Smith, 1997a;Smith & Murray, 2009), it is likely that there is a till layer present at the bed of both Rhonegletscher and the Rutford Ice Stream. This leads us to the interpretation either: that the icequake source is at the ice-till interface, therefore resulting in no reflections off a till layer or that the ice-till interface is poorly defined, with a very gradual change in seismic velocity gradient, again resulting in no reflections. The latter interpretation is deemed less likely given the length scales for such a gradual change in velocity constrained by the inversion results (<1 m). We therefore suggest that it is most likely that the icequakes originate at the ice-till interface. This interpretation agrees with the findings presented in section 3.3, since it allows for there to be a till layer present, as assumed in the results in section 3.3 and consistent with active source observations at Rutford Ice Stream (Smith, 1997b), while not requiring any ice-till reflections to be observed. Note. VR is the data-model variance reduction value, a measure of the quality of the fit. ΔBIC complex − simple is the difference between the highest likelihood complex and simple icequake source mechanism solutions. Details of how the seismic moments are calculated can be found in Text S6 and are based upon and elaborated upon in Aki and Richards (2002), Shearer (2009), Peters et al. (2012, and Hudson (2019). The half space that gives the highest variance reduction value is shown in brackets (e.g., ice-the ice only half space; gb-the ice with a granite bed half space).
It is worth noting that although we suggest that the icequake source is most likely at the ice-till interface, this does not necessarily imply that on the scale of the fault slip, the fault interface is that of either ice-till or ice-bedrock. It may be that at this scale, the seismicity is induced by contact between small rocks or other entrained sediment that are frozen into the glacier ice at the ice-bed boundary (Lipovsky et al., 2019).

Best Fitting Icequake Source Mechanisms
The best fitting source mechanisms from all the full waveform inversion results are shown in Figure 7a for the Rhonegletscher icequake and in Figure 7b for the Rutford Ice Stream icequake. Due to the depth of the Rutford icequake source (∼2 km below the surface), the P-S delay time for the Rutford icequake is sufficiently large that we split the P and S arrivals, with the P phase fitted on the Z component, and the S phase fitted on the R and T components. The apparent negative time offset of the S arrival relative to the P arrival in the observations in Figure 7b is therefore simply a result of where the waveforms are cut for each phase, with the Z component and the R and T components not aligned temporally with one another. All the modeled (red waveforms) phase polarities for P, SH, and SV phases are in agreement with the observed (black waveforms) polarities for both icequakes. Furthermore, the various modeled phase amplitudes are also in generally close agreement with the observed amplitudes, with significant later phase arrival complexity captured by the best source mechanism models for certain stations.
Since the simplest best fitting source mechanisms are DC, the slip vectors can be calculated, the directions of which are shown by the red arrows in Figure 7. An estimate of the uncertainties associated with each slip vector are shown by the red dashed lines. The slip vectors both approximately agree with the ice flow direction at each location (see Figure 1). While this might be expected, it should not be assumed as the ice slip direction associated with a single icequake is not required to match the overall slip direction of a glacier (Anandakrishnan & Bentley, 1993). Therefore, while our observed slip vectors are in agreement with the overall direction of glacial motion, all that can be interpreted from this result is that the icequake The observed waveforms at each seismometer are shown in black, for the Z, R, and T components, along with the modeled waveforms, shown by the red dashed waveforms. Note that the waveforms for the Z component for the Rutford data in (b) are not temporally aligned with the R and T components, for reasons given in the main text. The complete seismograms for each event can be found in Figure S2.
likely accommodates some of the overall motion of a glacier. A more significant result is that the vertical orientation of one of the fault planes for each icequake, and its associated slip vector, indicates slip along a sub-horizontal bed. The Rhonegletscher fault inclination is greater than the Rutford Ice Stream fault inclination, which is indeed the case in reality, as the alpine glacier has steeper bed topography than the Antarctic ice stream.
One potential issue with inverting for source mechanisms using our method is that the Green's functions used are effectively only generated in 1-D (Zhu & Rivera, 2002). In reality, 3-D source effects that could be caused by sudden local variations in bed topography, the presence of eroded material, basal crevassing introducing a local anisotropic ice structure, or accumulation of melt water (Walter et al., 2010), could have a detrimental impact upon our results. Indeed, 3-D source effects are shown to be important for near-bed tensile crack source mechanisms at Gornergletscher, another Swiss alpine glacier (Walter et al., 2010).  Figure 7. The complete seismograms for each event can be found in Figure S2.
To test whether 3-D effects affect our results, we perform the same inversions as used to obtain the results in Figure 7, but for the SH phase only (i.e., using the T component only). The SH phase is far less insensitive to 3-D effects for the geometry of our scenario, as it is parallel to the fault. If the SH inversions give similar results to the inversion using all body wave phases then one can assume that 3-D effects are of second order in our case. Such results are shown in Figure 8. It can be seen that these SH phase inversions are similar to the inversion results that use all body wave phases in Figure 7. The similarity of 3-D dependent (the P, SV, and SH joint phase inversions) and the 3-D independent SH phase inversions implies that our results are insensitive to any local topography, ice fabric damage, and the potential presence of meltwater. Our results are therefore robust and not substantially affected by 3-D effects.
A further possible source of uncertainty in the source mechanism inversion results could be caused by the vanishing traction condition at the free surface. If an earthquake source is sufficiently shallow, then the M xz and M yz terms of the moment tensor can approach zero and effectively vanish. If this is not accounted for when inverting for a shallow earthquake, then it can result in an inversion bias, for example, making shallow DC sources appear to a vertical dip-slip mechanism (Chiang et al., 2016) similar to the mechanisms we observe. However, any vanishing traction effects manifesting themselves in our results would be minimal, since although the icequakes are shallow in comparison to tectonic earthquakes, the source wavelengths are much shorter than the icequake depths below the surface, therefore not observing the vanishing traction effect (Chiang et al., 2016). For example, if one assumes a conservatively low dominant source frequency of 100 Hz for the Rhonegletscher icequake at a depth of 200 m below surface, the wavelength would be ≈36 m, which is much less than the source depth.
Assuming that the icequake source is located approximately at the ice-bed interface, the DC nature of the best fitting source mechanisms implies that the ice is mechanically coupled to the bed distally from the source. This is in contrast to the case where the best fitting source mechanism is a single-force source mechanism, where the overlying fault block (ice) would be free to accelerate relative to the underlying block (till or bedrock). Such a single-force mechanism would imply that the ice would be free to slide over the bed, mechanically decoupled from the bed. The significance of the DC source mechanism observation, and hence the implied mechanical coupling distally from the source, is that the net movement of the entire glacier is not attributed to a single micro-icequake. This has implications for how to understand glacial sliding on the spatial scale of an entire glacier. Here, we assume that this observation of distal mechanical coupling of the ice to the bed is either due to freezing of the ice to the bed (Christoffersen & Tulaczyk, 2003;Joughin et al., 2004) or due to a sufficiently high coefficient of friction at the ice-bed interface that is likely modulated by fluids (Tulaczyk et al., 2000). Figure 9. Plot of till shear modulus (μ) against density for the Rhonegletscher icequake. Black line is the inversion result, with the gray shaded region indicating the uncertainty. Scatter points show the shear modulus associated the ice and bedrock used in this study (Podolskiy & Walter, 2016;Walter et al., 2010), as well as values of till modulus calculated from Amplitude Vs. Offset (AVO) seismic observations for the Bindschadler Ice Stream, Antarctica (Peters et al., 2007), and laboratory derived measurements of till shear modulus from: Whillans Ice Stream, Antarctica (Dvorkin et al., 1999;Leeman et al., 2016); Storglaciaren, Sweden (Iverson et al., 1999); and the Laurentide paleo ice sheet (Rathbun et al., 2008). The uncertainties associated with the AVO observations are plotted as colored lines.

Investigating Till Properties
One of the most useful observations for constraining glacial sliding within numerical models is the strength of the interface between the ice and the bed, since this parameter governs the conditions under which ice will slide. If the bed is stiff and cannot deform then this bed strength is the frictional force per unit area of the interface. If the bed can deform then the strength of the interface is governed by the shear strength of the bed. Laboratory studies of till strength have been undertaken (e.g., Leeman et al., 2016;Tulaczyk et al., 2000). Since we have some confidence from previous studies that there is at least a thin till layer (of the order 10s cm to meters at Rhonegletscher) where our icequakes originate, in this section we attempt to estimate the till shear modulus using our icequake observations. As previously mentioned, the till shear modulus is an important parameter because it contains information regarding till properties that are required for estimating the shear strength of the till or ice-till interface.
The method we use to estimate the till shear modulus in this exploratory study is outlined in section 2.3. Unlike normal incidence active source seismic surveys, it is possible to obtain estimates for the till shear modulus since we have P and S phases with which to constrain our inversion results. The difference between the previously discussed results and the approach taken for the results in this section is primarily that we invert for the additional parameters  P and  S , the proportion of indirect P and S waves observed in addition to the direct phase arrivals. These values can then be used to derive the relationship between till density and till shear modulus, as described in section 2.3.
The results of the till properties inversion for the Rhonegletscher icequake are plotted in Figure 9 and given in Table 3. The source mechanism associated with the inversion is plotted in Figure S4. We do not invert for till thickness for the Rhonegletscher icequake, with the till layer being fixed at a thickness of 2 m, due to the required computational expense. Varying the till layer in an inversion scheme would change the waveform shape, as well as amplitude, and would accommodate potentially thicker till layers that are observed elsewhere Luthra et al., 2016). Table 3 shows that the variance reduction for the DC source mechanism when also inverting for till properties provides a better fit than the DC source mechanism found in the previous sections of this study, with VR DC, TPI equal to 0.51 compared to a VR DC value of 0.42. This implies that modeling for direct and indirect arrivals using our till properties inversion method is valid. The shear modulus is plotted against till density, with the till density range specified based on geophysical, field, and laboratory measurements (Halberstadt et al., 2018;Hausmann et al., 2007;Iverson & Iverson, 2001;Peters et al., 2007Peters et al., , 2008Truffer et al., 2000). Considering the assumptions made and the associated uncertainty of the till shear modulus result (see Figure 9), the shear modulus is similar to that predicted by the laboratory experiment results plotted in Figure 9 (Iverson et al., 1999;Leeman et al., 2016;Rathbun et al., 2008), with all the measured till shear moduli results except one falling within the uncertainty bounds of 10.1029/2020JF005627 our results. The outlier is the stiff till at the Bindschadler Ice Stream (Peters et al., 2007), which is not concerning as it simply implies that the till in our study is more likely soft than stiff. However, it is clear from Figure 9 that the uncertainty magnitude is significant. It is worth noting that the lower till shear modulus we find compared to that found for the Whillans Ice Stream could be a result of the lower effective normal stress at an alpine glacier compared to an Antarctic ice stream due to thinner ice (up to ∼16 MPa in our case, excluding basal water pressure effects), or because in situ till has a lower stiffness than that suggested by laboratory experiments (Winberry et al., 2009).
A limitation of the till inversion results is the spatial resolution of till layer thickness that one can resolve using an icequake. The spatial resolution is related to the spectrum of the icequake source and the observed spectrum at the receiver. The highest frequency component of the source spectrum provides a fundamental limit on the spatial scale that can be resolved by our till properties inversion method. In our study, our modeled source time function has a duration of 1 × 10 −3 s or less, potentially allowing for a till thickness sensitivity of 3.6 m for P waves and 1.8 m for S waves. The real source time function for an icequake could be of an even shorter duration than we assume in this study. However, we cannot resolve such high frequencies at the surface: partly due to attenuation in the medium; and partly due to the sampling rate and data processing, such as bandpass filtering to remove noise. Figure 10 presents a limited analysis of the ability to resolve a till layer 2 m thick when attenuation and receiver filtering for the Rhonegletcher icequake. The waveforms in Figure 10a show the observed waveforms arriving at reciever RA54 and the waveforms in Figures 10b-10d are for a modeled source with various different filters applied. For no attenuation of the medium and no filtering at the receiver, in Figure 10b, one can observe the fill complexity in the various arrivals due to the presence of the 2 m thick till layer. When realistic attenuation is introduced in Figure 10c, some of the complexity in the various arrivals is preserved, but some of the smaller amplitude, higher frequency components are lost. When realistic attenuation and signal filtering are applied at the receiver, Figure 10d, further complexity and more of the higher frequency features are lost. The latter data in Figure 10d are comparable to that used in our till properties inversion and that of the observed waveforms in Figure 10a. We therefore conclude that some critical information is lost through the natural attenuation characteristics of the ice medium but also due to the frequency response of the instruments rather than subsequent data processing. However, there is still some remaining phase information that is incorporated into the till properties inversion. We do not perform the till properties inversion with no filtering, since the noise conditions would result in potentially spurious inversion results, and in any case the dominant filtering is likely caused by the instrument response rather than our data processing. Unfortunately, the method we present here is therefore significantly limited by frequency filtering of the signal, and also to some extent by attenuation of the medium, so the results should be interpreted tentatively. One could attempt to remove the requirement for filtering by either finding events with less background noise present or by deploying instruments deeper into the ice, where the noise conditions are likely lower. Furthermore, the instrument sampling rates could be increased, increasing the Nyquist frequency limitations of the data. This would only be of benefit if the instrument response was sufficient at higher frequencies. One could also attempt to better understand the attenuation structure, again reducing the uncertainty associated with generating the Green's functions. Similarly, one could attempt to understand the source-time function characteristics better, possibly even inverting for the source-time function. Such an understanding of the source-time function would inform us of the maximum theoretical till thickness that we could resolve using a passive icequake source.
We also tried to invert for till properties for the Rutford Ice Stream icequake, since we are confident that there is also a till layer present where the icequake originates. Such an inversion would be expected to work if there is a strong seismic velocity contrast between the till and underlying bedrock, leading to strong reflections, like those observed at Rhonegletscher. However, we could not obtain realistic estimates for the bulk and shear moduli using our method, even when varying our till layer thickness over a range of 1 to 40 m. This differs from our previous experiments where the till layer thickness was fixed at 2 m (see Figure 6). The failure to obtain a realistic result from the inversion is likely to be because the seismic velocity contrast between the till and the bedrock is less distinct at this point on the bed of the Rutford Ice Stream than at the Rhonegletscher bed, with the bedrock at the Rutford Ice Stream being sedimentary (Smith, 1997a;Smith & Murray, 2009) compared to the higher seismic velocity bedrock in the Alps (Walter et al., 2010). If there is an insufficient impedance contrast, then any reflected energy will be weak and attenuated before reaching the surface, resulting in a null inversion result. This is a limitation of our method. However, if the seismic wave field were sampled at a higher spatial resolution, and/or a larger magnitude icequakes were observed, it may be possible to overcome this limitation.
Although we use passive seismic observations, which act as a P and S wave source, active seismology methods using a P wave source only can also be used to investigate glacier bed properties. The most useful active seismic method is amplitude-variation-with-offset/angle (AVO/AVA). This method involves using a near surface active source to generate P waves that then reflect off the ice-bedrock or ice-till interface and are measured at a number of surface receivers. If the source-receiver offset is varied, then the observed incidence angle of the P wave reflecting off the bed is varied. The reflectivity coefficient varies with P wave incidence angle and observational results can be compared to theoretical predictions in order to infer bed properties (Booth et al., 2016). AVO has been undertaken on glaciers and can be used to infer till properties such as whether the till is consolidated or unconsolidated and whether the bed is wet or dry (e.g., Christianson et al., 2014;Peters et al., 2007Peters et al., , 2008. We have plotted the till shear modulus calculated for AVO observations of S wave velocity and till density at the Bindschadler Ice Stream, Antarctica (Peters et al., 2007), in Figure 9. While the soft till result is in agreement with our observations, the uncertainties associated with AVO measurements are typically much smaller than using the passive seismic method outlined in this study. However, AVO studies are limited by the incidence angle that can be observed, important for deciphering between different bed models (Booth et al., 2016). Such studies are also limited by the practical challenges involved with the survey setup. For measuring incident angles of up to 40 • for ice 2 km thick, one would require a source-receiver spacing of over 3 km, with many receivers in between to provide adequate variation of incident angle. Such a survey would only provide a single point measurement at a certain location for one instant in time. Obtaining multiple measurements in a field campaign therefore is challenging. Theoretically, using passive seismic sources with the method we propose allows for a measurements of till properties at various locations within a seismic network over a period of time, as long as the icequake sources vary spatially and temporally. Our method could therefore complement active seismic methods for providing improved measurements of glacial bed conditions.
To our knowledge the Rhonegletscher till shear modulus result is the first remotely estimated value of shear modulus using passive observations, an important observational parameter for constraining the rheological properties of the till for ice dynamics models. The failure of our method to obtain a till shear modulus result for the Rutford Ice Stream further emphasizes that our method requires further development and greater sampling of the seismic wavefield, or larger magnitude icequakes. Nevertheless, our method of obtaining till shear modulus provides a valuable foundation for making observations of basal shear strength at glaciers. Figure 11 summarizes the key findings of this work. Firstly, a DC mechanism provides the best fit to the observations. Although this has been assumed in previous stick-slip icequake studies, we show that this is the best source mechanism model for such basal icequakes, using information from the full waveforms to constrain the results. Secondly, the icequake source mechanism appears to be independent of glacial scale, with an alpine stick-slip icequake at 200 m depth exhibiting the same properties as an icequake from a 2 km thick Antarctic ice stream. This result suggests that alpine icequakes could be used for studying basal sliding of bodies of ice at any scale. The significance of this result is that it is often far easier to access and observe phenomena on alpine glaciers due to their comparatively easy accessibility and the thinner layer of ice between the surface and the bed. Thirdly, stick-slip icequakes most likely originate at, or very near (<1 m), the ice-bed interface. The best fitting source mechanism results indicate that failure of the system during a sliding event most likely occurs at the ice-bed interface, with the waveforms being relatively simple suggesting few reflections off interfaces. The fourth result of this study is that our full waveform source mechanism results are approximately independent of 3-D effects, to first order. The fifth result is that the bed is coupled to the ice distally from the source. This result agrees with the theory that the bed has patches of higher friction, that is, it is mechanically coupled in multiple locations (e.g., Alley, 1993). The final result is that in certain circumstances it may be possible to use an icequake remotely estimate the till shear modulus, allowing for the possibility of constraining how ice dynamics models simulate basal sliding using real, remotely measured basal sliding observations with meaningful measured parameters. Although we only show this tentative observation for the alpine icequake, our method provides a foundation for future studies, where better constraint of the till shear modulus might be possible.

Data Availability Statement
The Rhonegletscher data for stations RA51-57, part of the 4-D local glacier seismology network (https://doi. org/10.12686/sed/networks/4d/) are archived at the Swiss Seismological Service and can be accessed via its web interface http://arclink.ethz.ch/webinterface/. The Rutford Ice Stream data are available through the IRIS data center, under network code YG (2009), Gauging Rutford Ice Stream Transients. The full waveform source mechanism inversion code used in this study is made available from Hudson (2020). Tom Hudson was funded by the Cambridge Earth System Science NERC Doctoral Training Partnership.