Highly variable friction and slip observed at Antarctic ice stream bed

The slip of glaciers over the underlying bed is the dominant mechanism governing the migration of ice from land into the oceans, with accelerating slip contributing to sea-level rise. Yet glacier slip remains poorly understood, and observational constraints are sparse. Here we use passive seismic observations to measure both frictional shear stress and slip at the bed of the Rutford Ice Stream in Antarctica using 100,000 repetitive stick-slip icequakes. We find that basal shear stresses and slip rates vary from 104 to 107 Pa and 0.2 to 1.5 m per day, respectively. Friction and slip vary temporally over the order of hours, and spatially over 10s of metres, due to corresponding variations in effective normal stress and ice–bed interface material. Our findings suggest that the bed is substantially more complex than currently assumed in ice stream models and that basal effective normal stresses may be significantly higher than previously thought. Our observations can provide constraints on the basal boundary conditions for ice-dynamics models. This is critical for constraining the primary contribution of ice mass loss in Antarctica and hence for reducing uncertainty in sea-level rise projections. Passive seismic observations from the Rutford Ice Stream in Antarctica reveal a highly complex bed and substantial variability in friction and slip rates at the ice–bed interface.

Glacier slip is the primary mechanism governing the migration of ice from land into the oceans, with accelerating slip providing a major contribution to sea-level rise 1,2 . Friction at the bed of a glacier fundamentally limits the speed at which the ice can slip. This friction is controlled by a number of factors, including bed material, the presence of debris in basal ice and hydrological systems that modulate effective normal stresses. However, basal friction and slip remain poorly understood or constrained by observations 1,3,4 . Such observational constraint of friction and slip is critical for the verification of ice-bed boundary condition assumptions in ice-dynamics models, which are required to reduce uncertainty in corresponding sea-level rise projections [5][6][7] .
Previous contributions to address this critical observational void come from laboratory-based experiments [8][9][10][11][12][13] , geophysical studies [14][15][16][17][18][19][20][21][22][23][24] and borehole measurements 25,26 . However, to date, there have been challenges with such approaches. Laboratory experiments provide insight into fundamental physical properties of the bed material (rock, till, debris-bearing or 'clean' ice) 10 and ice-bed interface interactions 12 but are limited by scale and the diversity of natural glacier beds. Geophysical studies have measured the in situ bed strength but with sparse spatial and temporal resolution 14 . Borehole measurements of slip are not only sparse but have not been accompanied by measurements of shear stress, making quantitative interpretations difficult. The ice streams and outlet glaciers that contribute the majority of ice flux into the oceans probably have active, spatially and temporally varying hydrological systems 27,28 , perturbing basal friction and slip over short time and length scales. An observational void therefore remains.
Here we address this observational void by using icequakes to provide spatially mapped, in situ observations of both frictional drag and slip rate at the bed of an ice stream. These icequakes are generated by the sudden release of strain at or near the ice-bed interface. The dataset analysed comprises 100,000 icequakes 29 from Rutford Ice Stream (RIS), Antarctica (Fig. 1). The icequakes nucleate approximately at the centre of the ice stream, where the dominant source of drag is postulated to originate from the bed rather than from the Article https://doi.org/10.1038/s41561-023-01204-4 within expected physical limits. Effective normal stresses remain below the maximum ice overburden pressure, which is the upper possible limit for the average effective normal stress over the entire fault. The observed shear stress ranges from ~10 4 to 10 7 Pa. If the icequake cluster locations, or sticky spots, contribute more drag than the surrounding bed, then sticky-spot shear stresses could theoretically have a much higher limit than the average bed shear stress. Although bed shear moduli vary substantially between clusters, the majority of the clusters' shear moduli agree with two of the few previous seismically derived in situ measurements, which range from 20 to 70 MPa, from Whillans Ice Stream 14,19 . Additionally, measurements do not exceed the shear modulus of ice.
Slip rates show smaller variations in amplitude, from ~0.2 to 1.5 m d −1 , but have higher associated uncertainties due to their dependence on both shear modulus and fault area. While a number of clusters exhibit time-averaged slip rates approximately equal to the steady-state surface velocity of RIS (dashed line, Fig. 2f) 33 , other clusters have significantly lower seismic slip rates. Figure 3a-c shows the variation in effective normal stress, shear stress and slip rate for the entire experiment duration. Histograms of the stress and slip-rate distributions are shown in Fig. 3d-f. The normal and shear-stress histograms show a bimodal distribution, with more than two-thirds of the icequakes having effective normal stresses lower than ~5 × 10 5 Pa and shear stresses lower than 2 × 10 5 Pa. Conversely, the slip rates exhibit a unimodal distribution, tailing off below 0.2 m d −1 and above 1.5 m d −1 .
Individual icequake clusters switch on and off, being active for the order of hours to days (Fig. 2). Within single clusters, bed friction and slip are modulated by signals with dominant periods of ~6 to 12 h (Fig. 2). shear margins, justified by the width of RIS being an order of magnitude greater than its thickness. These icequakes nucleate in clusters that are highly repetitive (Extended Data Fig. 1), with near-constant inter-event times of the order of 100s of seconds and icequake clusters active for hours to days 29 . These icequakes are inferred to be at the bed from: their hypocentral depths, the consistent flow-and bed-parallel orientation of their double-couple focal mechanism slip vectors and full-waveform modelling of a typical RIS icequake source 20,[29][30][31] . The tight spatial clustering and repetitive nature motivate our use of a rate-and-state friction law in combination with icequake observations to investigate the glacier-sliding process. This rate-and-state friction law 32 also enables the calculation of other basal parameters including bed shear moduli and insight into the modulation of glaciological effective normal pressures.

Observed ice-bed friction and slip rate
The icequake source properties and inter-event times are used in combination with a rate-and-state friction law to calculate: fault effective normal stress (σ ); total frictional shear stress, or drag per unit area (τ); shear modulus (G bed ); slip (d) and slip rate (v slip ) for seismically active regions of the bed of RIS. Note that fault effective normal stress is not equivalent to the effective normal pressure conventionally used in glaciology (Discussion and Methods provide more details). Moreover, our results are unaffiliated with a particular bed material. For example, measured bed shear moduli could apply to both bedrock or till. Figure 2 shows these results for a representative subset of icequake clusters. Fault effective normal stress, shear stress and shear modulus (Fig. 2a-c)  However, although this alludes to tidal modulation of basal friction, and indeed surface velocities are known to be modulated by tidal frequencies 33,34 , we cannot decipher a clear relationship between tidal signals propagated 40 km upstream from RIS's grounding line and our signals 29 . We therefore do not discuss any link with tidal signals further.
The spatial distribution of the temporally averaged basal shear stress, slip rate and fault radius for each cluster over a 7 × 6 km region are shown in Fig. 4. Shear stresses are largest at the clusters farthest upstream, approximately where the bed properties are inferred to transition from unconsolidated to consolidated till 35 (or weakly consolidated sedimentary bedrock 36 ) (pink dashed line, Fig. 1) and where the bed has shorter wavelength topography than upstream that probably inhibits ice flow. Average slip rate is spatially consistent across all clusters. This is expected, as our study site is located near the centre of the ice stream, with no spatial variation in surface slip rate 33 . Fault radius, defining the area of an icequake cluster sticky spot, is also measured (Fig. 4c). Fault radii indicate that individual seismically active sticky spots have areas <2,800 m 2 . Only a small number of sticky spots are active at any instant, with the seismically active proportion of the bed at any given time shown in Fig. 3g. The average seismically active percentage of bed directly beneath the network is ~0.07%. This suggests that regions of sufficiently high basal friction to generate seismicity are confined to the minority of the bed at a given point in time, yet support significant basal drag. Regions absent of seismicity between icequake clusters probably also contribute to the basal drag, presumably providing the dominant source of aseismic drag upstream of the unconsolidated-consolidated sediment boundary (pink dashed line, Fig. 1). The total drag accommodated by the icequakes averaged over the whole bed study area would equate to ~500 Pa, whereas basal drag inferred from models combined with surface observations is typically ~50 kPa for RIS, implying that much of the overall drag is accommodated aseismically.

Implications of friction and slip-rate observations
The most important, immediate finding of this work is the ability to observe in situ frictional shear stress and slip rate, the two critical parameters for constraining the basal drag boundary conditions of ice-dynamics models. Our approach could be applied to any glacier that generates icequakes. Most fast-moving glaciers probably generate such icequakes, with the majority of glaciers on which seismometers have been deployed exhibiting at least some basal seismicity 16,29,[37][38][39][40][41][42][43][44][45] . Seismic tremor associated with sliding can also occur 19,46 , thought to initiate at the boundary between the conditionally stable and unstable regimes of the rate-and-state friction model 19,32 . Indeed, the premise of this study was inspired by such observations 19 . However, due to the inability to extract both corner frequency and inter-event time information from tremor, it cannot be used to measure shear stress and slip using our approach.
Our confidence in the frictional shear-stress and slip-rate measurements is founded partially on the uncertainty amplitudes but also fundamentally on the agreement between the observed basal slip rates and Global Navigation Satellite System (GNSS) derived surface displacement 33 . This agreement validates assumptions of slip-dominant rather than deformation-dominant flow at RIS and the use of a rate-and-state model and assumptions of the icequake source properties. The small discrepancy between the surface and basal slip rates is primarily due to uncertainty, except for a minority of particularly sticky spots. These sticky spots exhibit particularly strong frictional drag that significantly inhibits local ice flow, albeit for short durations of the order of hours to days.
Observed basal shear stresses are of the order of 10 4 to 10 7 Pa, acting at sticky spots with diameters of the order of 10 to 60 m (Fig. 4). Basal shear stresses of the order 10 5 Pa are typical values used in ice-dynamics models 47 and laboratory experiments 8 for RIS's surface slip rate of ~1.1 m d −1 (ref. 33). Basal shear stresses of 10 6 to 10 7 Pa might initially appear inconsistently high compared to models and experiments 48 . However, these high friction sticky spots are spatially small compared to the total bed area. Our results therefore imply that certain icequake clusters accommodate a considerable proportion of the total basal drag.

Mechanisms that generate icequakes
We propose that the icequakes are generated by at least one of two mechanisms or sliding regimes. The presence of two sliding regimes is motivated by the physical system and the bimodal distributions observed in Fig. 3d,e. The regimes (Fig. 5) are: regime I, rock-on-rock friction between ice-entrained clasts and bedrock at the fault interface and regime II, where clasts plough through till, with failure accommodated by a till-on-till fault interface. Clasts are pieces of rock partially entrained into the ice (Fig. 5). The presence of such clasts is discussed in the literature 8,12,13,49,50 . The motivations for these clast-based icequake models are that they can explain the rate-weakening friction required to generate icequakes 13 , that clasts are required to generate icequakes in laboratory environments 12 and that such icequakes probably originate within one seismic wavelength of the ice-bed interface 20 . We suggest that the highest effective normal stress icequakes exhibit regime I sliding because this regime allows for the average effective normal stress over the entire fault area to be concentrated over much smaller clast-bedrock contact areas. Similarly, we postulate that the lower effective normal stress icequakes are associated with regime II sliding, although we cannot rule out that all icequakes are generated via regime I.

Effective normal stress versus effective fluid pressure
Our results imply significant temporal variation in basal effective normal stress. Full numerical rate-and-state simulations also corroborate the range of observed normal stresses (Supplementary Information Section 2). Such increases and decreases in effective normal stress are inferred to be caused by corresponding decreases and increases in basal water pressure 16,51,52 . However, while the icequake-derived effective normal stresses, σ , averaged over the entire fault are equivalent to the average glaciological effective pressure, P eff = P ice overburden − P water , within that same fault area, asperities and bed heterogeneity on length scales shorter than the fault diameter could significantly perturb local effective normal pressures. Although all our measured effective normal stresses remain below the ice overburden pressure, current hydrological models supported by other observations at ice streams 11,53 cannot reconcile glaciological effective pressures greater than ~0.5 MPa for expected till porosities, except potentially if they are transient 36 .
Sparse observations of effective normal pressures at RIS from borehole measurements find P eff ≈ 0.2 MPa (ref. 50), although till acoustic impedance measurements at RIS suggest that dewatering is possible 23 . Dewatered till would imply P eff = P ice overburden . Similarly, modelling 54 and seismic observations 55 elsewhere suggest that dilatational strengthening of till could also potentially provide a mechanism for generating such high effective normal pressures. Alternatively, if the consolidated material is actually approximately impermeable sedimentary rock, then melt cavities formed could be isolated, generating high effective pressures locally in between fracture-plucking events 56 . It is therefore physically possible to generate sufficiently high effective pressures transiently. Therefore, our highest observed effective normal stresses suggest either that our understanding of bed characteristics and associated physical models may have to be revisited, at the very least for RIS, if we are to adequately describe the transient phenomena that we observe at metres to tens of metres scales; or that the rate-and-state model does not adequately describe icequake physics. We suggest three possible explanations for resolving the discrepancy between observed and theoretical maximum effective normal stresses. First, prominent bedrock outcrops could significantly inhibit ice flow, allowing stoss-side effective normal stresses of the order of at least 1 MPa to develop 36 . A second explanation again lies with bedrock outcrops, whereby impermeable bedrock might inhibit the transport of fluids, facilitating dewatered regions. A third explanation is that till porosities are far lower than conceived in current models of bed properties, possibly supported by till impedance measurements 23 . These suggestions are not exhaustive. We suggest here only that our results motivate new models to explain such observations.
If effective normal stresses are related to P eff at sticky spots, then they provide an observational foundation for calibrating basal fluid pressures assumed in laboratory experiments 8 , ice-dynamics models 47 and glacier basal hydrology and tidal forcing 33,34 .

Enhanced knowledge of RIS bed conditions
The considerable variation in bed properties observed at RIS are presented as an example of the enhanced knowledge of the bed properties that our approach provides. First, for the unconsolidated till (label 1, Fig. 5) and much of the more-consolidated material (label 2, Fig. 5), the effective normal stresses are too low to generate the unstable stick-slip conditions required for icequake nucleation. Within consolidated regions (label 2, Fig. 5), there are small zones that become seismically active if the effective normal stress is sufficiently high (label 3, Fig. 5). These sticky spots turn on and off, modulated by changes in effective normal stress, bed strength and bed material.
Frictional shear stress at an individual sticky spot can vary temporally by up to an order of magnitude and spatially by several orders of magnitude. This variation occurs over the order of hours and 100s of metres. This implies that there are both till and bedrock outcrops, combined with an active hydrological system or permeable bed capable of such variable spatial and temporal variations in basal fluid pressure. The regions exhibiting the highest shear stresses are at the upstream edges of local topographic highs near the unconsolidated-consolidated bed boundary (label 5, Fig. 5). This is probably because resistive stresses of these materially stronger, topographic highs can accommodate more basal drag. Some regions of consolidated bed might contain pockets of melt water (label 4, Fig. 5). However, we cannot observe such a phenomenon using seismicity as these patches would have an approximately zero shear modulus.

Broader implications for ice-sheet sliding
Our results show that much of the basal drag at an ice stream can be accommodated within small zones of significantly higher-than-average friction. This could have a profound impact on how sliding is formulated in ice-dynamics models. However, although friction varies significantly, average basal slip rates remain predominantly stable at RIS. This is encouraging for current modelling efforts because if the temporally and spatially averaged slip rates are approximately constant, then perhaps such models are not required to be sensitive to small-scale, rapid variations in bed friction. Our observations quantify the highly variable bed properties over a sufficient duration required to test such a hypothesis.
Another important question to address is how our approach could be implemented at ice-sheet scale. One could deploy temporary seismic arrays on important ice streams and outlet glaciers for short durations. A number of targeted deployments would allow verification of the link between surface and basal velocity at ice-sheet scales 2 .
A further question that this study raises is could a rate-and-state friction model used for the icequake-sliding analysis also be used as a mathematical basis for informing sliding laws used in ice-dynamics models more generally. Such models have recently been proposed to describe surging glacier behaviour 54,57 , which have been validated at laboratory scale 12,13 and evidenced from seismic observations 55 . The rate-and-state model meets the conditional stability requirement, not allowing runaway acceleration of a glacier (Supplementary Text). Our results and the other recent works cited here suggest that rate-and-state friction models could be used to inform future glacier-sliding laws more generally.
Finally, these icequake observations can aid the understanding of earthquake mechanics more generally. Even the smallest magnitude icequakes (M w ≈ −1.5) have high signal-to-noise-ratios and so could elucidate any lower limits on the fundamental size of earthquake nucleation for given fault properties [58][59][60] . Additionally, icequakes in this study have stress drops that vary with magnitude, contrary to magnitude-invariant stress drops observed for larger earthquakes 61 .
Our findings show that icequakes can provide the critical observations required to constrain the highly variable friction at the bed of an Antarctic ice stream. Applying such observational constraint to ice-dynamics models would reduce uncertainty in corresponding sea-level rise projections.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions Residual topography (m) 10 20 30 Fault radius (m) 78

Fig. 5 | Schematic diagram summarizing the findings of this study in relation to basal friction and slip with bed characteristics.
Bed properties are labelled in the legend. Numbered labels are referred to in the text. Note that features are not to scale but arranged approximately according to spatial trends in Fig. 4. Regime I and regime II are shown schematically, with regime I being clast-on-rock icequake slip behaviour and regime II being till-on-till slip behaviour. Note also that consolidated till could instead be lithified sedimentary rock.

The icequake dataset
This study uses 100,000 icequakes at the bed of Rutford Ice Stream (RIS), Antarctica. An example of such an icequake arrival can be found in Extended Data Fig. 1. These data were collected over the period of November 2018 to February 2019. The icequakes are detected using QuakeMigrate 38,84 and relocated using NonLinLoc 64 . A full description of the detection, location and clustering analysis of this seismicity can be found in ref. 29. The hypocentral depths, orientation of focal mechanisms and full-waveform modelling provide us with confidence that these icequakes are associated with sliding within one seismic wavelength (~10 m) of the bed 20,29,30 .

Observable parameters from stick-slip icequakes
Earthquake source models can be used to calculate the size of the earthquake, its duration, the fault radius and the shear-stress drop associated with the release of seismic energy. These observable parameters are required for any analysis of frictional behaviour at the bed of glaciers using icequakes. The methods we use to obtain these parameters from the icequake signals are described below.
Seismic moment. The seismic moment, M 0 , of an earthquake is proportional to the energy released by that event. Seismic moment release can be defined as 65 where ρ is the density of the medium at the earthquake source, v i is the velocity of the seismic phase i (P or S), r is the hypocentre receiver distance, Ω 0 is the long-period spectral amplitude, A rad,i is the amplitude of radiation of seismic phase i for the particular source receiver azimuth and take-off angle and C free-surface is the free-surface correction term, which depends upon the angle of inclination of the seismic phase arrival at the surface. For this study, we assume typical ice values of ρ = 917 kg m −3 , v p,ice = 3,841 m s −1 , v s,ice = 1,970 m s −1 . A rad,i is calculated as described in ref. 20, based on the assumption that all the icequakes in this study are double-couple sources with strikes aligned with the ice flow direction. This assumption is based upon previous observations at Rutford Ice Stream 20,29,30 . Ω 0 is calculated by fitting a Brune source model to the noise-removed spectrum of the icequake 66 .
Corner frequency. The spectrum of an earthquake contains more information than just the long-period spectral amplitude. If one assumes that an earthquake's spectrum can be described by a Brune model 66 , then one can also measure the corner frequency, f c , of the earthquake. However, an earthquake's spectrum is also particularly sensitive to seismic attenuation. Seismic attenuation, often described by the quality factor, Q, reduces the amplitude of an earthquake spectrum nonlinearly across all frequencies. If path attenuation is poorly constrained, then it can lead to detrimental uncertainty in the measured corner frequency, as evidenced by the trade-off between Q and f c in the Brune model 66 where Ω(f) is the amplitude of the spectrum for a certain frequency f, and t is the travel time.
To obtain an accurate measurement of corner frequency, we therefore use a linearized spectral ratios method to constrain Q. This spectral ratios method isolates the path effects from the source effects. An example of the linearized Brune model fit and the observed spectrum for an example icequake is shown in Extended Data Fig. 1c. We obtain estimates of Q from this method of the order of 200 to 800 (Extended Data Fig. 2b), which are in agreement with other measurements for Antarctic ice 67 . This then allows equation (1) to be fit to the earthquake spectrum with only Ω 0 and f c as variables. We find that the icequake corner frequencies at RIS fall approximately within the range of 40-100 Hz (Extended Data Fig. 2c).

Fault radius and stress drop.
One can estimate the fault radius, R, and stress drop, Δτ, of an earthquake from the corner frequency.
The relationship between corner frequency and fault radius, R, is given by 68f where f c is the spherically averaged corner frequency for the earthquake, β is the shear-wave speed near the source and k i is a constant relating the spherically averaged corner frequencies for a specific fault model for the seismic phase i. Here we use the fault model of ref. 69, which gives k P = 0.38 and k S = 0.26 for a rupture speed of 0.98 β. We let β equal the shear velocity of ice (1,970 m s −1 (ref. 30)). For clast-on-bedrock slip (regime I, Fig. 5), this is valid as rupture will propagate through the bedrock and the ice that the clasts are embedded within, with us only observing the rupture propagation through the ice. For till-on-till slip (regime II, Fig. 5), our assumption of β is probably an overestimate, resulting in an overestimate of fault radius. As the seismic properties of the till are unknown, we are limited in assigning a lower value of β for any regime II events. We assume a symmetric circular fault for this analysis. We therefore calculate average corner frequencies for each event based on the corner frequencies observed at all receivers. The potential effects of the symmetric circular fault assumption are shown in ref. 70.
The uniform stress drop of an earthquake can then be found using the fault radius and the relationship given by ref. 71 We now have all the observable parameters required to constrain a friction model at an icequake source.

Using a rate-and-state friction law for deriving frictional shear stress and slip from icequakes
Calculating shear stress. Earthquakes are typically generated as the result of stick-slip frictional instabilities at a fault interface 32 . We hypothesize that icequakes associated with sliding at the bed of a glacier can be described by a similar model. For our investigation, we assume that the fault interface is at or near (<1 wavelength) the icebed interface. Schematic diagrams describing the model are given in Extended Data Fig. 3. Within this framework, we can apply the following rate-and-state friction law given by ref. 32 where τ is the total frictional shear stress, μ 0 is the steady-state friction coefficient at v = v 0 , v is the slip velocity, v 0 is a reference velocity, defined in this case to be the background slip rate, a and b are material properties, ℒ is the characteristic slip distance over which the system returns to steady state and renew surface contacts and θ is the state variable. The variation of the state variable, θ, through time can be defined by the ageing or the slip laws 72 , given by Article https://doi.org/10.1038/s41561-023-01204-4 The state variable, θ, represents the characteristic contact lifetime of a fault. To apply the rate-and-state model to the stick-slip icequake system in a mathematically tractable way, we assume that the state variable of the system is constant over the duration of an icequake cycle, that is, ∂θ ∂t = 0 through all time during a cycle. For a destructive frictional failure process, θ probably changes with time during earthquake nucleation and as the fault heals. However, for the icequake generation mechanisms proposed in this study (Fig. 5), damage at the fault interface that affects the frictional properties is probably less significant than at traditional earthquake fault interfaces. This lack of damage is evidenced to some extent by the highly repetitive nature of the icequakes 20,29 . We assume that at least part of the icequake patch is near steady state, or approximately at steady state if it slips sufficiently fast. A caveat to this is that some of the icequake patch could have remained below the steady-state sliding limit, which we do not explore this here. Overall, we deem the approximation of ∂θ ∂t = 0 between individual icequakes as acceptable in this case. This assumption can be used to find the state variable, θ, as an expression of v, ℒ from either the ageing law (equation (6)) or the slip law (equation (7)), which both yield Equation (5) can then be reduced to a rate-dependent friction law, given by The coefficient of friction in equation (9) can then be thought of as μ = μ 0 + Δμ, where μ 0 is the static friction component, and the dynamic friction component, Δμ, is given by which when multiplied by the effective normal stress, σ, can be assumed as equal to the icequake stress drop (equation (13)). One can then parametrize equation (9) in such a way so that it can be solved for individual icequakes. We take μ 0 = 0.4, a = 5 × 10 −3 and b = 15 × 10 −3 from ref. 19. We approximate the ratio of the instantaneous sliding velocity to the reference velocity, where d is the slip associated with an event (unknown), T is the slip duration, which we approximate to be equal to the inverse of the icequake corner frequency, f c 73 , and t inter-event is the time between two consecutive icequakes. The correspondence of these parameters to the stick-slip cycle is shown in Extended Data Fig. 3b. With this parameterization, the velocity ratio then becomes Assuming that the friction at the interface is velocity weakening and therefore unstable, one can then assume that the dynamic part of equation (9) is equal to the stress drop measured during an icequake, Δτ (ref. 74). One should note that this assumption implies that all the dynamic stress release during slip is accommodated seismically (red shaded region of Extended Data Fig. 3b). However, there is also frictional shear stress present that cannot be measured directly using stress-drop measurements. We also assume a seismic radiation efficiency of 1, which is obviously an approximation, with the actual seismic radiation efficiency unknown. Although the radiation efficiency will, in reality, be <1, due to thermal heating and the generation of additional surface area during abrasion, fracture tip energy and other phenomena such as off-fault cracking are probably insignificant in comparison to standard earthquakes 75 , so we deem our first-order approximation as reasonable in this case. For tectonic earthquakes, the seismic radiation efficiency typically might be of the order of 0.1 (for example, see ref. 76). If the icequake seismic radiation efficiencies were similarly low, then this would be approximately equivalent to reducing the magnitude of M 0 by a factor of 10. Sensitivity analysis in the Supplementary Text suggests that such a reduction in M 0 would reduce the shear stress, τ, by an order of magnitude, but the slip velocity, v slip , would be reduced only by a factor of 3. Assuming velocity-weakening friction and a radiation efficiency of 1 results in the definition of the effective normal stress at the fault interface, given by Once we know the effective normal stress, σ, we can find the overall shear stress on the fault, τ, from equation (9).
We emphasize that the effective normal stress, σ , is the normal stress on the fault, which is not necessarily equivalent to a traditionally defined glaciological effective pressure, P eff = P ice − P water . The fault effective normal stress, σ, is the effective normal stress that acts over the fault area, A fault , derived from the earthquake corner frequency (equation (3)). The actual normal stress acting through clasts in contact with the underlying contact surface might increase the normal stress acting through these clasts (sliding regime I, Fig. 5). However, fault-average normal stress, σ , must be equal to the average glaciological effective pressure, P eff , over the same area of the bed.
Calculating slip. The second glaciologically important parameter to measure at the bed is the slip and hence the basal slip rate. To calculate slip, we assume that while an individual icequake cluster is active, all (or at least the vast majority of) slip is accommodated seismically. This is probably the case for RIS, as evidenced by the close agreement between surface slip rate and seismically measured basal slip rates (Fig. 2f). Calculating the basal slip, d, from an icequake is challenging because one first has to determine a method of estimating the bed shear modulus, G bed , because the slip is given by where M 0 is the seismic moment released by an earthquake and A is the area of the fault. The bed shear modulus, G bed , is calculated by assuming a further behaviour of the rate-and-state friction law. This behaviour is that an earthquake can nucleate only if it is in the unstable regime. In this study, we assume that the temporally averaged driving shear stress at the fault varies over longer timescales than the icequake inter-event time, with the shear stress at which the fault fails governed by the effective normal stress acting on the fault, σ. The approximately constant inter-event time between individual consecutive icequake pairs (Fig. 2e) within a single cluster validates this assumption. The effective normal stress at which a fault becomes unstable is defined as the critical normal stress, σ c , with velocity-weakening behaviour prevailing above this stress. Although earthquakes can occur only in the unstable regime, the remarkable consistency in event time and stress drop between tens to hundreds of consecutive icequakes allows us to make the assumption that for a given normal stress, the icequakes probably fail at approximately the critical normal stress, σ c . This hypothesis is probably the https://doi.org/10.1038/s41561-023-01204-4 result of an approximately homogeneous distribution of till grain sizes that makes an ice-bed fault far more homogeneous than a traditional earthquake fault, which exhibits significantly higher randomness. If we make this assumption, then the equation governing σ c can be written as an equality rather than a lower bound. σ c is then given by ref. 32 where η v driving is the product of the radiation damping parameter and the load velocity and η v d /ℒ represents the inertial stability limit 19 . For our system we find η v d /ℒ «k/ℒ for all σ c . This term can therefore be dropped in equation (15). k is the spring constant of the system (Extended Data Fig. 3a), which is given by where G * is the effective shear modulus of the bimaterial interface 19 and R here is the radius of the fault, which can be found from the icequake corner frequency, if assuming a symmetric, circular fault 69,70 . However, this equation still has two unknowns: G * , the effective shear modulus that we require to calculate the slip, and ℒ, the critical slip distance, otherwise referred to as the state evolution distance. For the purposes of this study, we approximate ℒ to remain constant but allow G * to vary with effective normal stress, which from granular material theory [77][78][79][80][81] is assumed to take the generic empirical form where A, n and C are constants to invert for. We use a least-squares approach to minimize the function where σ c and R vary for each icequake and A, n, C and ℒ are varied to minimize the function. σ c is taken to be the effective normal stress for the first 100 icequakes when a cluster becomes active, as calculated using equation (13). These parameters are found to be A = 22,000, n = 0. 78, C = 8,200 Pa and ℒ = 7.7 × 10 −5 m, with the result of the minimization shown in Extended Data Fig. 4. Now ℒ can be substituted into equation (15) to find k, which can then be used in equation (16) to find the bimaterial shear modulus, G * . It should be noted that (a − b) and ℒ trade off in both equation (18) and equation (15), therefore if the assumed (a − b) is different to the real value, then it will affect the magnitude of ℒ but not the value of k and hence G * . The shear modulus of the bed, G bed can then be found using the Poisson ratios of ice (1/3) and till (0.49), which gives G * ≈ 3.5 G bed (ref. 19). Granular material theory, or at least the relationship of equation (17), is thought to still hold for clast-over-bedrock sliding because the shear modulus will still be related to some exponent, n, of σ , even if that exponent were ~0. Equation (14) can then be used to find the slip, d, associated with a single icequake, for the effective normal stress applied to the fault at that particular time. We also calculate the approximate slip rate associated with these highly repetitive icequakes. If one assumes that all the slip when an icequake cluster is active is accommodated seismically, then one can calculate the slip rate per day, v slip v slip = d t inter-event (19) The methods described above allow us to calculate the total shear stress, τ, and the slip, d, at the bed. These two parameters can provide observational constraint on ice-dynamics models of ice streams.
A note on assumptions. A number of assumptions are made to make the derivation of basal shear stress and slip from icequake observations and a rate-and-state friction model mathematically tractable. Extended Data Table 1 summarizes all the assumptions made. There are several assumptions that warrant particular emphasis. The first is the assumption that all slip at an individual sticky spot is accommodated seismically while that cluster is active. The highly repetitive nature of the icequakes (Extended Data Fig. 1 and ref. 29), with approximately constant inter-event times between consecutive icequakes in a cluster, is indicative of the stability of each sticky spot (Fig. 2), justifying this assumption. Second, a Brune model 66 is assumed to describe the earthquake source characteristics. Whereas such a model is probably an approximation for the complex physics of earthquake rupture, it is a common assumption for other earthquake studies that is probably also a valid approximation for the stick-slip icequakes presented here. Third, we approximate that the time derivative of the state variable in the rate-and-state friction model, ∂θ ∂t , equals zero during an individual icequake cycle. This approximation is valid if slip on the fault is sufficiently fast and if little damage occurs at the fault, compared to more complex earthquake faults. Obviously, this is only an approximation, as damage does probably occur at the fault, at least for the clasts-over-bedrock slip case (regime I, Fig. 5). We also assume that (a − b) is constant in space and time. This is probably not the case, but without this assumption, there are too many unknowns to solve for. We show in Supplementary Information Section 3 that assuming a fixed (a − b) and the approximated amplitude of (a − b) do not significantly affect our findings, especially the calculation of slip rate. Furthermore, an underestimation bias in slip may be introduced by the assumption of no fault frictional heating. Fault frictional heating would reduce the seismic radiation efficiency from our approximation of one 75 . The final assumption we emphasize here is that we assume that the icequakes at the beginning of an icequake cluster nucleate at approximately the critical normal stress for nucleation, σ c , rather than at some arbitrary value above it. The icequake slip calculations are dependent upon this assumption. This assumption would not be valid for sporadic earthquakes on complex faults, as shear stresses could build to different values before failure for each earthquake, even with constant effective normal stresses, due to fault heterogeneity. Nor would it necessarily be valid if the driving shear stress were perturbed over time scales shorter than the inter-event time, for example, by interactions with other icequake clusters. However, although icequake faults still exhibit a degree of heterogeneity due to an inhomogeneous distribution of clasts, this heterogeneity has negligible impact upon the consistency of both the inter-event times and shear stresses between consecutive icequakes at a given sticky spot (Fig. 2). Furthermore, there are only a small number of active icequake clusters at any given time, which are spatially isolated from one another. The consistency in inter-event times and shear stresses observed in our data, in agreement with similar, laboratory-generated icequakes 12 , provides us with confidence in our assumption of icequakes nucleating at the critical nucleation stress, σ c .

Data availability
All the seismic data used in this analysis will be deposited in the IRIS seismological data repository. The icequake catalogue used for this analysis is available from the UK Polar Data Centre 82 , with details on how this catalogue was constructed given in the peer-reviewed publication 29 . The results of this work are shared via an online repository (https://doi.org/10.5281/zenodo.7870307).