Melting temperature changes during slip across subglacial cavities drive basal mass exchange

Abstract

The importance of glacier sliding has motivated a rich literature describing the thermomechanical interactions between ice, liquid water and bed materials. Early recognition of the gradient in melting temperature across small bed obstacles led to focused studies of regelation. An appreciation for the limits on ice deformation rates downstream of larger obstacles highlighted a role for cavitation, which has subsequently gained prominence in descriptions of subglacial drainage. Here, we show that the changes in melting temperature that accompany changes in normal stress along a sliding ice interface near cavities and other macroscopic drainage elements cause appreciable supercooling and basal mass exchange. This provides the basis of a novel formation mechanism for widely observed laminated debris-rich basal ice layers.

1. Introduction

At an ice–liquid interface, the dependence of melting temperature on normal stress drives ice regelation (e.g. Bottomley, 1872; Telford and Turner, 1963; Nye, 1967; Drake and Shreve, 1973; Gilpin, 1979; Rempel and Meyer, 2019). Regelation facilitates glacier sliding by causing ice to melt on the upstream sides of small bumps, where elevated normal stresses lower the melting temperature, and subsequently refreeze on their downstream sides, where reduced normal stresses increase the melting temperature. Larger bed obstacles are surmounted more easily by ice deformation, and the combination of processes is parameterized in well-known sliding laws (e.g. Weertman, 1957; Nye, 1969; Kamb, 1970). An important complication develops when sliding transports ice more rapidly than deformation enables it to conform to the shapes of rough beds. In these locations, cavities develop and fill with pressurized meltwater (e.g. Lliboutry, 1968; Schoof, 2005; Fowler, 2010). Links between cavities enable water flow, making them important for subglacial drainage (e.g. Kamb, 1987; Flowers, 2015). Cavity formation is also associated with heterogeneity in the normal stress exerted on the basal ice interface, since the liquid pressure is typically lower than the ice overburden pressure (e.g. Iken and Bindschadler, 1986) and higher stresses elsewhere must support the remaining glacier weight. Such elevated stresses can be generated by a combination of elevated liquid pressures and intermolecular interactions between mineral (i.e. rock or till) and ice surfaces across microscopic premelted liquid films (Dash and others, 2006; Rempel, 2008; Rempel and Meyer, 2019). Whereas conventional regelation sliding, described first by Weertman (1957), relies on the steady conductive heat transfer that accompanies melting temperature contrasts as ice slides across bedrock bumps, our attention is drawn to consider transient phase change processes promoted by conductive heat transfer between the basal interface and overlying ice.

Here, we explore the consequences of changes in melting temperature that are caused by stress heterogeneities along the basal interface. We focus on the freezing that is induced as ice slides from regions of elevated normal stress, supported in part by ice–mineral interactions, onto regions of reduced normal stress balanced solely by the liquid pressure. Such transitions are expected as slip transports ice over macroscopic drainage elements (e.g. cavities), and as macroscopic drainage elements incised upwards into the ice (e.g. R-channels) are dragged across newly unloaded mineral exposures. We find that a simple model for ice–liquid exchange in a sliding regime predicts freeze-on thicknesses that are consistent with diverse basal observations. Where this freeze-on occurs in the presence of unconsolidated or suspended sediments, debris bands in basal ice can form. We suggest that mm-scale diffuse and sometimes laminated debris bands that are preserved in basal ice layers (e.g. Sugden and others, 1987; Knight, 1997; Hubbard and others, 2009) may represent the signatures of the freeze-on processes that we describe. This is noteworthy because the properties of basal ice, including the presence and concentration of entrained debris, can influence sliding behavior and erosion (e.g. Thompson and others, 2020).

2. Basal phase behavior: bulk melting and premelting

For pure water, the equilibrium temperature at an ice–liquid interface is determined to leading order by the pressure in the liquid P and the normal stress exerted on the solid ice surface σ_n; under typical glaciological conditions, the other components of the ice stress tensor at the interface contribute negligibly. The thermodynamic arguments leading to this result are given by a number of sources (e.g. Kamb, 1961; Paterson, 1973; Sekerka and Cahn, 2004). The offset of the equilibrium temperature T _eq from a constant reference temperature T ₀ can be written as (c.f. Rempel, 2008, Eq. 1)

(1)$$T_{\rm eq}-T_0\approx -C_0\left[P-P_0 + {\rho_l\over \rho_{\rm l}-\rho_{\rm i}}\left(\sigma_{\rm n}-P\right)\right],\; $$

where P ₀ is the constant reference pressure (defined so that T _eq = T ₀ when P = σ_n = P ₀), ρ _l and ρ _i are the densities of liquid water and ice, and the magnitude of the Clapeyron slope is

(2)$$C_0 = T_0{\rho_{\rm l}-\rho_{\rm i}\over \rho_{\rm l}\rho_{\rm i} {\cal L}}\approx 7.4\times 10^{-8}\, {\rm K/Pa},\; $$

where ${\cal L}$ is the latent heat of fusion (a list of symbols is given in Table 1). An expanded discussion of equilibrium melting conditions is provided in the Supplementary Information.

Table 1. List of symbols, in alphabetical order, with Greek letters below (value of physical constant in parentheses).

Along an interface between ice and a macroscopic drainage element like a cavity or channel, the normal stress in the ice balances the liquid pressure (i.e. σ_n = P, see Fig. 1). Hence, along the ice-walled surfaces of macroscopic drainage elements, the final term in Eqn (1) vanishes and variations in the melting temperature are directly proportional to variations in liquid pressure, decreasing by approximately 7.4 mK for each atmosphere (~ 10⁵ Pa) increase in P. More generally, the intermolecular forces that cause premelted liquid films to separate ice from mineral surfaces (i.e. bedrock or unconsolidated sediments) produce a difference between the normal stress σ_n and the liquid pressure P that is referred to as the thermomolecular, or disjoining, pressure (e.g. Dash, 1989; Dash and others, 2006; Wettlaufer and Worster, 2006). A significant literature has developed to describe the molecular scale mechanisms that facilitate stress transfer between solid surfaces across intervening premelted films (e.g. Dash and others, 2006; Israelachvili, 2011), but for our purposes it suffices to note that all of these mechanisms produce the same essential result. When homogenized over an area element of the (generally rough) glacier bed, the thermomolecular pressure, defined here as P _T = σ_n − P, is the component of ice normal stress supported by intermolecular interactions with mineral surfaces (discussed further in the Supplementary Information). As expressed by Eqn (1), changes in T _eq with P _T at constant P are an order of magnitude larger (i.e. ρ _l/(ρ _l − ρ _i) ≈ 12) than changes in T _eq with changes in P of the same size during which P _T is held constant (i.e. along a path with σ_n = P; discussed further in the Supplementary Information).

Fig. 1. Schematic of a representative cavity, with uniform liquid temperature T _drainage = T _premelt + ΔT, connected to premelted films that separate the bounding rock/till interface from the overlying glacier ice. Representative temperature profiles are shown (left to right): (i) at the upstream side of the cavity prior to unloading (constant at T _premelt), (ii) midway over the cavity with a conductive profile promoting freeze-on (decaying from T _drainage to T _premelt far above the bed), and (iii) at the downstream side of the cavity after the interface temperature has returned to T _premelt and melting has commenced (note the dissipating thermal wave above the bed). Note that the melt rate at (iii) is lower than the freezing rate at (ii), promoting net freeze-on as ice flows across the schematic. Lower schematics show the force balance in the premelted film (left; P _T > 0) and in the cavity (right; P _T = 0).

Consider the idealized case of a glacier with ice-equivalent thickness H sliding over a heterogeneous bed that contains macroscopic drainage elements with P _T = 0 (Fig. 1). For basal areas ${\cal A}$ over which the glacier weight can be regarded as locally supported (i.e. shear stresses on the boundaries of ${\cal A}$ are assumed to contribute negligibly), the vertical force balance requires

(3)$$\int_{{\cal A}} \rho_{\rm i} g H \, {\rm d}{\cal A} = \int_{{\cal A}} \sigma_{\rm n} \, {\rm d}{\cal A} = \int_{{\cal A}} \left(P + P_{\rm T}\right)\, {\rm d}{\cal A}\; .$$

(As discussed further in the Supplementary Information, the basal surface itself is generally expected to be rough, and measurements of its physical area hence depend on the scale at which its topography is resolved; in the integral expressions used here for evaluating the vertical force balance, the area element ${\rm d}{\cal A}$ is rigorously defined as the horizontal projection of the contorted basal surface, making it independent of the scale of observation.) We assume that: (i) liquid wets the entire basal interface (consistent with our focus on sliding glaciers), (ii) spatial variations in liquid pressure are small (i.e. the bed is ‘well-drained’), (iii) macroscopic drainage elements occupy fractional bed area ϕ < 1, and (iv) the thermomolecular pressure is uniform over the remainder of the bed where ice–mineral separation diminishes enough that intermolecular forces between them are significant (what we refer to as premelted films throughout). For an average ice-equivalent thickness H ₀, under these idealized conditions Eqn (3) simplifies to (see the Supplementary Information)

(4)$$\rho_{\rm i} g H_0 = P + \left(1-\phi\right)P_{\rm T},\; $$

which indicates that, for basal areas ${\cal A}$ over which the glacier weight is locally supported, part of that weight is supported by the liquid pressure acting everywhere, and the remainder is supported by the thermomolecular pressure acting on the fraction of the bed that is not occupied by macroscopic drainage elements. Assigning the reference pressure as the average overburden P ₀ = ρ _igH ₀, Eqn (1) implies that the equilibrium melting temperature over the macroscopic drainage elements is

(5)$$T_{\rm drainage} = T_0 + C_0 N,\; $$

where the effective stress is defined as N = ρ _igH ₀ − P. Over the remainder of the bed, where stress transfer takes place between the ice and mineral surfaces (i.e. the ice–mineral separation is reduced to the microscopic dimensions of premelted liquid films), the equilibrium temperature is

(6)$$T_{\rm premelt} = T_{\rm drainage}- {C_0\over 1-\phi}{\rho_{\rm l}\over \rho_{\rm l}-\rho_{\rm i}}N.$$

Hence, the difference between the equilibrium ice–liquid temperature over a drainage element and that over premelted regions is

(7)$$\Delta T = {C_0\over 1-\phi} {\rho_{\rm l}\over \rho_{\rm l}-\rho_{\rm i}}N.$$

When ϕ is vanishingly small, ΔT increases by approximately 90 mK (i.e. 7.4 mK × 12) for each atmosphere (~ 10⁵ Pa) increase in the effective pressure N, or nearly 1 °C MPa⁻¹.

3. Basal mass exchange

3.1. Freeze-on over a single cavity

Having established that appreciable differences in equilibrium temperature can be expected across the borders of drainage elements, we next examine the effects of these differences on basal mass exchange. For a simple illustration, we consider a scenario in which ice slides over its bed at constant effective stress N and constant velocity u _s (see Fig. 1). Treating the heat flow as perpendicular to the bed, with thermal diffusivity κ, and assuming that phase changes are sufficiently slow that their effects on thermal conditions can be neglected (i.e. bed-perpendicular advective heat transport is negligible), energy conservation is governed by the one-dimensional heat equation with an interface temperature that changes abruptly by ΔT as the boundary of a drainage element is crossed at time t = 0. The resulting perturbation to the temperature gradient along the sliding interface can be approximated as (e.g. Carslaw and Jaeger, 1959, §2.5; discussed further in the Supplementary Information)

(8)$$\nabla T\approx -{ \Delta T\over \sqrt{\pi\kappa t}},\; $$

where the adopted sign convention signals that the heat flow upwards, away from the bed, increases when the equilibrium temperature across the drainage element boundary rises by ΔT. Noting that dx = u _sdt and integrating over a sliding distance $\ell = u_{s} t $, the implied perturbation to conductive heat transport equates with the latent heat of fusion for a layer of ice-equivalent thickness

(9)$$h_0 = \int_0^{\ell}-{C_{\rm p}\kappa \nabla T\over u_{\rm s}{\cal L}}\, {\rm d}x = {2\over \sqrt{\pi}}\ell {C_{\rm p} \Delta T\over {\cal L}}\sqrt{{\kappa\over u_{\rm s} \ell}},\; $$

where C _p is the specific heat capacity, and we note that the ratio of h ₀ to a characteristic cavity size $ \ell $ is inversely proportional to the product of the Stefan number $S_{\rm T} = {\cal L}/( C_{\rm p} \Delta T)$ with the square root of a Peclet number $ P_{\rm e} = u_{\rm s} \ell / \kappa $. The energy needed to accommodate this phase change is dissipated by conduction into the overlying ice, which retains the thermal signature of having recently been adjacent to colder premelted basal regions.

Equation (9) indicates that for a given cavity size ℓ, the thickness h ₀ is greater if the sliding speed is lower because there is more time for conductive heat transport. Freeze-on layer thickness is also greater if the effective stress is higher, which promotes elevated ΔT through Eqn (7). For intuition, at typical glacial sliding velocities u _s ranging between 10 and 10³ m a⁻¹, with ΔT between 0.1 and 1 °C (N between 10⁵ and 10⁶ Pa), the freeze-on layer grows to achieve ice-equivalent thicknesses between h ₀ ≈ 0.1 mm (large u _s, small ΔT and N) and h ₀ ≈ 10 mm (small u _s, large ΔT and N) during the time taken to traverse a cavity of dimension ℓ = 1 m. However, we note that the characteristic sizes ℓ of drainage elements need not be fixed, often increasing with sliding speed and decreasing with effective stress. A simple, illustrative model for cavity size that displays this qualitative behavior can be constructed using Glen's flow law with softness A ≈ 6.8 × 10⁻²⁴ s⁻¹Pa ⁻ⁿ (Cuffey and Paterson, 2010) to estimate a characteristic creep rate (e.g. Creyts and Schoof, 2009)

(10)$$v\approx A N^n \ell,\; $$

so that the distance slipped during creep closure of a cavity in the lee side of an obstacle of height d is

(11)$$\ell\approx {d u_{\rm s}\over v}\approx\sqrt{d u_{\rm s}\over AN^n}.$$

Substituting this and the undercooling expression from Eqn (7) into Eqn (9), while adopting a flow exponent of n = 3, gives (nominal values of physical constants are provided in Table 1)

(12)$$\eqalign{ h_0 & \approx {2C_{\rm p}C_0\rho_{\rm l}\over \left(1-\phi\right)\left(\rho_{\rm l}-\rho_{\rm i}\right){\cal L}}\sqrt{{\kappa\over \pi}}\left({N d\over u_{\rm s} A}\right)^{1/4}\cr & \approx {\left(4.2\times 10^{-6}\; {\rm m\, Pa^{-1/4}s^{-1/4}}\right)\over 1-\phi}\left({N d\over u_{\rm s} }\right)^{1/4},\; }$$

which is a relatively weak function of the primary variables that characterize the basal environment near the cavity, namely: N, d, and u _s. For example, with ϕ ≪ 1, N = 10⁵ Pa and d = 0.1 m, h ₀ ≈ 2 mm when u _s = 10 m a⁻¹, and this thickness drops only slightly to h ₀ ≈ 0.6 mm when u _s = 10³ m a⁻¹ (with N = 10⁵ Pa and d = 0.1 m) and increases slightly to h ₀ ≈ 6 mm when N = 10⁶ Pa and d = 1 m (with u _s = 10 m a⁻¹). We emphasize that these values of h ₀ should be regarded as order of magnitude estimates rather than precise predictions, particularly given the approximate treatment of the characteristic creep rate v in Eqn (10) and ℓ in Eqn (11).

3.2. Downstream melting

Downstream of the cavity, some of the recently frozen-on water will melt. Upon reloading the basal interface on the opposite boundary of the drainage element, the temperature gradient is perturbed once again. Assuming a symmetrical unloading/reloading cycle so that the interface temperature drops abruptly by ΔT and thereafter remains fixed for x > ℓ (see Fig. 1), the perturbation to the temperature gradient in the ice becomes (see Supplementary Information)

(13)$$\nabla T\approx -{ \Delta T\over \sqrt{\pi\kappa }}\left(\sqrt{{u_{\rm s}\over x}}-\sqrt{{u_{\rm s}\over x-\ell}}\right),\; $$

which induces gradual melting so that the net ice-equivalent freeze-on thickness evolves according to

(14)$$h = h_0\left(\sqrt{{x\over \ell}} - \sqrt{{x\over \ell}-1}\right)\approx {h_0\over 2}\sqrt{{\ell\over x}},\; $$

where the approximation on the right is valid for distances x ≫ ℓ. Importantly, even though the changes in interface temperature for this simple scenario are symmetrical – first increasing by ΔT at x = 0, then decreasing by ΔT at x = ℓ – the freezing and melting rates are not symmetrical. For example, the freeze-on thickness remains at $h = ( \sqrt {2}-1) h_0\approx 0.4h_0$ after sliding to x = 2ℓ – a distance equivalent to the drainage element dimension beyond its downstream boundary. This asymmetry in phase change behavior arises because conductive transport ensures that the attenuated history of past temperature perturbations continues to exert an influence on the changes in heat flux imparted by each new jump in interface temperature – essentially, the thermal pulse produced by unloading continues to modify the heat transport even after reloading returns the interface temperature to the background level T _premelt. Our treatment assumes an initial steady-state profile that reaches T _premelt at the basal interface (consistent with cases in which most of the bed is mantled by films so that ϕ is small), and subsequent perturbations to the temperature field forced by brief episodes with slightly warmer boundary temperatures result in net freeze-on as ice flows across and beyond drainage elements. It is worth noting that the supercooling described here is associated with motion of a cold sliding interface into contact with comparatively warmer water rather than the motion of comparatively colder water into contact with a warmer interface, as occurs during glaciohydraulic supercooling (e.g. Alley and others, 1998).

3.3. Idealized cavity sequences

A natural extension to this idealized treatment can be made by considering slip over cavities of dimension ℓ that are uniformly spaced by ℓ/ϕ, leading to a predicted freeze-on thickness immediately prior to the J + 1st unloading of

(15)$$h_{\rm min} \approx {h_0\sqrt{\phi}\over 2} \sum_{\,j = 1}^J j^{-1/2}\approx {2C_{\rm p} C_0\rho_{\rm l} N \phi\over {\cal L}\left(1-\phi\right)\left(\rho_{\rm l}-\rho_{\rm i}\right)}\sqrt{{\kappa x\over \pi u_{\rm s}}},\; $$

where the second approximation is valid for J ≫ 1, or equivalently x ≫ ℓ/ϕ. The next freeze-on episode enables the total thickness to reach

(16)$$h_{\rm max} \approx h_{\rm min} + h_0.$$

Figure 2 shows the evolution of minimum (dashed) and maximum (solid) freeze-on thicknesses predicted by Eqns (15) and (16) respectively, as a function of the sliding distance scaled to correspond with the number of complete unloading/reloading cycles. Net freeze-on increases gradually with the number of cycles and it also increases with the proportion of the bed occupied by macroscopic drainage elements ϕ. For example, with N = 1 MPa and ℓ = 1 m, the total freeze-on thickness reaches approximately 10 cm after a sliding distance of 10 km when ϕ = 0.1 and u _s ≈ 10 m a⁻¹ so that h ₀ ≈ 1 cm. Irrespective of cavity dimension ℓ, Eqn (15) indicates that h _min is approximately proportional to the product Nϕ and the square root of the total sliding duration x/u _s, so two orders of magnitude more time would be required to grow h _min to 1 m, and one order of magnitude lower Nϕ would reduce h _min to 1 cm. However, the same dynamic considerations that can cause ℓ to vary with basal conditions also make ϕ sensitive to N and u _s. For example, estimating the characteristic cavity size using Eqn (11) while treating the obstacle spacing as fixed suggests that the drainage fraction varies in the vicinity of some reference level ϕ ₀ according to

(17)$$\phi = \phi_0\sqrt{{u_{\rm s}/u_{{\rm s}0} \over N^n/N_0^n }},\; $$

where N ₀ and u _s0 are the reference effective stress and sliding velocity for which ϕ = ϕ ₀. Substituting this into Eqn (15) while taking n = 3 and assuming ϕ ≪ 1 leads to the functional behavior

(18)$$h_{\rm min}\approx {2C_{\rm p} C_0\rho_{\rm l} \phi_0\over {\cal L}\left(\rho_{\rm l}-\rho_{\rm i}\right)}\sqrt{{\kappa N_0^3 x\over \pi u_{{\rm s}0} N}}\approx \left(0.4 \, {\rm m^{1/2}\, Pa^{1/2}}\right)\phi_{\rm ref}\sqrt{{x\over N}},\; $$

where the numerical factor on the right is valid when the reference drainage fraction ϕ _ref = ϕ ₀ is defined using N ₀ = 10⁵ Pa and u _s0 = 10 m a⁻¹. The value of h _min predicted by Eqn (18) is notably independent of sliding speed u _s since related potential changes in the time available for freeze-on over drainage elements are negated by increases in the drainage element fraction, according to the simple treatment leading to Eqn (17). Moreover, this particular model for the controls on drainage fraction implies that, despite the weak direct dependence of h ₀ on effective stress in Eqn (12), h _min actually decreases gradually with increased N since the nonlinearity in creep closure rate, described with Glen's flow law exponent n = 3, causes changes in N to affect ϕ more strongly than their linear influence on ΔT.

Fig. 2. Evolution of predicted freeze-on thickness with sliding over an evenly spaced sequence of identical cavities. Sliding distance is scaled by the cavity spacing ϕ ⁻¹ℓ, with the values of ϕ noted in the legend. For a bed that contains more extensive cavities (high ϕ), there is proportionately less time for melt-out so h is larger for the same number of unload/reload cycles. Here, h is scaled by the characteristic dimension h ₀ from Eqn (9). Dashed lines depict h _min/h ₀, corresponding to predicted thicknesses on the upstream sides of cavities, while solid lines depict thicknesses h _max/h ₀ on the downstream sides of cavities. At other locations, h is expected to fall between these limits.

4. Discussion and conclusions

The simple cases we consider here, with P and P _T in Eqn (4) both constant, represent a considerable idealization. In reality, both of these fields are likely to be heterogeneous, as are the temperature variations ΔT that drive freeze-on and melt-out during each unloading/reloading cycle. Commonly observed temporal variations in basal effective stress further complicate the balance of freezing and melting at the bed (e.g. Iken and Bindschadler, 1986; Huss and others, 2007; Meierbachtol and others, 2013; Andrews and others, 2014; Rada and Schoof, 2018). All of these factors have the potential to leave their imprint on the basal ice and debris record. However, the relatively small variations that we predict for h ₀ and h _min in comparison with the much larger proportionate changes that we consider in characteristic basal parameters such as N and u _s, suggest that the dominant features of basal exchange may nevertheless be captured adequately by our idealized treatment. Extensions to consider more complicated and realistic scenarios with P and P _T both functions of time and space are not expected to introduce any significant conceptual hurdles.

In situations where unloading occurs, liquid water at T _drainage is invariably exposed to a supercooled surface, initially at T _premelt < T _drainage. Solidification in a dendritic growth habit may occur – whether downwards from the sliding glacial ice of a cavity roof (hinted at by reports of ‘regelation spicules’ and other similar features seen in cavities beneath relatively shallow ice, e.g. Kamb and LaChapelle, 1964; LaChapelle, 1968; Vivian and Bocquet, 1973; Theakstone, 1979; Rea and Whalley, 1994, see Fig. 3A), or upwards as anchor ice growing on the mineral surfaces of a cold channel floor. The liquid water that flows through macroscopic drainage elements is expected to often contain suspended debris and indeed is sometimes observed to be quite turbid (e.g. Rada and Schoof, 2018). These freezing styles may enhance the potential to incorporate mineral debris into the basal ice, particularly if mechanical disruption upon reloading (e.g. encapsulation in compacted dendritic layers or upwards particle displacement driven by debris–bed contact forces) acts to limit the efficiency of debris removal with small amounts of melt.

Fig. 3. Photos from the subglacial environment. (a) Dendritic ice crystal recovered in 1967 from the wall of a drained subglacial cavity in marginal ice beneath Blue Glacier, WA (reprinted with permission from LaChapelle, 1968). (b) Laminated facies photographed by KLR in 2014 beneath interior regions of Larsbreen glacier in Svalbard, accessed through the englacial drainage network and approximately 2.2 km from the glacier terminus.

The approximate perturbation to the temperature gradient described by Eqn (8) does not include the effects of background heat flow (discussed briefly in the Supplementary Information), and neither have we accounted for bed-normal deformation (e.g. Sugden and others, 1987; Knight, 1997). Thicker basal ice sequences may develop if the average rate of conductive transport away from the basal interface is greater than the combined effects of geothermal heating and frictional work so that net freeze-on (i.e. congelation) takes place. In such circumstances, the basal mass exchange caused by the unloading/reloading cycles described here are expected to add and then subtract from the total evolving thickness of basal ice layers. Other factors worthy of consideration include: (i) the enhanced thermal conductivity of most common bed materials in comparison with liquid water (tending to focus geothermal heating away from drainage elements and toward premelted regions), (ii) the absence of frictional work along the stress-free ice boundaries of macroscopic drainage elements (tending to favor slightly larger net freeze-on locally), (iii) advective transport of heat with moving meltwater (potentially counteracting or even overwhelming the reduced heat input associated with an absence of frictional work), (iv) longitudinal stress gradients that are often compressive in marginal regions (causing basal layers to thicken), and (v) the potential for intergranular phase changes in temperate ice to reduce conductive heat transport (discussed further in the Supplementary Information). Without accounting for such complications, the simplest scenario to consider is one in which the background rate of freeze-on is constant over the entire glacier bed. To this is added the supercooling-induced freeze-on upon basal unloading and interaction with potentially turbid cavity waters, followed by subsequent melt-out upon re-establishment of ice–mineral stress transmission, as outlined above. Different rates and characters of debris entrainment are expected to accompany slow freezing at the equilibrium melting temperature from premelted films of water over mineral surfaces in comparison with that entrained during more rapid freezing along a supercooled and potentially dendritic interface with macroscopic drainage elements. The effects of subtle differences in freezing environment beneath different portions of the glacial bed are difficult to quantify (e.g. drainage system turbidity, bedrock debris cover and particle size distribution, background freeze-on rate, spatially and temporally varying effective stress), though some effects may be amenable to laboratory investigation. These differences may be preserved as distinct mm-scale layering involving diffuse debris concentrations of a similar nature to those termed ‘clotted ice’, ‘dispersed facies’ or ‘laminated facies’, in reports of basal ice layers (e.g. Sugden and others, 1987; Knight, 1997; Hubbard and others, 2009). Further efforts to develop and test quantitative models that can reproduce specific features of observed basal debris distributions in particular settings hold promise for providing a window on the basal conditions sampled during flow.

The mechanism of debris entrainment introduced here is notably distinct from the ice-lensing phenomena that may sometimes accompany the growth of ‘frozen fringes’ with much higher (i.e. $> 50\percnt$ by volume) debris concentrations beneath soft-bedded glaciers (e.g. Rempel, 2008; Christoffersen and Tulaczyk, 2003; Meyer and others, 2018, 2019), or any of the other most commonly recognized debris entrainment mechanisms (e.g. Alley and other, 1997). As noted previously (e.g. Rempel, 2011), the thicknesses of sediment-rich entrained layers produced by frozen-fringe modeling depends on lensing criteria that can be related to a characteristic temperature offset, and for the low heat fluxes that are typical of subglacial environments this makes it very difficult to use that mechanism to explain mm-scale layering, such as that illustrated in Figure 3B. Particularly intriguing observations of entrained basal debris deep in the interior of Kamb Ice Stream are captured in borehole video accessible through Engelhardt (2013, https://nsidc.org/data/nsidc-0528/versions/1/documentation). In principal, entrainment by glaciohydraulic supercooling can result in diffuse debris contents and fine-scale layering, but this mechanism is normally invoked to describe behavior on adverse bed slopes (e.g. Alley and others, 1998), making it difficult to reconcile with the ubiquity of debris observations. While noting that several other potential mechanisms have been proposed and likely dominate in some environments, we suggest that remnants of the basal mass exchange processes introduced here may often be preserved as dispersed and laminated facies within basal ice layers.