Experimental constraints on subglacial rock friction

Abstract

Subglacial rock friction is an important control on the sliding dynamics and erosive potential of hard-bedded glaciers, yet it remains largely unconstrained. To explore the relative influence of basal melt rate, effective stress and ice temperature on frictional resistance, we conducted abrasion experiments in which limestone beds were slid beneath a fixed slab of ice laden with granitic rock fragments. Shear stress scales linearly with melt rate and cryostatic stress, confirming that both viscous drag and effective stress are first-order controls on the contact force in drained conditions. Furthermore, temperature gradients in the ice increase the contribution of viscous drag on basal shear stress. In all experiments, the relationship between melt rate and shear stress is best explained by a model that accounts for the effects of regelation and viscous creep on the bed-normal drag force. We interpret this to mean fluid flow around entrained clasts contributed to basal drag even at subfreezing temperatures. Incorporating premelting dynamics into the Watts/Hallet model for subglacial rock friction, we find that the predicted debris-bed drag decreases by approximately an order of magnitude, with a corresponding ~3.5 × increase in the transition radius. This is lower than we observe for ice slightly below the pressure melting point.

1. Introduction

Friction between entrained basal debris and the bed affects the erosion rates (Schweizer and Iken, 1992; Hallet and others, 1996; Koppes and Hallet, 2006; Herman and others, 2015; Koppes and others, 2015), sliding speed (Hallet, 1981; Shoemaker, 1988) and slip stability of glaciers (Zoet and others, 2013; McCarthy and others, 2017; Lipovsky and others, 2019). Basal ice is typically assumed to be debris-free in glacier flow models, meaning basal resistance arises solely from ice deformation around bedrock obstacles (Lliboutry, 1968; Nye, 1969; Iken, 1981; Zoet and Iverson, 2015, 2016). However, observation shows that glaciers transport substantial bedload (Kirkbride, 2002), and theory and field measurements indicate subglacial rock friction may be significant, albeit poorly constrained (Iverson and others, 2003; Cohen and others, 2005).

A central assumption in theories describing subglacial rock friction is that inclusions in ice are surrounded entirely by a thin water film (Boulton, 1974; Hallet, 1979b, 1981; Hindmarsh, 1996; Cohen and others, 2005; Emerson and Rempel, 2007). For temperate ice, classic regelation and creep theory predicts the presence of a film due to ice pressure melting along grain boundaries (Nye, 1969, 1973; Kamb, 1970). However, considerations of premelting indicate a liquid film exists even at subfreezing temperatures (Gilpin, 1979; Dash and others, 1995; Rempel and Worster, 1999; Rempel and others, 2001; Rempel and others, 2004; Dash and others, 2006). Since ice is nowhere in contact with the clast, stresses in the ice must be transmitted to the rock through the fluid surrounding the clast. Fundamentally, discrepancies between interfacial models stem from their assumptions regarding the distribution of water pressure along the clast boundary and its controlling mechanisms.

To date, the competing models of Boulton (1974) and Hallet (1979b, 1981) have driven theoretical framing of subglacial rock friction (Fig. 1). Boulton (1974) first proposed that the bed-normal contact force beneath an abrading clast, F _c, is the product of the effective pressure in the film, N (i.e. cryostatic pressure, P _i, minus the water pressure, P _l) and the surface area of the particle, A, such that

(1)$$F_{\rm c} = F_{{\rm eff}} = AN{\rm,} $$

where F _eff is the contact force due to effective stress acting on the clast. For a spherical particle in contact with the bed, A = π(2tr−t ²), where r is the radius and t is the thickness of a water film that separates the ice from the bedrock. However, as Hallet (1979b, 1981) recognized, this model neglects the contribution of the buoyant force due to gravity, F _b, and the drag force, F _d, induced as ice creeps and/or regelates around suspended particles. Accounting for these additional components, the bed-normal contact force is more completely expressed:

(2)$$F_{\rm c} = F_{\rm b} + F_{\rm d} + F_{{\rm eff}}.$$

Fig. 1. Conceptual models developed by (a) Boulton (1974) and (b) Hallet (1979b, 1981) to describe the bed-normal contact-force beneath an abrading clast. A water film of thickness t separates ice from bed, and a second film envelopes the clast with thickness δ in Hallet's (1979b, 1981) model. In Boulton's (1974) model (a), the film around the clast is assumed to be sufficiently thin so that liquid pressure in it is effectively zero. F _eff is the force due to effective stress, F _d is the viscous drag force, F _b is the buoyant force due to gravity, P _l is the water pressure in the films, P _i is the cryostatic stress, r is the radius of the clast, u _v is the vertical ice velocity and A is the cross-sectional area of the clast at height t. The schematic is adapted from Byers and others (2012).

In this later formulation, F _eff = 0 is implicitly assumed, implying that deviatoric stresses arise only from the combined influence of F _b and F _d and that cryostatic stress (and by extension glacier thickness) has no bearing on the erosive potential or frictional properties of basal debris. For small clasts (r < 0.5 m), F _b is negligible (Hallet, 1979b). To quantify the drag force, Hallet adapted Watts's (1974) analytic solution for a spherical inclusion moving through temperate ‘Newtonian’ ice by regelation and creep to a half space, such that

(3)$$F_{\rm d} = \displaystyle{{\Phi 4\pi \eta _{\rm i}{\rm \;} u_{\rm v}r^3} \over {r_{^\ast}^2 + r^2}},$$

where Φ is a bed influence factor, η _i is the viscosity of ice, u _v is the bed-normal ice velocity, r is the radius of the particle, and r _* is a transition radius analogous to the controlling wavelength in glacier sliding. Watts (1974) defines the transition radius as

(4)$$r_{^\ast} = \left( {\displaystyle{{3\eta_{\rm i}{\rm \;} CK} \over {\rho_{\rm i}{\rm {\cal L}}}}} \right)^{1/2},$$

where C is the Clapeyron slope, K is the mean thermal conductivity of ice and rock, and $\rho _{\rm i}{\rm {\cal L}}$ is the volumetric latent heat of fusion for ice.

The prefactor Φ in Eqn (3) corrects for the influence of a bed on the viscous flow field around a clast. Hallet (1981) based his initial estimate of Φ = 2.4 on the analysis of Goren (1970). However, the boundary conditions prescribed therein are not characteristic of the subglacial environment (Cohen and others, 2005), and empirical evidence regarding the magnitude of this parameter is limited and often contradictory. For example, Iverson (1990) observed Φ = 0–5 for non-idealized granite clasts during laboratory abrasion experiments. Byers and others (2012) found that for the case of a stationary spherical inclusion in contact with the bed, F _d was approximately twice the Stokes drag force induced by the flow of Newtonian ice around a slippery sphere (i.e. Φ = 1.8). Lee and Rutter (2004) found Φ ≈ 10 was necessary to reconcile their observations with theory when comparing experimental rock-on-rock wear data with the suspended sediment yield in a proglacial stream.

In subsequent years, Hallet's model (1979b, 1981) has been widely implemented in glacier erosion models (Herman and others, 2011; Anderson, 2014; Beaud and others, 2014; Egholm and others, 2017; Ugelvig and Egholm, 2018), and available experimental data support a linear relationship between F _c and u _v as predicted (Iverson, 1990; Byers and others, 2012). Although Hallet's model identified key controls on debris-bed friction, numerous lines of evidence indicate the problem remains unresolved. As Hindmarsh (1996) and Morris (1979) address in considerable detail, mathematical inconsistencies arise when Watts's (1974) analytic solution is considered in the context of common subglacial boundary conditions. Myriad complicating factors, including the presence of a far-field lubricating film, bed transmissivity, the formation of cavities along clast boundaries, impurities in the ice, incomplete refreezing during regelation, high debris concentrations, temperature fluctuations in basal ice and interfacial melting are pervasive in nature and likely preclude this model from being universally applicable. Moreover, not all available data agree with this model's predictions. Experiments conducted by Emerson and Rempel (2007) show a statistically significant relationship between cryostatic stress and basal drag for sand-sized particles in a melt-dominated regime, and pressure sensors emplaced beneath Engabreen Glacier recorded anomalously high shear stresses that also covaried with effective stress (Iverson and others, 2003; Cohen and others, 2005). Thus, it seems that a more nuanced treatment of subglacial rock friction is needed in order to capture the complexities of this system.

This lingering uncertainty provides an impetus to revisit this classic problem and better constrain the mechanisms controlling debris-bed drag. To this end, we designed laboratory experiments to simulate glacial abrasion under a range of realistic subglacial conditions, such as variable cryostatic stresses, basal melt rates and ice temperatures. We tested the relative influence of these controlling parameters on basal shear stress using non-idealized clasts. In the following sections, we incorporate the effects of premelting into the existing Watts/Hallet model for subglacial abrasion and then present our methodology and experimental results. Finally, we weigh our results against theory.

1.1. The effect of premelt on debris-bed contact forces

Recently, Rempel and Meyer (2019) postulated that premelting increases both the rate of regelation by an order of magnitude and the controlling obstacle size in glacier sliding by a factor of ~3.5. Given that basal drag in the Watts/Hallet model is a function of the regelation rate, this may have important implications for the drag force imparted against an entrained clast. Incorporating the same premelting condition into Watts' (1974) original analysis, we find that it predicts a significant decrease in F _d. The steps of this derivation are detailed in the Appendix, Eqns (A1–A10). We obtain a modified expression for drag force, F _d, imparted against the sphere,

(5)$$F_{\rm d} = \displaystyle{{4\pi \eta _{\rm i}u_{\rm v}r^3} \over {r_{^\ast}^2 + r^2}}$$

where

(6)$$r_{^\ast} = \left( {\displaystyle{{\rho_{\rm l}} \over {\rho_{\rm l}-\rho_{\rm i}}}\displaystyle{{3\eta_{\rm i}CK} \over {\rho_{\rm i}{\rm {\cal L}}}}} \right)^{1/2}.$$

Here ρ _l represents the density of liquid water, and ρ _i is the density of ice.

Notably, we find that the transition radius, r _*, is approximately $\sqrt {\rho _{\rm l}/\lpar {\rho_{\rm l}-\; \rho_{\rm i}} \rpar } \approx 3.5$ times larger than the transition radius predicted by Watts's (1974) original analysis (Eqn (4)). Figure 2 juxtaposes the relationship between clast radius and drag force predicted for both solutions, assuming η _i = 3.15 × 10¹² Pa s, u _v = 500 mm a⁻¹, ρ _l = 1000 kg m⁻³, ρ _i = 917 kg m⁻³, ${\rm {\cal L}}$ = 3.3 × 10⁵ J kg⁻¹, C = 0.074 K MPa⁻¹, K = 2.5 W m⁻¹ K⁻¹.

Fig. 2. A comparison of clast radius versus contact force for the Watts/Hallet model and a solution that accounts for the effects of premelting on regelation. Incorporating premelting dynamics into the Watts/Hallet model significantly decreases the expected viscous drag force imparted against an abrading clast below the transition size radius.

A decrease in the drag force due to premelting may have major implications for the capacity of entrained debris to resist glacier sliding or erode bedrock, especially if F _eff = 0. Hallet's model predicts relatively low shear stresses for typical bed loads and subglacial conditions (Hallet, 1979b, 1981; Cohen and others, 2005); however, as seen in Figure 2, debris-bed drag decreases significantly when the scaling factor is introduced. Furthermore, a ~3.5 ×  increase in the transition radius would limit the number of available bed tools capable of imparting significant force to the bed. For example, the transition radius for a granite clast increases from r _* ≈ 11 cm to r _* ≈ 38.5 cm. Given a typical size distribution, a large percentage of basal debris would fall below this size class (Haldorsen, 1981; Baltrūnas and others, 2009).

2. Methods

We simulated glacial abrasion with a direct shear device (Fig. 3a) by sliding a 10 cm × 10 cm limestone bed beneath a fixed slab of ice laden with 12 granitic rock fragments. To construct the ice slab, de-ionized (DI) water was first frozen, crushed and sieved to obtain a mixture of ice granules ranging in size from 2.0 to 6.3 mm. Encasing these granules in an ice matrix (Fig. 3b), we approximated a randomly oriented agglomeration of interlocking ice crystals, analogous to glacier ice. Granitic rock fragments were selected based on size and overall shape (Fig. 3e). Standard operating procedure dictates that particles must be sufficiently small to inhibit wall effects, which typically corresponds to a diameter one-tenth the width of the sample chamber (Iverson, 1990) and at least half its height (Gudmundsson, 1997; Zoet and Iverson, 2016). We chose angular clasts with a relatively high sphericity and mean diameters between 9 and 11 mm. The rock fragments were arrayed in a square insulating fiberglass sample chamber wrapped in foam tape; ice granules were layered on top; chilled DI water was added to fill in the void space; and the arrangement was frozen in a deep freezer. In the bottom half of the sample chamber, a smooth, flat limestone bed was emplaced, so that its top surface was flush with the adjacent edges and the overriding section can slide freely over it (Fig. 3c). The height of these beds varied slightly (generally less than ±0.5 mm) between experiments, due to their construction or as a result of sanding their surface between runs. If the height of the bed was lower than the walls of the sample chamber, paper spacers were placed beneath the limestone bed to make its edges flush with the adjacent walls. Five small holes (<4 mm) were drilled into the bed to facilitate drainage of meltwater during the experiment.

Fig. 3. (a) A schematic of the modified direct shear device, (b) debris-laden ice within the insulating sample chamber, (c) a typical limestone bed at the conclusion of an experiment (note that multiple striations are common for a single abrading rock), (d) striations, and (e) 12 representative granite clasts used in our experiments.

Once the sample chamber was assembled, it was loaded into the direct shear device. Vertical stress was applied to the top pressure plate by weights hung from the lever arm. The apparatus was then left to sit overnight, allowing the clasts to fully indent the bed and the ice to adjust under the load. To set the melt rate at the ice–bed interface (and by extension the vertical flow of ice toward the bed), we varied the power output to a heating plate installed in the bed. A glass bead thermistor, installed in the ice ~1 cm from the ice–bed interface, recorded the temperature of the ice as it warmed, and a linear variable displacement transducer measured the downward vertical displacement of the pressure plate as the ice melted at the interface. Because the heat flux was concentrated at the base of the ice and the ambient temperature of the freezer was held below the pressure melting point (PMP), we assume the total vertical displacement to approximately equal the thickness of the ice melted at the interface. Once the melt rate reached a relative steady state (>3 h), the shear ram was engaged at a prescribed velocity of 0.1 mm min⁻¹ for a distance of ~1 cm.

The relative influence of vertical ice velocity and effective stress on basal drag was assessed by conducting experiments at four normal stresses, σ _N, (four runs at 118 kPa, five runs at 294 kPa, five runs at 588 kPa and three runs at 784 kPa) over a range of melt rates. During the experiments conducted at 294, 588 and 784, the ambient freezer temperature was set to ~−0.9°C (±0.7°C). However, for the runs conducted at 118 kPa, an insulating case was built to house the direct shear in order to achieve smaller temperature gradients in the ice. The freezer temperature was set slightly colder at ~−1.5°C (±0.7°C) to compensate for heat retention within the insulating box. We set the melt rates between ~580 and 3000 mm a⁻¹, on the order (O(10² mm a⁻¹) to O(10³ mm a⁻¹)) expected for high rates of convergence expected on the stoss side of small bedrock obstacles (Hallet, 1979b; Iverson, 1990).

The shear stresses measured in our experiments represent friction between granite clasts and the bed and friction between clean ice and the bed. To assess the relative contribution of these two factors to the overall drag, we conducted experiments using only clean ice or only granite clasts, respectively. The background drag exerted by clean ice was determined by running experiments over a range of normal stresses, following the experimental protocol outlined above. We note that all σ _N = 118 kPa and σ _N = 784 kPa runs used the same limestone slabs as their respective background stress test, whereas the background stress for σ _N = 294 kPa and σ _N = 588 kPa was estimated using different beds. To determine the rock-on-rock frictional coefficient, 12 granite clasts were cast in epoxy resin and slid over a limestone bed at three normal stresses. Only the individual points of the indenting clasts contacted the bed in the rock-on-rock runs.

3. Results

During each experiment, shear stresses increased rapidly to a maximum before reaching a relative steady state or gradually declining for the duration of the run (Fig. 4). Fluctuations in shear stress occurred at timescales ranging from seconds to minutes, likely representing individual stick slip cycles or transient changes in bed-normal contact forces, F _c (Emerson and Rempel, 2007; Zoet and others, 2013). Similarly, the rate of vertical displacement of the top platen could vary on comparable timescales. Shear stresses presented in this study represent mean values calculated from 1500 s to the conclusion of the experiment when possible (Fig. 4a). Melt rates are calculated as the slope of an ordinary linear regression through vertical displacement data collected over the same time window (Fig. 4b). The number of independent data points used for each calculation varied between runs, but generally the sample size was >3500, which ensured small standard errors (Table 1).

Fig. 4. (a) Sample shear stress data for a representative experiment. The applied normal load was σ _N = 118 kPa. We calculate the average shear stress from 1500 s to the end of the experiment for most runs. The sampling window here consists of ~3800 unique datum. Corresponding displacement data for the top platen of the sample chamber is shown in (b). Melt rate is calculated as the average vertical rate of displacement over the same sampling window using an ordinary least squares regression.

Table 1. Data table showing the applied normal load (σ _N), basal melt rate (u _v), mean debris-bed shear stress (τ), the sample standard deviation of the measured shear stresses (STD τ), the vertical temperature gradients in the ice (dT/dz), the standard error of the mean for τ and dT/dz (SEM τ and SEM dT/dz, respectively) and the standard error of the model for u _v (SEM u _v)

N/A denotes a lack of data due to equipment malfunction.

Cavities formed sporadically in the lee of clasts in numerous experiments (0–3 per run; Fig. 5), but they were not a dominant feature and did not appear to affect the drag appreciably under our specific experimental conditions. Iverson (1990) reached a similar conclusion after noting comparable features in his experiments. However, this does not preclude the possibility that cavities may be more significant when they occur with greater frequency – merely that their impact on the force balance on a couple clasts did not dominate the overall drag. At higher melt rates, warmer ice temperatures or longer displacements, we anticipate these features would be more pervasive and could potentially have a greater impact on basal shear stress. In several runs at very high melt rates, we observed large cavities that appeared to have a more significant impact; however, we do not include them in our results as the experimental conditions during those runs are outside the scope of the present study.

Fig. 5. The morphology of a typical cavity that formed along the lee side of an abrading clast, cast in clay. In the experiments presented herein, similar cavities formed sporadically but did not significantly impact basal drag.

3.1 Friction tests

We obtained the frictional coefficients for the clean ice-on-limestone and granite-on-limestone tests by computing the slope of an ordinary linear regression through each dataset in σ _N–τ space. For both sets of experiments, shear stress increased linearly with increasing normal stress in a manner characteristic of Coulomb friction. For the granite-on-limestone runs, the frictional coefficient is μ _G−L = 0.72 (R ² = 0.99), which reflects the large hardness contrast between the two lithologies (Drewry, 1986, p. 52–53). For the clean ice-on-limestone experiments, the background coefficient of friction is taken to be μ _I−L = 0.14 (R ² = 0.86) for the runs conducted at σ _N = 294–784 kPa and μ _I−L = 0.19 (R ² = 0.99) for the runs at σ _N = 118 kPa. These values are high given that drag between clean, temperate ice and a level bed is typically assumed to be zero. However, a number of factors likely contributed to these elevated values: (1) friction between sample chamber walls, (2) edge effects, (3) micro-roughness elements on the rock surface and (4) the relatively short displacements in our experiments. Care was taken to ensure that the edges of the limestone bed were flush with the sample chamber. However, in some experiments, the bed remained slightly lower than the surrounding walls (<~0.5 mm), introducing a component of viscous drag to the resistive stresses as ice overrode this lip. Furthermore, as ice sat overnight under the load, it could creep into available void spaces around the edges of the slabs such that when the ram was engaged, ice in the void space would initially shear, adding another small viscous-drag component.

The maximum height of the micro-roughness elements on the surface of the limestone bed, imaged using a white-light interferometer, was on the order of ~10² µm. If the amplitude of these features exceeded the thickness of the lubricating film at the ice–bed interface, they would act as asperities to resist sliding. We estimate an upper bound for film thickness by calculating the thickness of a film separating temperate ice from an impermeable disc at the melt rates and normal stresses used in our experiments (Eqn (A7) in Emerson and Rempel (2007)). Film thicknesses in our experiments were likely thinner due to the permeability of the rock. We find that for all experimental conditions, the estimated thickness of the lubricating film (O(10⁻⁶–10⁻⁷ m) is likely smaller than the amplitude of the larger micro-roughness elements. Therefore, regelation past these features likely contributed to the measured background shear stress. Finally, we note that in the clean ice runs, the shear stress occasionally decreased from a maximum value toward a lower steady state, causing the calculated averages to potentially represent a value slightly higher than a true steady state. We attempt to account for the influence of all these mechanisms in all our experiments by subtracting a background drag from the observed shear stresses to isolate the drag, τ, imparted by the entrained clasts against the bed (i.e. τ = τ _observed − μ _I−Lσ _N).

3.2 Temperature gradients in the ice

As stated, experiments were conducted under two thermal regimes, in which (1) the sample chamber was exposed to the ambient freezer temperature (T _amb = −0.9 ± 0.7°C for σ _N = 294–784 kPa), or (2) the direct shear was housed in an insulating box (σ _N = 118 kPa) to minimize temperature gradients in the ice sample (σ _N = 118 kPa). Figure 6 shows how basal melt rate changes in response to the vertical temperature gradient in the ice, dT/dz, for runs with available temperature data. The thermal gradient reflects the change in temperature, dT, per unit distance, dz, from the PMP at a given σ _N according to the Clapeyron slope for air-free ice, C = 0.074 K MPa⁻¹. At the ice–bed interface, the ice temperature was at PMP in our experiments, as evidenced by the visible melt observed along the boundaries of the sample chamber and the continuous vertical displacement of the top pressure plate once the heating plate warmed. We calculate the gradient as the difference between PMP and the mean temperature recorded over the duration of the experiment at a known position in the ice divided by the distance of the thermistor from the interface. We achieved shallower thermal gradients in the ice by insulating the sample chamber. Equipment malfunctions with the thermistor prevented us from collecting temperature data during all experiments, but 65% of runs were successfully recorded. Regression models show that the melt rates recorded in our experiments scale linearly with dT/dz for both thermal regimes (Fig. 6).

Fig. 6. Vertical temperature gradients, dT/dz, observed in the ice slab scale linearly with basal melt rate, u _v. When the sample chamber was exposed to the ambient freezer temperature (σ _N = 294–784 kPa), the two parameters covary according to the relationship u _v = 11 194 dT/dz + 3018 (R ² = 0.82). With the addition of the insulating box (σ _N = 118 kPa), a linear regression model predicts u _v = 45 621 dT/dz + 2068 (R ² = 0.77). Standard errors for both parameters are smaller than the bounds of the plotted data points (Table 1).

3.3 Basal melt rate and effective stress influence the contact force

In all experiments, shear stress, τ, appears to scale linearly with basal melt rate, u _v (Fig. 7a), indicating the viscous drag force, F _d, is a control on debris-bed friction. For runs conducted at σ _N = 294–784 kPa (steeper dT/dz), data conform to a single linear trend, where τ = 0.057u _v + 35 (R ² = 0.85). However, the slope of the linear regression computed for the σ _N = 118 kPa data (shallower dT/dz) is notably smaller (τ = 0.016u _v + 20; R ² = 0.90). Table 1 lists our measurements and associated uncertainty. Because the calculated standard errors for u _v and τ are ~2–3 orders of magnitude smaller than the actual reported measurements, we assume their influence on the fit is negligible and therefore neglect to propagate these uncertainties in the regression.

Fig. 7. (a) Shear stress, τ, scales linearly with basal melt rate, u _v, but the rate of change is different between the two respective thermal regimes. (b) A statistically significant correlation is observed between the applied normal stress, σ _N, and τ as well, indicating that effective stress influences basal drag in our experiments. In (b), the σ _N = 118 kPa runs are excluded from the σ _N–τ regression, as they were conducted at different dT/dz. Standard errors for both u _v and τ are significantly smaller than the bounds of the plotted data points (Table 1) and are not included in the regression analysis.

Given that the contact force, F _c represents the contribution of both effective stress, F _eff, and the drag force, F _d (Eqn (2)), one anticipates that if F _eff is measurable, (1) the y-intercept of a regression in u _v–τ space at a given overburden will represent the contribution of F _eff to the observed shear stress (hereafter referred to as τ _eff) and (2) that this intercept should scale with increasing σ _N. This trend, however, is not obvious in our data. Table 2 lists the y-intercepts (given as τ _eff) calculated for linear regressions at each cryostatic stress. Though we cannot be certain of the statistical significance given the small sample sizes (n ⩽ 5), the data imply effective stress contributed meaningfully to basal drag. We also observe a weak relationship (τ = 0.14σ _N + 37; R ² = 0.49) between the applied normal load and basal drag (Fig. 7b), further indicating that effective stress influences debris-bed friction in our drained experiments. Scatter in this plot reflects the influence of basal melt rate for a given σ _N.

Table 2. Estimated contribution of effective stress to the observed debris-bed drag, τ _eff, based on the y-intercepts of linear regressions calculated for each applied normal stress σ _N in u _v–τ space

The values given for dτ/du _v are the slopes of the linear fit. We neglect the uncertainties for these estimates due to small sample sizes (n ⩽ 5).

4. Discussion

To frame our experimental results against theory, we first consider the viscosity and rheology of the ice used in our experiments. These parameters exert first-order controls on the viscous drag force (Eqns (3) and (4)) yet represent significant sources of uncertainty in the experiments. Initial treatments of abrasion assumed ice to be a linear viscous fluid with a power law exponent n = 1 (Watts, 1974; Hallet, 1979b, 1981; Shoemaker, 1988; Iverson, 1990). However, the deformation behavior of ice is variable depending on the stress state, timescale, temperature, water content and dominant deformation mechanism (Duval, 1979; Weertman, 1983; De La Chapelle and others, 1998; Goldsby and Kohlstedt, 2001). Therefore, it is instructive to explore how debris-bed contact forces vary with melt rate over a range of possible stress exponents (n = 1–3) and viscosity values.

Determining the creep exponent, n, and effective viscosity directly from the data is not possible due to uncertainty regarding the stress field around an abrading clast. To calculate effective viscosity, the effective shear stress acting in the immediate vicinity of a clast must be reasonably well-constrained (Iverson, 1990). However, this depends sensitively on the amount of melt produced along clast boundaries (see Eqn (A7)), which is unknown in the experiments. It is not obvious if the difference in basal drag observed between the two temperature regimes corresponds to a change in ice viscosity, a change in regelation rate or both. According to classic regelation and creep theory, regelation will not occur at temperatures below PMP. However, the presence of a premelt film may allow fluid flow to influence the drag force at subfreezing temperatures.

To address this, we compare our experimental results against the predictions of three different models using a range of possible viscosity values and creep exponents (n = 1–3). First, the drag force induced by viscous creep around a slippery sphere is calculated (Lliboutry and Ritz, 1978). In this formulation, slip is assumed to occur along clast boundaries (i.e. shear stresses are not supported on the surface). This model (henceforth referred to as ‘L & R’) serves as an analog for ice deforming around a suspended clast encased in a premelt film with negligible regelation. Treating ice as a power law fluid and neglecting regelation, the drag force, F _d, can be expressed

(7)$$F_{\rm d} = \; \alpha r^2\displaystyle{{\eta _{\rm i}} \over 2}\left( {\displaystyle{{u_{\rm v}} \over r}} \right)^{1/n},$$

where n is the creep exponent and α is a prefactor that varies with both n and particle geometry (Lliboutry and Ritz, 1978; Byers and others, 2012). When n = 1, Lliboutry and Ritz (1978) give the prefactor α = 4π and for n = 3, α = 8π.

Secondly, we implement the Watts/Hallet model (Watts, 1974; Hallet, 1979b, 1981) that accounts for the effects of both regelation and viscous creep on the drag force at PMP. As stated previously, this formulation (referred to as ‘W/H’) assumes boundary conditions typical of classic regelation and creep theory. To account for non-linear effective viscosities with the W/H, we use the extension derived by Watts (1974) and Cohen and others (2005; see Eqns (16–20)). Finally, we compare our experimental results against our modified Watts/Hallet solution (Section 1.1; Appendix) that scales the drag force to account for premelting effects at PMP. We denote this model as ‘PM’.

For each model, we calculate how τ, exerted by 12 idealized granite spheres (r = 5 mm) against a limestone bed varies with u _v. Shear stress is calculated as

(8)$$\tau = \mu D_{\rm c}\Phi F_{\rm d},$$

where D _c = 12 clasts per 0.01 m² is the areal debris concentration, μ = 0.72 is the coefficient of friction, and F _d is the drag force given by Eqns (3), (5) or (7), respectively. (Note the prefactor Φ in Eqn (3) is neglected since it is included in Eqn (8).) To calculate τ, the following nominal values are used: ρ _l = 1000 kg m⁻³, ρ _i = 917 kg m⁻³, ${\rm {\cal L}}$ = 3.3 × 10⁵ J kg⁻¹, C = 0.074 K MPa⁻¹, K = 2.5 W m⁻¹ K⁻¹. We also set the bed influence factor to Φ = 1.8 (Byers and others, 2012). For calculations involving a linear viscous rheology (n = 1), η _i = 5.9 × 10¹⁰ Pa s was used as a lower bound (Cohen, 2000) and 3.1 × 10¹² Pa s was used as the upper bound (Watts, 1974; Hallet, 1979b, 1981; Shoemaker, 1988; Iverson, 1990). For non-linear models, the viscosity values span an order of magnitude (i.e. O(10⁸ Pa s^1/2) and O(10⁷ Pa s^1/3) for n = 2 and n = 3, respectively), which generally encompasses the range of published values (review in Cohen, 2000).

4. 1 Model comparisons

Figures 8a–d compare the data collected at the colder temperature regime against the predictions of the L & R and W/H models with an assumed bed influence factor Φ = 1.8. For non-linear W/H, we consider how the form of the curve calculated for the upper and lower viscosity bounds fits the data, using the y-intercept as a free variable (Table 3). For non-linear L & R, we find the best fit for each exponent. W/H models require a y-intercept >0 to fit the data, whereas the best-fit non-linear L & R curves intercept the origin. For linear viscosities, we compare how the slope of the calculated u _v–τ relationship compares to the slope of a linear regression through the data: τ = 0.057u _v + 35 (R ² = 0.85; Table 4).

Fig. 8. A comparison of how shear stress, τ, varies with basal melt rate, u _v, in regimes dominated by viscous creep (a) or a combination of regelation and viscous creep (b–d) for different power law exponents (n = 1, 2, 3). Plot (a) displays best fit lines using Lliboutry and Ritz's (1978) model for viscous creep around a slippery sphere for n = 2 and n = 3, as well as a line calculated using a linear viscosity of 5.9 × 10¹⁰ Pa s. In (b–d), shaded areas represent the range of values calculated using the Watts/Hallet model (Watts, 1974; Hallet, 1979b) and an assumed bed influence factor Φ = 1.8. Lower and upper bounds correspond to the smallest and largest viscosity values used in the calculation. In plots (b) and (d), the darker lines represent curves calculated using preferred viscosity values published for similar experimental ice (Byers and others, 2012 for n = 1 and Zoet and Iverson, 2015 for n = 3, respectively). Black dashed lines in all plots represent a linear regression model (τ = 0.057u _v + 35; R ² = 0.85) through the larger dT/dz data (gray circles). The y-intercept is assumed to be 35 in plots (b–d) to compare the Watts/Hallet predictions with a best fit linear model.

Table 3. Goodness-of-fit statistics for the various models. We calculate a correlation coefficient, R ², and the root mean squared error of the residuals (RMSE) for each curve with respect to the data.

Lin. Reg. St refers to the linear regression models for the steeper dT/dz experiments. Models designated ‘lower’ or ‘upper’ refer to calculations using the lower and upper viscosity bounds.

Table 4. Calculated slopes in u _v–τ space and the corresponding linear viscosities, η _i, used for the various models

Lin. Reg. St. and Lin. Reg. Sh. refer to the linear regression models for the steeper dT/dz and the shallower dT/dz experiments, respectively. Models designated ‘lower’ or ‘upper’ refer to calculations using the lower and upper viscosity bounds. Slopes were calculated assuming a bed influence factor Φ = 1.8.

Estimating basal drag using non-linear exponents (Fig. 8a) produces a suboptimal fit compared to the linear model in our experiments (see Table 3 for goodness-of-fit statistics) indicating the ice in our experiments can be assumed to behave linearly. Although the magnitudes of observed τ lie within the range of the possible computed outcomes for the viscous creep model, the form of the L & R curves for n = 2 and n = 3 qualitatively do not conform to the data using realistic physical parameters (Watts, 1974; Lliboutry and Ritz, 1978; Hallet, 1979b, 1981; Cohen and others, 2005). Furthermore, using the lower bound for linear viscosity (B = 5.9 × 1¹⁰ Pa s; Cohen, 2000), the Newtonian case is still greater than the slope of the modelled regression by a factor of ~3.2 (Fig. 8a). Pronounced misfit of the viscous model with our observations suggests the drag force is not controlled solely by viscous creep.

Accounting for the effects of regelation in conjunction with viscous creep improves model fits overall. A stress exponent of n = 1 best approximates the observed relationship between u _v and τ (see Table 3). For linear ice viscosities using the W/H model, a bed influence factor of Φ ≈ 2.5–3 was found to match the slope of the modeled regression (Fig. 8b). These values closely align with Hallet's (1981) initial estimate of Φ = 2.4 and prior experimental results. Furthermore, using the viscosity reported by Byers and others (2012) for Newtonian ice results in a bed influence factor of Φ = 2.7 for our data. We favor this linear viscosity, as their ice was created using a similar protocol to ours and likely shared similar physical properties. However, it is important to note that the W/H model assumes fully temperate ice, whereas we observed vertical temperature gradients (O(−10⁻²°C cm⁻¹)) during all experiments with available temperature data (Fig. 6), indicating a direct comparison may not be appropriate.

For higher order stress exponents (n = 2–3) in the W/H model, the form of the curve transitions from a power law with exponent <1 to a straight line as ice viscosity increases (Figs 8c and d). This shift reflects a trade-off between the two competing flow mechanisms – viscous creep and regelation – on the overall drag force, as regelation and shear stress are linearly related (Eqn (A7); Figs 8c and d). The dominance of regelation is further illustrated in our modified version of the Watts/Hallet model, which incorporates premelting effects (Fig. 9a). For the physical parameters implemented in this calculation, the rate of regelation is amplified to the point where the influence of ice rheology on the drag force is effectively negated. For all three stress exponents, n = 1–3, the u _v–τ relationship for the premelt model can be approximated as a straight line with the slope of dτ/du _v ≈ 0.0034 kPa a mm⁻¹ – nearly ten times smaller than what the original W/H model would predict (Fig. 9a).

Fig. 9. Relationships between basal melt rate, u _v, and shear stress, τ. Figure (a) juxtaposes the data against the Watts solution with premelt assuming a y-intercept of 20. In (b), the blue-shaded area represents the bounds calculated for the Lliboutry and Ritz (1978) model with n = 1, using the viscosities in listed Table 4. Yellow bounds are the values calculated using the Watts/Hallet model with exponent n = 1 and red is n = 3. The green line is the modified Watts/Hallet model that accounts for premelting for all three stress exponents. Dashed lines are linear regressions through data collected at the two different dT/dz.

Because the curves in the W/H and PM model converge to a straight line for non-linear stress exponents at high viscosities, an n value cannot be confidently stated based on observed trends (Table 3; Figs 8 and 9). However, n = 1 is likely a fair approximation. For example, if n = 3 with a viscosity of B = 6.3 × 10⁷ Pa s^1/3 is chosen, as Zoet and Iverson (2015) determined for similarly constructed ice, the expected u _v–τ relationship can be estimated for comparison (Fig. 8d). Using these values and fitting the result to the data collected under the steeper dT/dz experiments also results in a nearly linear power law function (τ = 0.076u _v^0.91 + 60; R ² = 0.72). This leads us to infer that the ice in the experiment can be represented as a linear viscous fluid, in agreement with other ice deformation and subglacial contact-force experiments (Mellor and Testa, 1969; Dahl-Jensen and Gundestrup, 1987; Duval and Castelnau, 1995; Byers and others, 2012).

On timescales and stresses comparable to our experiments, ice crystals creeping past an entrained clast are exposed to a constantly changing stress field, giving rise to a transient creep regime in which ice behaves as a Newtonian fluid (Jacka, 1984; Cohen and others, 2005). In the context of prior abrasion/contact force experiments, Iverson (1990) observed a direct relationship between u _v–τ, and Byers and others (2012) reasoned a linearly proportional relationship between stress and strain rate for their ice. This lends credence to our interpretation that the y-intercepts of regressions through the data in u _v–τ space represent a contribution of effective stress to basal drag.

4.2 The effect of temperature gradients on debris-bed drag

The results at shallower dT/dz fall below the lower bound of the W/H curves and the upper bound of the PM solution (Fig. 9a). A linear u _v–τ relationship is likely for these experiments as well because they were conducted using the same ice as the colder runs, but with even lower deviatoric stresses and warmer temperatures, which would facilitate more melt along clast boundaries. Importantly, the slope of a linear regression through the data, dτ/du _v = 0.016 kPa a mm⁻¹ is ~2.1 times smaller than the lower bound of the W/H model for n = 1 (dτ/du _v = 0.033 kPa a mm⁻¹) and ~4.7 times larger than the PM version (dτ/du _v = 0.0034 kPa a mm⁻¹).

Melt along clast boundaries appear to influence basal drag in our experiments – even when ice temperatures are below PMP – since a purely viscous model overpredicts our data. We infer that the observed temperature gradients in the ice influence the u _v–τ relationship in part by changing the rate of fluid flow in the interfacial film surrounding the clast. This process likely happens in concert with temperature-dependent changes in viscosity. However, current models are unable to accurately capture the complex interplay between these two mechanisms at subfreezing temperatures. Basal melting occurred when heat flow to the ice–bed interface exceeded heat flow into the ice and out the boundaries. This flux produced larger temperature gradients in the sample than would occur in a thick, fully temperate section of ice beneath a glacier, but it might simulate local deviations in temperature within basal ice. Essentially, heat transfer along the boundary of the clast – whether produced by the heating plate, frictional heating, or refreezing – was sufficient to melt ice at a rate similar to the original W/H model for the experiments conducted at larger dT/dz. However, at smaller dT/dz, the melt rate along the clast boundaries exceeded the W/H model but not PM.

Although the W/H model fits the data well to the order of magnitude, it predicts slightly larger shear stresses than what we observe for experiments near PMP. This suggests that decreasing the vertical temperature gradient in the ice may further reduce τ for a given u _v. Potentially, the u _v–τ relationship may even converge toward the premelt solution when clasts are fully encased in temperate ice. However, achieving similar melt rates with fully temperate ice is beyond the capabilities of the experimental setup.

In the W/H and scaled PM models, the primary heat source along the clast boundary is assumed to be latent heat released when ice refreezes along zones of low pressure. However, the presence of a bed may inhibit refreezing and initiate cavity growth (Morris, 1979). Incomplete refreezing was observed to a degree in our experiments when cavities were present and in prior work (Iverson, 1990; Byers and others, 2012), but whether or not refreezing occurred in our experiments is not readily apparent. The fact that most clasts were fully encased in ice at the conclusion of each run speaks to its likelihood. However, even if refreezing is prevented in the vicinity of a clast, the process of regelation on a larger scale may still contribute meaningfully to the heat budget of entrained debris in a glacier. If debris-laden ice regelates around a bedrock obstacle, it releases heat on the lee side, which will diffuse through the bump into any abrading clasts in contact with the stoss face. Given that most abrasion is predicted to occur on the stoss side of bumps (Hallet, 1981; Shoemaker, 1988; Hallet and others, 1996), the process of regelation is likely an important source of heat for abrading clasts lodged in temperate ice.

4.3 Contribution of effective stress to debris-bed contact forces

The influence of effective stress on debris-bed contact forces relates to the presence of a lubricating film at the ice–bed interface, assumed to exist everywhere ice is in contact with the bed (Nye, 1973; Dash and others, 2006; Emerson and Rempel, 2007). If such a film exists and is hydrologically connected to the water film around the clast, a potential drainage pathway exists for melt produced along the clast's surface. Though Hallet's (1979b, 1981) model implicitly assumes that N = 0 in the film surrounding the clast, it does not take into account the presence of a far-field lubricating film. For this theory to be consistent, the pressure in both films must equal the cryostatic stress and remain constant in space and time. This assumption holds when the transmissivity at the bed is zero or if the ice/clast film is completely isolated from the far-field film. However, in drained conditions, such as our experiments, this assumption is likely not appropriate.

Although hydraulically isolated regions are common along glacial beds (Nye, 1973; Hallet, 1979a; Fountain, 1994; Engelhardt and Kamb, 1997; Andrews and others, 2014; Hoffman and others, 2016), fluctuating basal water pressure is routinely observed at hard-bedded glaciers (Fountain, 1994; Hubbard and others, 1995; Sugiyama and Gudmundsson, 2004). For example, Fountain (1994) observed that while water pressures in ~50% of the boreholes drilled at South Cascade glacier were at floatation, a subset was near atmospheric pressure in the vicinity of a conduit. If the hydraulic transmissivity of the bed is >0, as we designed it to be in our experiments, the magnitude of effective stress around the clast will then be controlled by (1) fluid potential in the far-field film, (2) the vertical stress components driving flow around the clast (i.e. N, F _d and F _b), (3) the relative thickness of the two films (Emerson and Rempel, 2007), and (4) the geometry of the clast.

The component of drag we attribute to effective stress in our experiments ranges from ~6% to 17% of the overburden (Table 2). For comparison, we calculate the expected drag due to effective stress, τ _eff, using Boulton's (1974) theory for 12 idealized spheres (Eqn (1)). (We neglect the viscous drag here to focus on τ _eff.) Shear stress is given

(9)$$\tau _{{\rm eff}} = \mu D_{\rm c}F_{{\rm eff}}\; $$

where D _c = 12 clasts per 0.01 m² is the areal debris concentration, and μ = 0.72 is the coefficient of friction determined empirically for granite against limestone. Even when considering an unrealistic, end-member scenario that has the thick far-field film (t = 1 mm) and negligible water pressure (P _l = 0), the maximum expected shear stress is ~2.4% of the ice load; for a more realistic film thickness on the order of microns, the expected drag is <1%. The shear stresses predicted by Boulton's theory are much smaller than we observe.

We propose two possible explanations for this discrepancy: (1) non-idealized clast geometry and (2) an uneven distribution of normal stress at the ice–bed interface. For a spherical particle, cryostatic stress is distributed uniformly across its surface due to symmetry, with vertical stresses balanced above and below the clast. A more angular or blocky clast, however, may have a greater contact area with the bed, and vertical stresses applied to the top surface will outweigh opposing, upward-facing stresses, leading to higher contact forces. Given that single clasts in our experiments commonly left multiple striations, this scenario is feasible (Figs 3c–e). In the context of the subglacial environment, it implies that (1) blocky clasts in transport along drained areas of a glacier's bed can impart significant force to the bed, and (2) the capacity for debris to impart force to the bed will decrease as particles get rounded in transport.

Another mechanism that could potentially increase the contribution of effective stress is an uneven distribution of normal stress at the ice–bed interface. Boulton's model (1974) assumes that pressure is evenly distributed across the bed due to the presence of the water film and that fluid pressure in that film bears a large portion of the normal load. If the thickness of this film was smaller than the dominant roughness elements along the bed's surface, stress would be transferred to the ice-on-rock contacts and the abrading clasts, which would act as isolated asperities (Zoet and Iverson, 2018). Decreasing the real area of contact corresponds to an increase in the normal force exerted on each clast to maintain the same shear stress, resulting in greater contact forces than Boulton's model predicts.

5. Conclusion

We conducted drained laboratory abrasion experiments to investigate the relative influence of basal melt rate and effective stress on subglacial rock friction in two thermal regimes. Basal shear stress scales linearly with both melt rate at the ice–bed interface and applied normal stress, and slight depressions in ice temperature increase the contribution of bed-normal viscous drag. We estimate the contribution of effective stress to basal drag to be between 6% and 17% of the overburden in our experiments, larger than predicted by Boulton's (1974) model for subglacial abrasion. These results suggest effective stress may be a dominant control on debris-bed contact forces in drained conditions for temperate ice, in contrast to the common assumption that it is negligible.

Comparing the observed relationships between melt rate and drag to the predictions of three friction models, we find the data are best described by a model that accounts for the effects of viscous creep and regelation on the bed-normal drag force. Observed shear stresses are lower than expected if viscous creep alone set the drag force, and the Watts/Hallet model (Watts, 1974; Hallet, 1979b, 1981) provides the best fit to our data overall, even for experiments conducted at temperatures below the PMP. We interpret this to mean that melt along clast boundaries influenced the contact force even at subfreezing temperatures due to the presence of an interfacial film. Incorporating premelting dynamics into the Watts/Hallet model (W/H) for debris-bed friction predicts nearly an order-of-magnitude decrease in shear stress and a corresponding ~3.5 ×  increase in the transition radius with respect to the original model. This model underpredicts our data by at least a factor of ~4, but it may provide a better fit for fully temperate ice.