2. Clean ice

In sliding theories, no shear traction acts locally at the bed surface, owing to the water film that divides wet-based glaciers from their beds. Ideally, a local shear traction commensurate with debris-bed friction could be included in such models so that stresses, ice flow and cavity sizes could reflect that friction. Although aspects of this problem have been studied with the finite-element method (Schweizer and Iken, Reference Schweizer and Iken1992), no analytical model exists that includes local shear tractions and ice–bed separation. Therefore, our tactic is to investigate what insights are possible using existing sliding theory to estimate debris-bed friction, neglecting its feedback on ice flow – a reasonable assumption if shear stress from debris-bed friction is a small fraction of the total basal shear stress.

To analyze sliding with ice–bed separation, we use the theory of Lliboutry (Reference Lliboutry1968) because it enables, in an approximate way, use of a Glen-type, non-linear ice rheology, and because aspects of the theory have been tested experimentally (Zoet and Iverson, Reference Zoet and Iverson2015). Like most other sliding theories that include ice–bed separation, the theory neglects regelation. This simplification is justifiable because on stoss surfaces where ice is in contact with rock, glacial abrasion generally smooths surfaces (e.g. Hooke, Reference Hooke2005), causing roughness elements as small as the controlling obstacle size (e.g. Weertman, Reference Weertman1957) to be relatively rare. Slip velocity is imposed as an independent variable in the far-field above bumps. No shear traction is allowed locally on the bed surface, owing to the water film, as noted, and the shear traction between water in cavities and cavity roofs is negligibly small. All drag therefore, for the case of clean ice, stems from the distribution of normal stress on the glacier sole exerted by rock in zones of ice–bed contact and by water under pressure in lee side cavities.

Lliboutry (Reference Lliboutry1968) considered a 2-D, rigid bed with a longitudinal profile,

(1)$$z_{\rm b} = \displaystyle{1 \over 2}a\cos \displaystyle{{2\pi x} \over \lambda}, $$

where x is measured parallel to average slope of the bed, a is the bump height from trough to crest, and λ is the wavelength (Fig. 1a). The roughness is ${\rm {\cal R}}$ = a/λ. Making various approximations, Lliboutry used the bump geometry, slip velocity u, effective pressure, N and the ice flow law parameters to estimate the cavity length, L _c (Fig. 2):

(2)$$L_{\rm c} = \lambda (1-s),$$

where s is the fraction of the bed in contact with ice:

(3)$$s = \displaystyle{1 \over {6\pi {\rm {\cal R}}}}\left( {\displaystyle{{aAN^n} \over u}} \right)^{1/2}.$$

Fig. 1. (a) Model parameters and coordinates for the wavy bed of Lliboutry (Reference Lliboutry1968). (b) Model parameters and coordinates for a particle in frictional contact with the bed over the zone of ice–bed contact, sλ, in panel (a).

Fig. 2. Cavity length, L _c, scaled by the bed wavelength, λ, as a function of slip velocity scaled by the velocity at which ice separates from lee surfaces, as indicated by the theories of Lliboutry (Reference Lliboutry1968) and Kamb (Reference Kamb1987, see his Eqns (8) and (13)).

The flow law of ice is represented by the pre-factor, A, and the creep exponent n (Cuffey and Paterson, Reference Cuffey and Paterson2010, p. 55). Effective pressure N = ρgH cos α − p _w, where ρ is the ice density, g is the gravitational acceleration, H is the ice thickness, α is the glacier surface slope and p _w is the cavity water pressure. Neglected in deriving Eqn (3) is heat dissipated by water flowing through cavities and resultant melting of cavity roofs; such dissipation is thought to be small in cavities (Walder, Reference Walder1986). Confidence in this estimate of cavity size comes from Kamb (Reference Kamb1987), who considered the same bed geometry and ice rheology and arrived at a similar estimate using a different analysis (Fig. 2).

Also required for estimating basal shear stress is the location, x _d, where the ice detaches from the bed:

(4)$$x_{\rm d} = \displaystyle{\lambda \over {2\pi}} \cot ^{-1}\left[ {\displaystyle{{2\pi \lpar {1-s} \rpar + \; \sin \lpar {2\pi s} \rpar } \over {1-\cos \lpar {2\pi s} \rpar }}} \right]$$

(Lliboutry, Reference Lliboutry1968), where x _d = 0 would correspond to ice detachment at the bump crest (Fig. 1a). For increasing sliding speeds, the point of detachment is increasingly close to the bump crest upstream, and the reattachment point, x _r = x _d + L _c (Fig. 1a) approaches the downstream crest (Fig. 3). These results are qualitatively similar to those of Fowler (Reference Fowler1986) who considered a linear ice rheology, although in his analysis x _d extends slightly up-glacier beyond the bump crest, and the ice reattachment point extends farther down-glacier for a given slip velocity scaled by N and ${\rm {\cal R}}$.

Fig. 3. Points of detachment, x _d, and reattachment, x _r, of ice, as a function of slip velocity scaled by the velocity at which ice separates from lee surfaces, indicated by the theories of Lliboutry (Reference Lliboutry1968) and Kamb (Reference Kamb1987). Bump crests are at x/λ = 0 and x/λ = 1.0.

The basal drag depends on the normal stress, σ _n, that ice exerts on the bed over zones of ice–bed contact. Lliboutry (Reference Lliboutry1968) assumed that this stress varies symmetrically about the centers of such zones where the stress reaches a maximum, to arrive at

(5)$${\; \sigma _{\rm n} = p_{\rm w}-\; \displaystyle{{\pi \; N} \over {\sin \lpar {\pi s} \rpar -\pi s\cos \lpar {\pi s} \rpar }}\left\{ {\cos \lpar {\pi s} \rpar + \cos \left[ {\displaystyle{{2\pi} \over \lambda} \left( {x-\displaystyle{{x_{\rm d} + x_{\rm r}} \over 2}} \right)} \right]} \right\}}$$

(see his Eqns (10 and 12)). From this stress distribution, he estimated the basal shear stress, τ, for clean ice:

(6)$$\tau = \displaystyle{1 \over 2}\pi {\rm {\cal R}}N\left\{ {\displaystyle{{\lsqb {\pi s-(1/2)\sin \lpar {2\pi s} \rpar } \rsqb \sin \lpar {\pi s-(2\pi x_{\rm d}/\lambda )} \rpar } \over {\sin \lpar {\pi s} \rpar -\pi s\cos (\pi s)}}} \right\}$$

(Lliboutry, Reference Lliboutry1968, his Eqns (13 and 14)), where s depends on sliding speed through Eqn (3). The resultant rate-weakening drag (Fig. 4) is qualitatively similar to that determined in subsequent analyses (Fowler, Reference Fowler1986; Schoof, Reference Schoof2005; Gagliardini and others, Reference Gagliardini, Cohen, Råback and Zwinger2007). A drawback of Lliboutry's model is that he derived separate relationships for sliding velocities below and above the velocity value at which ice–bed separation occurs, resulting in a non-smooth first derivative of the stress relationship at the separation velocity. Thus, Figures 2–4, with sliding velocity scaled by the separation velocity, do not include values <1.0. This scaling is acceptable because we are concerned with the effect of debris on rate-weakening drag that accompanies cavity growth. The separation velocity can be found from Eqn (3) by setting s = 1.

Fig. 4. Sliding rule of Lliboutry (Reference Lliboutry1968) for sliding velocities above those required for ice–bed separation. Shear stress, τ, is scaled following Gagliardini and others (Reference Gagliardini, Cohen, Råback and Zwinger2007) where $C = \pi {\rm {\cal R}\;}$ is the maximum slope of the bed. The factor, 1.51, adjusts τ to bring it into accord with Iken's (Reference Iken1981) bound. Solving for τ in each of her Eqns (3) and (4) and dividing the results yields this factor; it reflects different assumptions made by Lliboutry (Reference Lliboutry1968) and Iken (Reference Iken1981) regarding the distribution of normal stress on zones of ice–bed contact.

3. Debris-bed friction

Lliboutry's (Reference Lliboutry1968) sliding model allows estimation of rates of ice–bed convergence with stoss surfaces and cavity size as functions of slip velocity. The former controls contact forces between particles and the bed, whereas the latter controls the area of the bed over which debris-bed friction occurs. We wish to explore how the tradeoff between these two variables for different cavity sizes affects the form of the relationship between shear stress and slip velocity.

A major component of ice–bed convergence on stoss surfaces is from extension of ice parallel to the bed as ice flows past bumps (Hallet, Reference Hallet1981). Sliding theory for a sinusoidal bed, in the absence of ice–bed separation and neglecting regelation, yields a component of ice velocity toward the bed everywhere on stoss surfaces at rates that peak at the point midway up stoss surfaces (Nye, Reference Nye1969, see his Eqn (32) in the limit of zero regelation past bumps). Iken's (Reference Iken1981, her Fig. 6a) numerical calculations of the velocity field in Newtonian ice near a sinusoidal bed for the case of steady cavities similarly demonstrate that a streamline near stoss surfaces converges with the bed and that ice speeds up as it moves down-glacier along zones of ice–bed contact, indicative of bed-parallel extension on stoss surfaces. These facts are also confirmed by recent finite-element simulations of the same process with a Glen-type ice rheology (Helanow and others, Reference Helanow, Iverson and Zoet2018). The ubiquity of striations on stoss surfaces (e.g. Hallet, Reference Hallet1979) of large bumps where leeward cavities have been mapped (e.g. Hallet and Anderson, Reference Hallet and Anderson1980) provides further confirmation that extension of ice parallel to the bed during flow past bumps presses debris particles against stoss surfaces – a process that was first recognized by Gilbert (Reference Gilbert1906).

This ice extension causes a component of ice velocity toward the bed that is zero at the bed surface but increases away from the bed in the zone of extension. To estimate the bed-normal ice velocity that results from this ice deformation, v _d, we consider the more compressive principal stress, σ ₁₁, to be normal to areas of ice–bed contact, following Lliboutry (Reference Lliboutry1968), so σ ₁₁ = σ _n (Fig. 1a). A reasonable choice for the less compressive principal stress is σ ₂₂ = p _w, the water pressure in cavities, which confines the ice on each side of the zone of ice–bed contact (Fig. 1a). Making these assumptions, defining the mean stress as p _m = 1/2(σ ₁₁ + σ ₂₂), and using Glen's flow law yields the effective strain rate, $\dot{\varepsilon} _{\rm e}$:

(7)$$\dot{\varepsilon} _{\rm e} = {\rm \;} A\left\{ {\displaystyle{1 \over {\sqrt 2}} {\lsqb {{\lpar {\overline {\sigma_{\rm n}} -p_{\rm m}} \rpar }^2 + {\lpar {\,p_{\rm w}-p_{\rm m}} \rpar }^2} \rsqb }^{1/2}} \right\}^n,$$

where the term in the outermost brackets is the effective shear stress, (1/2σ _ijσ _ij)^1/2 (e.g. Hooke, Reference Hooke2005), and $\overline {\sigma _{\rm n}} \;$ is the normal stress (Eqn (5)) averaged over the zone of ice–bed contact. We idealize this strain rate as uniform in the thin zone of ice adjacent to stoss surfaces where debris particles in contact with the bed reside. Noting that $\dot{\varepsilon} _{\rm e} = (1/\sqrt 2 )\lpar {\dot{\epsilon}_{11}^2 + \; \dot{\epsilon}_{22}^2} \rpar ^{1/2}$ and that for the plane strain of this problem, continuity requires $\dot{\varepsilon} _{11} = -\dot{\varepsilon} _{22}\;$(Fig. 1b), then $\dot{\varepsilon} _{\rm e} = \dot{\varepsilon} _{11}$. Integrating $\dot{\varepsilon} _{11}$ from the bed surface along a coordinate, ν, normal to the bed indicates that v _d = $\dot{\varepsilon} _{11}$ν. The value of v _d that is most relevant to particle–bed contact forces is the velocity at the midpoint of a particle (Hallet, Reference Hallet1981), which corresponds to its radius, R (Fig. 1b). Thus, taking ν = R, and using Eqn (7) yields the component of ice velocity toward the bed due to the bed-parallel extension of ice:

(8)$$v_{\rm d} = \displaystyle{{AR} \over {2^{n/2}}}\lsqb {{\lpar {\overline {\sigma_{\rm n}} -p_{\rm m}} \rpar }^2 + {\lpar {\,p_{\rm w}-p_{\rm m}} \rpar }^2} \rsqb ^{n/2}.$$

The other component of bed-normal ice velocity results from basal melting. The heat flux dissipated by sliding ice is q _s = τu. We assume that all of this heat flux and the geothermal heat flux, q _G, melt ice only along the ice–bed contacts where ice is coldest, so the rate of basal melting, v _m, depends on cavity size through s:

(9)$$v_{\rm m} = \; \displaystyle{1 \over {s\rho L}}\lpar {q_{\rm G} + \tau u} \rpar.$$

As in the derivation of Eqn (3), heat dissipated by water flowing through cavities is neglected. Equations (8) and (9) provide the total bed-normal ice velocity, v _n (Fig. 1b) in zones of ice–bed contact: v _n = v _d + v _m.

Contact force, F, exerted normal to the bed by particles in zones of ice–bed contact, depends on v _n because as ice moves toward the bed it exerts a bed-normal drag on particles (Fig. 1b). Watts (Reference Watts1974) determined the drag on an isolated sphere exerted by temperate ice moving by regelation and enhanced creep. Hallet (Reference Hallet1979, Reference Hallet1981) adapted this result and added to it the sphere's buoyant weight in ice to obtain

(10)$$F = \displaystyle{{\,f4\pi \eta _{\rm e}R^3} \over {R_{\rm {^\ast}}^2 + R^2}}v_{\rm n} + \displaystyle{4 \over 3}\pi R^3g\lpar {\rho_{\rm r}-\rho} \rpar \cos \theta, $$

where η _e is the effective viscosity of ice, R _* is the controlling particle size, f is an empirical factor >1.0 that accounts for the accentuating effect of the bed proximity on drag, ρ _r is the particle density and θ is the angle between the local bed slope and horizontal (Fig. 1b). Particle buoyant weight, as indicated by the right-hand term, is a significant fraction of F only for R > 0.1 m (Hallet, Reference Hallet1979). Experiments with inclusions in ice in which v _n has been measured indicate that drag indeed depends on v _n (Iverson, Reference Iverson1990; Byers and others, Reference Byers, Cohen and Iverson2012), with f ≈ 2.0 (Byers and others, Reference Byers, Cohen and Iverson2012). Cavities that form in these spaces between particles and the bed led Boulton (Reference Boulton and Coates1974) to suggest that the difference between the ice pressure and water pressure in these cavities controlled contact forces. However, water pressure in cavities of steady size beneath isolated particles on the bed must be at least slightly higher than the mean ice pressure to enable drainage from cavities through the water film along the bed surface to ‘connected’ parts of the subglacial hydraulic system. We thus use Eqn (10) to estimate contact forces but acknowledge that it is contingent on the assumption that water produced by melting along particle surfaces drains through the water film at the ice–rock interface. Equally important is that Eqn (10) applies only to debris particles that are sufficiently isolated from each other that their flow fields do not interact.

Owing to the power-law rheology of ice, η _e and R _* in Eqn (10) depend on v _n. The effective viscosity depends on the flow-law parameters:

(11)$$\eta _{\rm e} = \displaystyle{{\dot{\varepsilon} _{r\nu} ^{\lpar {1/n-1} \rpar }} \over {2A^{1/n}}},$$

where r is the distance from particle centers parallel to the local bed slope and remembering that ν is the coordinate perpendicular to the local bed slope (Fig. 1b). Approximating η _e therefore requires choosing a reference value of the shear strain rate, $\dot{\varepsilon} _{r\nu}$. Neglecting bed-normal gradients in velocity parallel to the local bed slope, $\dot{\varepsilon} _{r\nu} = (1/2)(\; \delta v_{\rm n}/\delta r)$. The radial variation in bed-normal velocity, δv _n, is at most v _n. Taking δv _n = v _n, and choosing δr as the particle diameter, 2R, then yields

(12)$$\dot{\varepsilon} _{r\nu} = \displaystyle{{v_{\rm n}} \over {4R}}.$$

The resultant value of effective viscosity from Eqn (11) determines in Eqn (10) the value of R _*:

(13)$$R_{\rm {^\ast}} = \left( {\displaystyle{{3\eta_{\rm e}C_{\rm m}K} \over {\rho L}}} \right)^{1/2},$$

where C _m is the melting point depression with pressure, and K is a bulk thermal conductivity of particles and ice (Watts, Reference Watts1974).

From the contact force, F, and fractional area of ice–bed contact, s, the component of bed shear stress due to debris-bed friction, τ _d, is

(14)$$\tau _{\rm d} = sC\mu F\cos \theta, $$

where C is the concentration of particles in contact with the bed (particle number per unit area), μ is the coefficient friction between rock particles and the bed and cos θ isolates the component of stress parallel to the average bed slope, assumed to be zero. Eqn (14) may overestimate τ _d because particles may rotate in ice under a torque smaller than the torque, μFR, caused by full mobilization of clast-bed friction. We neglect this complication to maximize the possible effect of debris-bed friction on the sliding rule. For the same reason we consider only maximum values of C, limited by the requirement that debris particles be widely enough spaced so that flow fields around them do not interact. Velocity distributions in ice during its power-law flow past an isolated sphere (Lliboutry and Ritz, Reference Lliboutry and Ritz1978) indicate that particles separated by at least ~1.5 particle diameters satisfy this criterion. This spacing of particles at the bed surface implies that the number of particles per unit area of the bed is

(15)$$C = \displaystyle{1 \over {6.25\pi R^2}}.$$

If the same particle spacing applies normal to the bed, then the associated concentration of debris in ice by mass is ~15%.

Equations (3–5) and (8–13) allow estimation of contact forces and areas of ice–bed contact as a function of sliding speed, and Eqn (14) provides the resulting estimate of debris-bed friction, maximized for each particle size through Eqn (15). Table 1 shows relevant parameter values.

Table 1. Parameter values

^a Values from Cuffey and Paterson (Reference Cuffey and Paterson2010) unless noted otherwise here or in the text.

^b Watts (Reference Watts1974).

^c Jaeger and Cook (Reference Jaeger and Cook1979).

^d Pollack and others (Reference Pollack, Hurter and Johnson1993).