1. Introduction
Friction between entrained basal debris and the bed affects the erosion rates (Schweizer and Iken, Reference Schweizer and Iken1992; Hallet and others, Reference Hallet, Hunter and Bogen1996; Koppes and Hallet, Reference Koppes and Hallet2006; Herman and others, Reference Herman2015; Koppes and others, Reference Koppes2015), sliding speed (Hallet, Reference Hallet1981; Shoemaker, Reference Shoemaker1988) and slip stability of glaciers (Zoet and others, Reference Zoet2013; McCarthy and others, Reference McCarthy, Savage and Nettles2017; Lipovsky and others, Reference Lipovsky2019). 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, Reference Lliboutry1968; Nye, Reference Nye1969; Iken, Reference Iken1981; Zoet and Iverson, Reference Zoet and Iverson2015, Reference Zoet and Iverson2016). However, observation shows that glaciers transport substantial bedload (Kirkbride, Reference Kirkbride and Menzies2002), and theory and field measurements indicate subglacial rock friction may be significant, albeit poorly constrained (Iverson and others, Reference Iverson2003; Cohen and others, Reference Cohen2005).
A central assumption in theories describing subglacial rock friction is that inclusions in ice are surrounded entirely by a thin water film (Boulton, Reference Boulton and Coates1974; Hallet, Reference Hallet1979b, Reference Hallet1981; Hindmarsh, Reference Hindmarsh1996; Cohen and others, Reference Cohen2005; Emerson and Rempel, Reference Emerson and Rempel2007). For temperate ice, classic regelation and creep theory predicts the presence of a film due to ice pressure melting along grain boundaries (Nye, Reference Nye1969, Reference Nye, Glen, Adie and Johnson1973; Kamb, Reference Kamb1970). However, considerations of premelting indicate a liquid film exists even at subfreezing temperatures (Gilpin, Reference Gilpin1979; Dash and others, Reference Dash, Fu and Wettlaufer1995; Rempel and Worster, Reference Rempel and Worster1999; Rempel and others, Reference Rempel, Wettlaufer and Worster2001; Rempel and others, Reference Rempel, Wettlaufer and Worster2004; Dash and others, Reference Dash, Rempel and Wettlaufer2006). 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 (Reference Boulton and Coates1974) and Hallet (Reference Hallet1979b, Reference Hallet1981) have driven theoretical framing of subglacial rock friction (Fig. 1). Boulton (Reference Boulton and Coates1974) 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
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 2), where r is the radius and t is the thickness of a water film that separates the ice from the bedrock. However, as Hallet (Reference Hallet1979b, Reference Hallet1981) 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:
Fig. 1. Conceptual models developed by (a) Boulton (Reference Boulton and Coates1974) and (b) Hallet (Reference Hallet1979b, Reference Hallet1981) 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 (Reference Hallet1979b, Reference Hallet1981) model. In Boulton's (Reference Boulton and Coates1974) 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 (Reference Byers, Cohen and Iverson2012).
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, Reference Hallet1979b). To quantify the drag force, Hallet adapted Watts's (Reference Watts1974) analytic solution for a spherical inclusion moving through temperate ‘Newtonian’ ice by regelation and creep to a half space, such that
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 (Reference Watts1974) defines the transition radius as
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 (Reference Hallet1981) based his initial estimate of Φ = 2.4 on the analysis of Goren (Reference Goren1970). However, the boundary conditions prescribed therein are not characteristic of the subglacial environment (Cohen and others, Reference Cohen2005), and empirical evidence regarding the magnitude of this parameter is limited and often contradictory. For example, Iverson (Reference Iverson1990) observed Φ = 0–5 for non-idealized granite clasts during laboratory abrasion experiments. Byers and others (Reference Byers, Cohen and Iverson2012) 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 (Reference Lee and Rutter2004) 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 (Reference Hallet1979b, Reference Hallet1981) has been widely implemented in glacier erosion models (Herman and others, Reference Herman, Beaud, Champagnac, Lemieux and Sternai2011; Anderson, Reference Anderson2014; Beaud and others, Reference Beaud, Flowers and Pimentel2014; Egholm and others, Reference Egholm2017; Ugelvig and Egholm, Reference Ugelvig and Egholm2018), and available experimental data support a linear relationship between F c and u v as predicted (Iverson, Reference Iverson1990; Byers and others, Reference Byers, Cohen and Iverson2012). Although Hallet's model identified key controls on debris-bed friction, numerous lines of evidence indicate the problem remains unresolved. As Hindmarsh (Reference Hindmarsh1996) and Morris (Reference Morris1979) address in considerable detail, mathematical inconsistencies arise when Watts's (Reference Watts1974) 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 (Reference Emerson and Rempel2007) 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, Reference Iverson2003; Cohen and others, Reference Cohen2005). 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 (Reference Rempel and Meyer2019) 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' (Reference Watts1974) 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,
where
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 (Reference Watts1974) original analysis (Eqn (4)). Figure 2 juxtaposes the relationship between clast radius and drag force predicted for both solutions, assuming η i = 3.15 × 1012 Pa s, u v = 500 mm a−1, ρ l = 1000 kg m−3, ρ i = 917 kg m−3, ${\rm {\cal L}}$
= 3.3 × 105 J kg−1, C = 0.074 K MPa−1, K = 2.5 W m−1 K−1.
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, Reference Hallet1979b, Reference Hallet1981; Cohen and others, Reference Cohen2005); 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, Reference Haldorsen1981; Baltrūnas and others, Reference Baltrūnas, Šinkūnas, Karmaza, Česnulevičius and Činkūnė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, Reference Iverson1990) and at least half its height (Gudmundsson, Reference Gudmundsson1997; Zoet and Iverson, Reference Zoet and Iverson2016). 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−1 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−1, on the order (O(102 mm a−1) to O(103 mm a−1)) expected for high rates of convergence expected on the stoss side of small bedrock obstacles (Hallet, Reference Hallet1979b; Iverson, Reference Iverson1990).
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, Reference Emerson and Rempel2007; Zoet and others, Reference Zoet2013). 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 (Reference Iverson1990) 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 2 = 0.99), which reflects the large hardness contrast between the two lithologies (Drewry, Reference Drewry1986, p. 52–53). For the clean ice-on-limestone experiments, the background coefficient of friction is taken to be μ I−L = 0.14 (R 2 = 0.86) for the runs conducted at σ N = 294–784 kPa and μ I−L = 0.19 (R 2 = 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 ~102 µ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 (Reference Emerson and Rempel2007)). 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−6–10−7 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−1. 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 2 = 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 2 = 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 2 = 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 2 = 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 2 = 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, Reference Watts1974; Hallet, Reference Hallet1979b, Reference Hallet1981; Shoemaker, Reference Shoemaker1988; Iverson, Reference Iverson1990). However, the deformation behavior of ice is variable depending on the stress state, timescale, temperature, water content and dominant deformation mechanism (Duval, Reference Duval1979; Weertman, Reference Weertman1983; De La Chapelle and others, Reference De La Chapelle, Castelnau, Lipenkov and Duval1998; Goldsby and Kohlstedt, Reference Goldsby and Kohlstedt2001). 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, Reference Iverson1990). 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, Reference Lliboutry and Ritz1978). 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
where n is the creep exponent and α is a prefactor that varies with both n and particle geometry (Lliboutry and Ritz, Reference Lliboutry and Ritz1978; Byers and others, Reference Byers, Cohen and Iverson2012). When n = 1, Lliboutry and Ritz (Reference Lliboutry and Ritz1978) give the prefactor α = 4π and for n = 3, α = 8π.
Secondly, we implement the Watts/Hallet model (Watts, Reference Watts1974; Hallet, Reference Hallet1979b, Reference Hallet1981) 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 (Reference Watts1974) and Cohen and others (Reference Cohen2005; 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
where D c = 12 clasts per 0.01 m2 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−3, ρ i = 917 kg m−3, ${\rm {\cal L}}$
= 3.3 × 105 J kg−1, C = 0.074 K MPa−1, K = 2.5 W m−1 K−1. We also set the bed influence factor to Φ = 1.8 (Byers and others, Reference Byers, Cohen and Iverson2012). For calculations involving a linear viscous rheology (n = 1), η i = 5.9 × 1010 Pa s was used as a lower bound (Cohen, Reference Cohen2000) and 3.1 × 1012 Pa s was used as the upper bound (Watts, Reference Watts1974; Hallet, Reference Hallet1979b, Reference Hallet1981; Shoemaker, Reference Shoemaker1988; Iverson, Reference Iverson1990). For non-linear models, the viscosity values span an order of magnitude (i.e. O(108 Pa s1/2) and O(107 Pa s1/3) for n = 2 and n = 3, respectively), which generally encompasses the range of published values (review in Cohen, Reference Cohen2000).
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 2 = 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 (Reference Lliboutry and Ritz1978) 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 × 1010 Pa s. In (b–d), shaded areas represent the range of values calculated using the Watts/Hallet model (Watts, Reference Watts1974; Hallet, Reference Hallet1979b) 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, Reference Byers, Cohen and Iverson2012 for n = 1 and Zoet and Iverson, Reference Zoet and Iverson2015 for n = 3, respectively). Black dashed lines in all plots represent a linear regression model (τ = 0.057u v + 35; R 2 = 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 2, 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, Reference Watts1974; Lliboutry and Ritz, Reference Lliboutry and Ritz1978; Hallet, Reference Hallet1979b, Reference Hallet1981; Cohen and others, Reference Cohen2005). Furthermore, using the lower bound for linear viscosity (B = 5.9 × 110 Pa s; Cohen, Reference Cohen2000), 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 (Reference Hallet1981) initial estimate of Φ = 2.4 and prior experimental results. Furthermore, using the viscosity reported by Byers and others (Reference Byers, Cohen and Iverson2012) 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−2°C cm−1)) 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−1 – 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 (Reference Lliboutry and Ritz1978) 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 × 107 Pa s1/3 is chosen, as Zoet and Iverson (Reference Zoet and Iverson2015) 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 v0.91 + 60; R 2 = 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, Reference Mellor and Testa1969; Dahl-Jensen and Gundestrup, Reference Dahl-Jensen and Gundestrup1987; Duval and Castelnau, Reference Duval and Castelnau1995; Byers and others, Reference Byers, Cohen and Iverson2012).
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, Reference Jacka1984; Cohen and others, Reference Cohen2005). In the context of prior abrasion/contact force experiments, Iverson (Reference Iverson1990) observed a direct relationship between u v–τ, and Byers and others (Reference Byers, Cohen and Iverson2012) 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−1 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−1) and ~4.7 times larger than the PM version (dτ/du v = 0.0034 kPa a mm−1).
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, Reference Morris1979). Incomplete refreezing was observed to a degree in our experiments when cavities were present and in prior work (Iverson, Reference Iverson1990; Byers and others, Reference Byers, Cohen and Iverson2012), 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, Reference Hallet1981; Shoemaker, Reference Shoemaker1988; Hallet and others, Reference Hallet, Hunter and Bogen1996), 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, Reference Nye, Glen, Adie and Johnson1973; Dash and others, Reference Dash, Rempel and Wettlaufer2006; Emerson and Rempel, Reference Emerson and Rempel2007). 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 (Reference Hallet1979b, Reference Hallet1981) 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, Reference Nye, Glen, Adie and Johnson1973; Hallet, Reference Hallet1979a; Fountain, Reference Fountain1994; Engelhardt and Kamb, Reference Engelhardt and Kamb1997; Andrews and others, Reference Andrews2014; Hoffman and others, Reference Hoffman2016), fluctuating basal water pressure is routinely observed at hard-bedded glaciers (Fountain, Reference Fountain1994; Hubbard and others, Reference Hubbard and Sharp1995; Sugiyama and Gudmundsson, Reference Sugiyama and Gudmundsson2004). For example, Fountain (Reference Fountain1994) 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, Reference Emerson and Rempel2007), 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 (Reference Boulton and Coates1974) theory for 12 idealized spheres (Eqn (1)). (We neglect the viscous drag here to focus on τ eff.) Shear stress is given
where D c = 12 clasts per 0.01 m2 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 (Reference Boulton and Coates1974) 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, Reference Zoet and Iverson2018). 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 (Reference Boulton and Coates1974) 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, Reference Watts1974; Hallet, Reference Hallet1979b, Reference Hallet1981) 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.
Acknowledgements
The authors thank the reviewers for their enlightening comments and the following people for their valuable insight and conversation: N. Lord, P. Sobol, C. Meyer, N. Iverson, N. Stevens, J. Woodard, and K. Feigl. N. Lord and P. Sobol were instrumental in the success of this project, constructing key pieces of our experimental setup and providing invaluable technical support and advice. E. Parrish provided the photograph of the clast. D. D. Hansen was supported by funding from Wisconsin Alumni Research Association awarded to L. K. Zoet.
Appendix
Premelting effects on debris-bed contact forces
Closely following the rationale outlined in Rempel and Meyer (Reference Rempel and Meyer2019), we incorporate the same premelting condition into Watts's (Reference Watts1974) original analysis of temperate ice moving past a sphere through a combination of regelation and viscous creep. For brevity's sake, we touch on the major points of their argument and refer the reader to Eqns (1–9) in Rempel and Meyer (Reference Rempel and Meyer2019) for an in-depth discussion of their justifications. Equations from Rempel and Meyer (Reference Rempel and Meyer2019) are referenced with ‘RM’ preceding the number; equations from Watts are preceded by ‘W’.
Watts (Reference Watts1974) considers the case of a spherical inclusion forced through an infinite body of temperate ice, in which motion is accommodated through a combination of regelation and viscous creep. A thin film of water surrounds the clast, and the temperature at the ice–liquid interface is approximated as
where ΔT = T − T m, T m is the bulk melting temperature at ambient pressure, T is the temperature at the ice–liquid interface, $\rho _{\rm i}{\rm {\cal L}\;}$
is the volumetric latent heat of fusion for ice, K is the mean thermal conductivity of ice and rock, r is the radius of the sphere, and v r is the regelation rate (see W15−16). Normal stresses imparted against the clast through the water film are assumed to covary with temperature (see W16), such that
where C is the change in melting temperature with pressure, and Δσ rr is the corresponding change in radial stress in polar coordinates relative to ambient hydrostatic pressure, σ amb or P amb (i.e. Δσ rr = σ rr – σ amb). Though Watts (Reference Watts1974) expresses the radial stress in polar coordinates, Cauchy's law shows that Δσ rr is equivalent to the change in water pressure in the film, ΔP l = P l – P amb, which is assumed to equal the stresses imparted by the ice, ΔP i = P i – P amb. Therefore, by eliminating ΔT from equating Eqns (6) and (7), we obtain an expression for v r in terms of the pressure in the film (see W17), whereupon
Water movement through the thin film surrounding the clast can be approximated as Poiseuille flow between two parallel plates. Likewise, conservation of mass dictates that the pressure gradient driving fluid flow through the film is a function of the regelation rate (see RM5), expressed as
where δ is the thickness of the film surrounding the clast, and η w is the viscosity of liquid water. Combining Eqn (A3) with (A4) and solving for δ using textbook values for the various physical constants (see RM5-RM6 in Rempel and Meyer (Reference Rempel and Meyer2019) and parameter values listed therein) implies a film thickness of ~3.1 nm for the case where ΔP i = ΔP l. This result essentially forms the crux of Rempel and Meyer's (Reference Rempel and Meyer2019) argument: (1) A nm-scale film thickness is unlikely to occur in nature; (2) on that scale, premelting dynamics enable significant changes in effective stress within the film (i.e. ΔN ≠ 0); (3) therefore, regelation by pressure melting alone is unlikely for the case of water flowing through a thin film bounded by impermeable solids. Instead, it is necessary to invoke the generalized Clapeyron equation to describe temperature changes with respect to both ΔP l and ΔN (Rempel, Reference Rempel2008; RM7 in Rempel and Meyer, Reference Rempel and Meyer2019), in which
and ΔN = ΔP i − ΔP l.
Substituting the expression for ΔT in Eqn (A5) into (A1) and considering the typical case where ΔP i >> ΔP l (see RM8 and RM9) leads to
This expression for the regelation rate replaces W5 in Watts's (Reference Watts1974) analysis. Watts gives the radial pressure imparted by the ice through the film (see W18) as
where η i is the viscosity of ice. Following Watts's (Reference Watts1974) steps, we substitute the right side of Eqn (A7) into (A6) and solve for v r to find
Finally, by replacing v r in Eqn (A7) with the equivalent expression in Eqn (A8) and then integrating across the sphere's surface, we obtain the drag force, F d, imparted against the sphere, given as
where
Notably, we find that 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 (Reference Watts1974) original analysis
Open access
